VMD-L Mailing List

From: politr_at_huji.ac.il
Date: Thu Apr 24 2008 - 03:31:22 CDT

Dear all,
I'm trying to run my script. This script reads a trajectory and rotates
through all the frames. For each frame it perform the following: saves
all the water atoms from each frame in some variable (total_wat), after
that, for each OH2 atom I find the 4 closest OH2 atoms and measure the
angle between them. I don't want to use atoms on the boundaries of my
box so I have defined some max radius (max_length) and I use only the
water molecules within this radius. The problem is that the script works
well only for 25 frames and after that it stuck the computer. I have
tried to monitor the script performance with top and I have discovered
that the script is using to much memory and at some point I'm out of the
memory. My script is attached. I will appreciate it very much if anyone
could help. Any help is more than welcome.
Thank you in advance.
Regina.

proc angle {D H A} {
# cosine of the angle between three points
# get Pi
# global M_PI
   set M_PI 3.1415927

## setup vectors hd and ha
set hd [vecsub $D $H]
set ha [vecsub $A $H]
set cosine [expr [vecdot $hd $ha] / ( [veclength $hd] * [veclength $ha])]

# convert cosine to angle in degrees
#return [expr acos($cosine)*(180.0/$M_PI)]
return $cosine
}

set q_name [lindex $argv 2]

set q_file [open /usr/people/daniel/politr/namd/peptide/analysis/$q_name w]

set psf [lindex $argv 0]
set dcd [lindex $argv 1]
set psf_filename $psf
set dcd_filename $dcd
mol load psf $psf_filename
animate read dcd $dcd_filename waitfor all
set n [molinfo top get numframes]
set all [atomselect top all]

set total_wat [atomselect top "resname TIP3 and name OH2 "]
 
#rotate through all the frames
for { set i 0 } { $i < $n } { incr i 10} {

      $all frame $i
     $all update
       set mol_ID [molinfo top]
       molinfo $mol_ID set frame $i
       set box [molinfo $mol_ID get {a b c}]
       set length [lindex $box 0]

      set max_length [expr $length/2.0-3.0]

       set all [atomselect top all]
       set center [measure center $all]
       set cen_x [lindex $center 0]
       set cen_y [lindex $center 1]
       set cen_z [lindex $center 2]
       $total_wat frame $i
       $total_wat update
    #to create list of resids. after that we can go through all coord
and through all resids i
n foreach loop (if we need to print resid)
       set wat_resid [$total_wat get resid]

   foreach water $wat_resid {
        set wat_coord [atomselect top "name OH2 and resid $water"]
       set xyz [$wat_coord get {x y z}]
        set wat_x [lindex [lindex $xyz 0] 0]
        set wat_y [lindex [lindex $xyz 0] 1]
        set wat_z [lindex [lindex $xyz 0] 2]
       set x [expr $wat_x-$cen_x]
       set y [expr $wat_y-$cen_y]
       set z [expr $wat_z-$cen_z]
       set distance [expr { sqrt( $x * $x + $y * $y + $z*$z) }]
       if {$distance <= $max_length} {
       #this part is used to find 4 closest water molecules
       set max 6
       set min [expr $max/2.0]
       set smalest 0
       set close_wat [atomselect top "name OH2 and within $min of (resid
$water and name OH2)"]
       set ntot [$close_wat num]
       while {$ntot != 5} {

           if {$ntot > 5} {
               set max $min
               set min [expr $smalest + ($max-$smalest)/2.0]
      }
           if {$ntot < 5} {
               set smalest $min
               set min [expr $min + ($max-$min)/2.0]
           }
           set close_wat [atomselect top "name OH2 and within $min of
(resid $water and name OH2)
"]
           set ntot [$close_wat num]
       }
#find the center atom in list of atoms

       for { set j 0 } { $j < 5 } { incr j } {

           set resid [lindex [$close_wat get resid] $j]
           if {$resid == $water} {
               set c_atm $j
           }

       }
       set sum 0.0
       for { set k 0 } { $k < 4 } { incr k } {
           for { set l [expr $k+1] } { $l < 5 } { incr l } {
           if {$k != $c_atm && $l != $c_atm} {
               set atom1 [lindex [$close_wat get {x y z}] $k]
               set atom2 [lindex [$close_wat get {x y z}] $c_atm]
               set atom3 [lindex [$close_wat get {x y z}] $l]
               set angle [angle $atom1 $atom2 $atom3]
               set sum [expr {(($angle+1.0/3.0)*($angle+1.0/3.0))+$sum}]

           }
           }
       }

       set q [expr {1.0-((3.0/8.0)*$sum)}]
       puts $q_file "$q"

}

   }
   puts $q_file "ENDMDL"
}
close $q_file