VMD-L Mailing List

From: Leonardo Sepulveda Durán (leonardosepulveda_at_gmail.com)
Date: Mon Jan 08 2007 - 09:05:01 CST

Hello. I have 3 systems simutated at different temperatures. For each
System-Temperature pair I have 20 simulations started from different random
seeds. I would like to calculate, for each 20-simulation group, the rmsds
between all trajectory pairs. I wrote a script to do this:

# cycle ovel all simulation temperatures
foreach temp { 300 400 450 500} {

   # Cycle over all systems
   foreach name { paa pac pan} {

      # load the 20 trajectories for each temperature and type of system
      for {set i 1} {$i <= 20} {incr i} {
         mol load psf ./$temp/din-$temp-$name-$i.psf dcd
./$temp/din-$temp-$name-$i.dcd
      }

      # Obtain the index of the top molecule (the last loaded)
      set molnum [molinfo top]

      # define a counter for the traj-traj comparisons
      set counter 0

      # Main calculation cycle. For each temperature and system, this code
calcule all the trajectory-
      # trajectory rmsds, and then print the information to independent
files. The "for" loops decrease their
      # indexes as the provided trajectory number is the last one
      for {set i $molnum} {$i >= [ expr $molnum - 20 + 1 ] } {incr i -1} {
         for {set j [ expr $i - 1] } {$j >= [ expr $molnum - 20 + 1 ] }
{incr j -1} {

            incr counter

            set num_frames [molinfo top get numframes]

            # make selection with all the atoms
            set all [atomselect $i all]

            # Calculate the rmsd between frames k of 2 trajectories
            for {set k 0 } {$k < $num_frames} {incr k } {

               # select the structures to be compared
               set sim1 [atomselect $i "protein" frame $k]
               set sim2 [atomselect $j "protein" frame $k]

               # do the alignment
               set trans_mat [measure fit $sim1 $sim2]
               $all move $trans_mat

               # measure RMSD between frames
               set rms [measure rmsd $sim1 $sim2]

               # define rmsd array numbering
               set h [expr $k + 1]

               set rmsd($h) $rms
            }

            puts "\t$counter $i $j"

            # open output file
            set outfile [open ./rmsd/$name-$temp-$counter.dat w]
            for {set m 1} {$m <= $num_frames} {incr m} {
                  puts $outfile [format "%5d%7.3f" $m $rmsd($m)]
            }
            close $outfile
         }
      }

      unset sim1
      unset sim2
      unset all
      unset rmsd
      unset rms
      unset h

      # delete all the loaded trajectories
      mol delete all

   }
}

I put the script above in a file called rmsdv1.0.tcl, and ran the script
with the command "vmd -dispdev text -e rmsdv1.0.tcl".

It works without any trouble when I use 5 simulations instead of 20. In the
later case, the calculation dies out without any explanation besides writing
"killed" in the console. During the calculation, each iteration becomes
slower, so I think that the problem is asociated to memory usage. I dont
know how the memory is used during the computation, but is reasonable that
all the iterations where calculated at the same rate. It seems that after a
pairwise comparison the compared simulations remain in memory, I tried to
solve the problem using the "unset" statements at the end of the script, but
it seems to have no effect at all. Maybe the problem is another, any help
from a tcl programmer could help me.

Thanks in advance.

Leonardo Sepulveda
Universidad de Chile