From: Guo Zhi (gzgzgz_at_gmail.com)
Date: Sun Apr 15 2007 - 00:54:48 CDT

Hi all, I am trying to use a script to do statistics on sidechain
contact in the whole trajectory. Since the MD trajectory were yielded
by a home-made MD program which uses only one united atom for each
sidechain, the selection code might looks a little strange. Well, the
script seem to produce the right result, but the memory usage
coninuously grow up which implies certain memory leakage. I check the
code for several times(mainly the "atomselect" sentences) and still
can not find any where that cause this problem. Would anyone of you
bother to look into this and give me some help?

########### Code for sidechain contact statistics ########

set frame_no [molinfo top get numframes]
set res_set [lsort -increasing {"ASN" "HIS" "VAL" "THR" "LEU" "SER" "GLN"}]
set file_id [open "contact.dat" w]

set sel_all [atomselect top all]
set res_leng [llength [lsort -unique [$sel_all get resid]]]

for {set m 0} {$m<$frame_no} {incr m} {
puts $file_id "Frame $m:"
set all_contact {}
for {set i 0} {$i<$res_leng} {incr i} {
        set tmp_sel [atomselect top "(sidechain or (not protein)) and
(exwithin 6.5 of ((resid $i) and (sidechain or (not protein))))" frame
$m]
        
        set contact_res [$tmp_sel get resid]
        foreach ele $contact_res {
                set ele_contact $i
                lappend ele_contact $ele
                set ele_contact [lsort -increasing $ele_contact]
                lappend all_contact $ele_contact
        }

        $tmp_sel delete
}
        
        set short_contact [lsort -unique $all_contact]
        foreach ele_short $short_contact {
                set ele_name [atomselect top "(resid $ele_short) and (name CA)"]
                set pos [lsearch -exact $short_contact $ele_short]
                set short_contact [lreplace $short_contact $pos $pos [lsort
-increasing [$ele_name get resname]]]
                $ele_name delete
        }
        
        # initialize matrix
        
        set matrix {}
        
        for {set j 0} {$j<28} {incr j} {
                lappend matrix 0
        }
        
        foreach ele_short $short_contact {
                set count 0
                for {set j 0} {$j<7} {incr j} {
                        for {set k $j} {$k<7} {incr k} {
                                if {[lindex $ele_short 0]==[lindex $res_set $j]} {
                                        if {[lindex $ele_short 1]==[lindex $res_set $k]} {
                                                # Now set matrix element
                                                set new_ele [expr [lindex $matrix $count]+1]
                                                set matrix [lreplace $matrix $count $count $new_ele]
                                        }
                                }
                                incr count
                        }
                }
        }
        
        puts $file_id $matrix
}

$sel_all delete
close $file_id