From: Axel Kohlmeyer (akohlmey_at_cmm.chem.upenn.edu)
Date: Wed Jun 03 2009 - 13:57:40 CDT

On Wed, 2009-06-03 at 14:54 -0300, Alberto Sergio Garay wrote:
> Dear Axel Kohlmeyer

dear alberto,

> I have checked my script following your recommendation:
>

ok. main loop starts here.

> for {set n 0} {$n <= $num_steps} {incr n 1} {
>
> puts "calculation on frame $n"
>
> #Taking the box dimension
> set Xbox [lindex [pbc get] 0 0]
> set Ybox [lindex [pbc get] 0 1]
> set Zbox [lindex [pbc get] 0 2]
>
> #Setting the coordinates of the target residue
> set sx [[atomselect top "resid $resid" frame $n] get {x} ]
> set sy [[atomselect top "resid $resid" frame $n] get {y} ]
> set sz [[atomselect top "resid $resid" frame $n] get {z} ]

... and here is already a memory leak.

you create selections without deleting them due to the nesting
of the square brackets.

> #Calculating the min-max dimention of the selected molecule

have you seen that VMD has a 'measure minmax' command?

> set minx [lindex $sx 0]
> set miny [lindex $sy 0]
> set minz [lindex $sz 0]
> set maxx $minx
> set maxy $miny
> set maxz $minz
> foreach x $sx y $sy z $sz {
> if {$x < $minx} {set minx $x} else {if {$x > $maxx} {set maxx $x}}
> if {$y < $miny} {set miny $y} else {if {$y > $maxy} {set maxy $y}}
> if {$z < $minz} {set minz $z} else {if {$z > $maxz} {set maxz $z}}
> }
>
> #Calculation of geometric centre

there also is 'measure center'...

> set centre_X [expr ($minx + $maxx) / 2]
> set centre_Y [expr ($miny + $maxy) / 2]
> set centre_Z [expr ($minz + $maxz) / 2]
> set cg "$centre_X $centre_Y $centre_Z"
>
> #Moving the geometric centre to each periodic box and###
> #Making a list of water molecules inside the defined sphere
> for {set i 1} {$i <= 3} {incr i} {
> if {$i == 1} {set X $centre_X}
> if {$i == 2} {set X [expr $centre_X - $Xbox]}
> if {$i == 3} {set X [expr $centre_X + $Xbox]}
> for {set j 1} {$j <= 3} {incr j} {
> if {$j == 1} {set Y $centre_Y}
> if {$j == 2} {set Y [expr $centre_Y - $Ybox]}
> if {$j == 3} {set Y [expr $centre_Y + $Ybox]}
> for {set k 1} {$k <= 3} {incr k} {
> if {$k == 1} {set Z $centre_Z}
> if {$k == 2} {set Z [expr $centre_Z - $Zbox]}
> if {$k == 3} {set Z [expr $centre_Z + $Zbox]}

please have a closer look at the VMD-1.8.7 beta documentation.
there is a new 'pbwithin' selection flag. this may help to
reduce this monster construct in combination with the updated
PBCtools to just a few lines...

> set water [atomselect top "same residue as (resname SOL) and
> ((x-($X))*(x-($X)) + (y-($Y))*(y-($Y)) + (z-($Z))*(z-($Z))) <
> ($radius*$radius)" frame $n]
> set lipids [atomselect top "same residue as (resname DMP) and
> ((x-($X))*(x-($X)) + (y-($Y))*(y-($Y)) + (z-($Z))*(z-($Z))) <
> ($radius*$radius)" frame $n]

hmmmm.... creating selection deep inside nested loops.

>
> #Moving the selected molecules close to the target
>
> if {$i == 2} {$water moveby "$Xbox 0 0"
> $lipids moveby "$Xbox 0 0"}
> if {$i == 3} {$water moveby "-$Xbox 0 0"
> $lipids moveby "-$Xbox 0 0"}
> if {$j == 2} {$water moveby "0 $Ybox 0"
> $lipids moveby "0 $Ybox 0"}
> if {$j == 3} {$water moveby "0 -$Ybox 0"
> $lipids moveby "0 -$Ybox 0"}
> if {$k == 2} {$water moveby "0 0 $Zbox"
> $lipids moveby "0 0 $Zbox"}
> if {$k == 3} {$water moveby "0 0 -$Zbox"
> $lipids moveby "0 0 -$Zbox"}
>
> set Nwater [llength [lsort -integer -unique [$water get residue]]]
> set Nlipids [llength [lsort -integer -unique [$lipids get residue]]]
> if {$Nwater > 0} { puts "water molec found: $Nwater" }
> if {$Nlipids > 0} { puts "lipids molec found: $Nlipids" }
>
> }
> }
> }
>
> $water delete
> $lipids delete

...and deleting them _outside_ the loop.
this is the next memory leak. you keep overwriting
the variable with new selection function names and
you delete only the last one.

>
> #Making a new selection around the selected residue
> set water [atomselect top "same residue as (resname SOL) and
> ((x-($centre_X))*(x-($centre_X)) + (y-($centre_Y))*(y-($centre_Y)) +
> (z-($centre_Z))*(z-($centre_Z))) < ($radius*$radius)" frame $n]
> set lipids [atomselect top "same residue as (resname DMP) and
> ((x-($centre_X))*(x-($centre_X)) + (y-($centre_Y))*(y-($centre_Y)) +
> (z-($centre_Z))*(z-($centre_Z))) < ($radius*$radius)" frame $n]
>
> #Here should be start the Hbonds calculation
>
> $water delete
> $lipids delete

this is proper. so no leakage.

>
> }
> puts "End of calculation..."
> }

there you go. you _do_ leak atom selections
and you have to fix your script.

cheers,
   axel.

-- 
=======================================================================
Axel Kohlmeyer   akohlmey_at_cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.