VMD-L Mailing List

From: Axel Kohlmeyer (akohlmey_at_vitae.cmm.upenn.edu)
Date: Thu Feb 02 2006 - 22:03:15 CST

On Thu, 2 Feb 2006, Eliud Oloo wrote:

EO>
EO>
EO> I need to calculate pairwise distances between two sets of residues in
EO> my trajectory. The idea is to determine which two residues come within
EO> a specified cut-off distance relative to each other during the
EO> simulation. The Tcl analysis script below runs well if the two sets are
EO> small (contain less than 10 residues). If the sets are large (~230),
EO> vmd hangs up midway through the calculation. Is there a trick I can
EO> invoke in my script to prevent this from hapenning.?

eliud,

your script has a memory leak. each atomselect command will create
a new selection and thus allocate new memory.

the brute force method would be to delete them after each use, but
you can actually advance a selection by calling: $sel frame $i
or alike and thus just need to define the selections up front and
store them in a list. see below...

EO>
EO> Thanks
EO> Eliud
EO>
EO> # Given two groups eg. two different chains, calculate the distance
EO> #between the geometric centres of all residue pairs and print out pair
EO> #distances that are less than a certain specified cutoff
EO>
EO> #select top molecule
EO> set mol 0
EO>
EO> #set cut-off distance
EO> set cutoff 10
EO>
EO> #specify residue id numbers of first and last residue in each group
EO>
EO> #group one
EO> set minA 1
EO> set maxA 200
EO> #group two
EO> set minB 9853
EO> set maxB 10084
EO>
EO> #open file for writing
EO> set distfile [open pairdistances.dat w]
EO>
EO> #get total number of frames
EO> set nf [ molinfo $mol get numframes ]
EO> #loop over frames
EO> for {set i 0} {$i < $nf} {incr i 100} {

add here:

for {set j $minA} {$j <= $maxA} {incr j} {
set selA($j) = [atomselect $mol "resid $i"]
}

for {set j $minB} {$j <= $maxB} {incr j} {
set selB($j) = [atomselect $mol "resid $i"]
}

EO> #loop over group one
EO> for {set j $minA} {$j<=$maxA} {incr j} {

and replace:

EO> set group_a [atomselect $mol "resid $j" frame $i]
EO> set com_a [measure center $group_a]

by:

$selA($j) frame $i
set com_a [measure center $selA($j)]

EO> #loop over group two
EO> for {set k $minB} {$k<=$maxB} {incr k} {

and again:

EO> set group_b [atomselect $mol "resid $k" frame $i]
EO> set com_b [measure center $group_b]

$selB($k) frame $i
set com_b [measure center $selB($k)]

that should do the trick.

regards,
     axel.

EO> #calculate distance and check if value is within cutoff
EO> set dist [veclength [vecsub $com_a $com_b]]
EO> if {[expr $dist <= $cutoff]} {
EO> # puts "frame $i respairs $j $k dist $dist
EO> centers $com_a $com_b"
EO> puts $distfile "$i $j $k $dist"
EO>
EO> }
EO>
EO> }
EO> }
EO> }
EO> close $distfile
EO>
EO>
EO> 

-- 
=======================================================================
Axel Kohlmeyer   akohlmey_at_cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.