VMD-L Mailing List

From: David Poger (d.poger_at_uq.edu.au)
Date: Mon Aug 04 2008 - 20:31:10 CDT

Hello

I have written a Tcl script that basically reads a trajectory and for
each frame, finds the shortest distance between every single residue of
a first domain and every single residue of another domain in a protein.
The final matrix with all the minimal distances is written in an output
file. The script works fine... for 2 frames. In fact, after a few
minutes (~10-20 frames), I have to kill it as it fills up the memory of
the PC (32GB of RAM!!!). Is there any way to optimize this script
(release memory for example?). I know I have several nested loops but I
don't how to do differentlty and I have absolutely no knowledge in code
optimization!

Thanks a lot!
David

# In order to run it start vmd without graphical interface by typing
# $> vmd -dispdev text -e script.tcl

#======================================================#
# TO BE EDITED #

set grofile protein.gro
set trajfile traj.xtc
set outputfile traj_contact.dat

# #
#======================================================#

set datafile [open $outputfile w]
puts $datafile [format "%1s%7s%8s%9s%9s" "#" "Res1" "Res2" "<d>" "rmsd"]

mol new $grofile

# Residues in domains 1D and 2D
set Dom1_1 20
set Dom1_n 155
set Dom2_1 126
set Dom2_n 244

for {set i $Dom1_1 } { $i<=$Dom1_n } {incr i} {
    for {set j $Dom2_1 } { $j<=$Dom2_n } {incr j} {
        set sumdistij($i,$j) 0
        set sumdistij2($i,$j) 0
    }
}

#animate style once
mol addfile $trajfile waitfor all
set nframes [molinfo top get numframes]
puts "nframes=$nframes"

# The 1st frame is the $grofile used for the topology so I don't use it
# That's why k starts from 1 and not 0
for {set k 1} {$k<$nframes} {incr k} {
    animate goto $k
    display update
    puts "frame #$k"
   
    for {set i $Dom1_1 } { $i<=$Dom1_n } {incr i} {
        # Number of atoms in the residues of Dom1 which is selected
        set isel [atomselect top "resid $i"]

        for {set j $Dom2_1 } { $j<=$Dom2_n } {incr j} {
            set distm 1000

            foreach coordi [$isel get {x y z}] {
                # Number of atoms in the residues of Dom2 which is selected
                set jsel [atomselect top "resid $j"]

                foreach coordj [$jsel get {x y z}] {
                    set dist [vecdist $coordi $coordj]
                    if {$dist<$distm} {set distm $dist}
                }
            }
            # The array is in angstroem!!
            set sumdistij($i,$j) [expr $sumdistij($i,$j)+$distm]
            set sumdistij2($i,$j) [expr $sumdistij2($i,$j)+$distm*$distm]

        }
    }

}

for {set i $Dom1_1 } { $i<=$Dom1_n } {incr i} {
    for {set j $Dom2_1 } { $j<=$Dom2_n } {incr j} {
        set mean [expr $sumdistij($i,$j)/($nframes-1)]
        set mean2 [expr $sumdistij2($i,$j)/($nframes-1)]
        set msd [expr $mean2-$mean*$mean]
        set rmsd [expr sqrt($msd)]

        puts $datafile [format "%8d%8d%9.3f%9.3f" $i $j [expr $mean/10]
[expr $rmsd/10]]
    }
}
 
close $datafile

puts "Finished"

exit 

-- 
_________________________________________________
David POGER, Ph.D.
Postdoctoral Research Fellow
Molecular Dynamics Group
School of Molecular and Microbial Sciences (SMMS)
Chemistry Building (#68)
The University of Queensland
Brisbane, QLD 4067
Australia
Email: d.poger_at_uq.edu.au
Tel: +61 7 3365 7562
Fax: +61 7 3365 3872
** Notice: Unless stated otherwise, this e-mail represents only the views of the Sender and not the views of The University of Queensland.