From: Charles McAnany (
Date: Wed Apr 01 2015 - 14:42:37 CDT

I'm analyzing a set of trajectories based on nearest-atom contact maps: for
each pair of residues (the x and y axes), what is the distance between the
closest atoms on those residues? (hydrogens excluded, of course) I'm doing
this because I'm interested in the side-chain interactions, so just looking
at CA distances might miss something cool.

My contact map has two parts: one part shows the closest two residues get
during the entire trajectory, essentially representing what parts of the
protein can interact. This is easy to calculate quickly by using 'measure
bond' and specifying the entire trajectory for the frame range, then
getting the minimum value of the returned list.

The other part shows the average closest-atom distance over the trajectory,
representing the parts of the protein that spend lots of time interacting.

For the life of me, I can't figure out how to efficiently write the second
analysis. The awkward part is that the outer loop needs to be the frame
number, since I'm averaging a property of each frame, not of each atom.
This means lots and lots of calls to 'measure bond' and painfully long
analysis times.

Here's the problem area of what I have so far:
for {set frameNumber 0} {$frameNumber < $numFrames} {incr frameNumber 1} {
    #Just a very large number, this will go away when compared with a real
    set closestAtomsThisFrame 1E30
    foreach atom1 [$residue1 list] {
        foreach atom2 [$residue2 list] {
            set thisDist [measure bond [list $atom1 $atom2 ] frame
            if {$thisDist < $closestAtomsThisFrame} {
                 set closestAtomsThisFrame $thisDist
     set averageContact [expr { $averageContact + $closestAtomsThisFrame /
$numFrames }]
     if {$closestAtomsThisFrame < $minContact} {
         set minContact $closestAtomsThisFrame

Does anyone have advice on getting tcl to perform better on this sort of