VMD-L Mailing List

From: Justin Gullingsrud (justin_at_ks.uiuc.edu)
Date: Sun Aug 19 2001 - 18:48:02 CDT

Hi,

On Wed, Aug 15, 2001 at 02:03:11PM -0600, Elizabeth K. White wrote:
> I modified the RMSD calculator to align parts of a protein over multiple
> frames, in order to see what regions might be folding. The calculations
> are going suspiciously slowly, and VMD often crashes without finishing
> them. I am writing the files to a huge disk, so space shouldn't be a
> problem. Can anyone suggest what is wrong? I am running VMD 1.6 on a DEC
> Alpha using Linux. The script follows.

What's probably happening is that you're running out of memory because the
atom selections in the inner loop of your script allocate resources that
aren't freed until Tcl gets out of the for loop. Try running top in a
separate window while you run your script to be sure.

Here are two ways you can modify your script to keep this from happening:
        - Add the following lines to the end of your inner loop:
                unset upproc_var_$reference
                unset upproc_var_$compare

        or:

        - Put the inner loop of your script inside a proc and call it from
          the inner loop.

In both cases, the atom selections will be deleted each time through the
inner loop, so you won't run out of system memory. The first method will
definitely be faster, while the second method is more robust in that it
doesn't mess with VMD internals that could possibly change in a future
revision of VMD.

Incidentally, you could probably speed up your script quite a bit by
reversing the order of the loops, i.e. move the the line where you set up
the loop over frames to just after the line where you create the compare
selection.
 
Hope this helps,

Justin

>
> Thanks, Elizabeth
>
> # Prints the RMSD of the protein atoms between each timestep
> # and the first timestep for the given molecule id (default: top)
> proc overlap_rmsd {aminos dir seglength {mol top}} {
> puts "$aminos aminos, $dir subdir, $seglength seglength"
> set file1 [molinfo top get {source filename}]
> set file2 [molinfo top get {source filename2}]
> # note the filename of the molecule
> exec echo ".parm file is $file1" > thisframe.out
> exec cat thisframe.out >> header.out
> exec echo ".crd file is $file2" > thisframe.out
> exec cat thisframe.out >> header.out
> set num_steps [molinfo $mol get numframes]
> set maxk [expr $aminos - $seglength + 2]
> for {set frame 0} {$frame < $num_steps} {incr frame} {
> for {set k 1} {$k < $maxk} {incr k} {
> set m [expr $k + $seglength - 1]
> # use frame 0 for the reference
> set reference [atomselect $mol "protein and resid
> $k to $m" frame 0]
> # the frame being compared
> set compare [atomselect $mol "protein and resid $k
> to $m"]
> # get the correct frame
> $compare frame $frame
> # compute the transformation
> set trans_mat [measure fit $compare $reference]
> # do the alignment
> $compare move $trans_mat
> # compute the RMSD
> set rmsd [measure rmsd $compare $reference]
> set spc "_"
> set output "$seglength$spc$frame"
> # puts "output is $output"
> # print the RMSD
> puts "frame=$frame k=$k m=$m $rmsd=rmsd"
> set x [expr ($k + $m) / 2.0]
> exec echo "$x $rmsd" > thisframe.out
> exec cat thisframe.out >>
> /home/ri/ekwhite/rmsd/$dir/$output
> }
> }
> puts "Okay, done"
> }
> 

-- 
Justin Gullingsrud      3111 Beckman Institute
H: (217) 384-4220       I got a million ideas that I ain't even rocked yet...
W: (217) 244-8946       -- Mike D