VMD-L Mailing List

From: John Stone (johns_at_ks.uiuc.edu)
Date: Thu Apr 24 2008 - 09:20:17 CDT

Hi,
  Your script has an atom selection leak in it. You are
setting the "all" atom selection in your loop, but you never
call "$all delete", so you're leaking an atom selection on every frame
of your loop. There may be other problems I didn't notice, but that's
an obvious one.

Cheers,
  John Stone
  vmd_at_ks.uiuc.edu

On Thu, Apr 24, 2008 at 11:30:19AM +0300, Regina Politi wrote:
> Dear VMD users and developers,
> 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 stucks 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
>
>
> 

-- 
NIH Resource for Macromolecular Modeling and Bioinformatics
Beckman Institute for Advanced Science and Technology
University of Illinois, 405 N. Mathews Ave, Urbana, IL 61801
Email: johns_at_ks.uiuc.edu                 Phone: 217-244-3349
  WWW: http://www.ks.uiuc.edu/~johns/      Fax: 217-244-6078