• ## Outreach

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} {

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.
```