From: Eliud Oloo (elo712_at_mail.usask.ca)
Date: Thu Feb 02 2006 - 19:26:35 CST

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

Thanks
Eliud

# Given two groups eg. two different chains, calculate the distance
#between the geometric centres of all residue pairs and print out pair
#distances that are less than a certain specified cutoff

#select top molecule
set mol 0

#set cut-off distance
set cutoff 10

#specify residue id numbers of first and last residue in each group

#group one
set minA 1
set maxA 200
#group two
set minB 9853
set maxB 10084

#open file for writing
set distfile [open pairdistances.dat w]

#get total number of frames
set nf [ molinfo $mol get numframes ]
#loop over frames
for {set i 0} {$i < $nf} {incr i 100} {
#loop over group one
        for {set j $minA} {$j<=$maxA} {incr j} {
        set group_a [atomselect $mol "resid $j" frame $i]
        set com_a [measure center $group_a]
        #loop over group two
                for {set k $minB} {$k<=$maxB} {incr k} {
                set group_b [atomselect $mol "resid $k" frame $i]
                set com_b [measure center $group_b]
                #calculate distance and check if value is within cutoff
                set dist [veclength [vecsub $com_a $com_b]]
                    if {[expr $dist <= $cutoff]} {
# puts "frame $i respairs $j $k dist $dist
centers $com_a $com_b"
                        puts $distfile "$i $j $k $dist"

                    }

                }
        }
}
close $distfile