VMD-L Mailing List

From: Alberto Sergio Garay (sgaray_at_fbcb.unl.edu.ar)
Date: Wed Jun 03 2009 - 12:54:53 CDT

Dear Axel Kohlmeyer

I have checked my script following your recommendation:

> you may have mixed up loop nesting!
> memory leaks in this kind of scripting almost
> always (= 99%) are due to forgotten or
> overwritten(!) selections.

But I could not find any mistake so I'm pasting a part of the script
which when it is run takes some of the ram memory and never leaves it until
I close vmd.

BTW, I'm waiting for my password in order to be able to up my md
gromacs trayectory files(gro and xtc ) for testing the script.

Thank you in advance

#Script to look for some molecules around (radius) a particular
residue (resid) in all the
#periodic replicas of the principal box of simulation and move them
close to it.
#The script is thought to work on lipids bilayers.

proc rsel {resid radius} {
set num_steps [molinfo top get numframes]

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} ]

#Calculating the min-max dimention of the selected molecule
   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
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]}

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]

#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

#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

}
puts "End of calculation..."
} 

-- 
Dr. Sergio Garay
Facultad de Bioquimica y Cs. Biológicas
Universidad Nacional del Litoral
Santa Fe - Argentina
C.C. 242 - Ciudad Universitaria - C.P. S3000ZAA
Argentina
Ph. +54 (342) 4575-213
Fax. +54 (342) 4575-221