From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Fri Jan 27 2012 - 17:24:23 CST

On Fri, Jan 27, 2012 at 5:26 PM, Mustafa Tekpinar <tekpinar_at_buffalo.edu> wrote:
> Hi,
> I have a molecular dynamics trajectory of (protein+ligand+water+ions) box. I
> want to cut a sphere whose center is center of mass of (protein+ligand).
> After determining the vector from center of mass to the farthest atom, I
> want add some buffer distance to that radius. Then, I want to save all atoms
> (except ions) in that radius. I have a script but it seems like there is a
> memory issue. I'd be glad if some experts points me my mistake(s).

most likely you are leaking atom selections (like about 95% of the
people that have reported memory problems in their scripts).

>
> Here is my script:
> //=====================================================================
> set hydra_shell 7.0
> set frm_beg 0
> set frm_end 5
>
> #Load protein box simulation
>
>
> mol load psf MD_protein_in_water/build/solvated.psf dcd
> MD_protein_in_water/prod/equil.dcd

this is not a bug, but "mol load" is a deprecated API and
will go away in the future. you should use "mol new" and
"mol addfile" instead.

> #Select protein and (if exists) ligands
>
>
> set protAndLigand [atomselect top "not water and not ions"]
>
> for {set frame $frm_beg} {$frame < $frm_end} {incr frame} {
>
>     # get the correct frame
>
>
>     $protAndLigand frame $frame
>
>     ## Determine the center of mass of the (protein+ligands) and store the
> coordinates
>
>     set CoM [measure center $protAndLigand weight mass]
>     set x1 [lindex $CoM 0]
>     set y1 [lindex $CoM 1]
>     set z1 [lindex $CoM 2]
>
>     ## Determine the distance of the farthest atom from the center of mass
>
>
>     foreach atom [ $protAndLigand get index] {
>         set pos [lindex [[atomselect top "index $atom" frame $frame ] get {x
> y z}] 0]

yup. bingo! here you create atom selection and don't delete it.
furthermore, you also create atom selections (which are costly)
without need. the following code does the same thing:

foreach pos [$protAndLigang get {x y z}] {

>         set x2 [lindex $pos 0]
>         set y2 [lindex $pos 1]
>         set z2 [lindex $pos 2]
>
>         set dist [expr (($x2-$x1)*($x2-$x1) + ($y2-$y1)*($y2-$y1) +
> ($z2-$z1)*($z2-$z1))]
>         if {$dist > $max} {
>             set max $dist
>         }
>     }
>
>    set max [expr sqrt($max)]
>
>    #Add some buffer of water to max distance
>    set max [expr $max + $hydra_shell]
>    set maxSqrd [expr $max*$max]
>    puts "Radius squared is $maxSqrd]"
>
>    #Select protein and (if exists) ligand +plus water within the radius!
>
>
>    set pwSphr [atomselect top "(not water and not ions) or (same residue as
> {water and ((x-$x1)*(x-$x1) + (y-$y1)*(y-$y1) +
> (z-$z1)*(z-$z1))<($maxSqrd)})" frame $frame]
>
>    #write coordinates.
>
>
>    $pwSphr writepdb sphere${frame}.pdb
>
>    $pwSphr delete
> }
> $protAndLigand delete
> exit
>
> //=====================================================================
> Thanks in advance,
>
> mtekpinar

that should fix it.

axel.

-- 
Dr. Axel Kohlmeyer
akohlmey_at_gmail.com  http://goo.gl/1wk0
College of Science and Technology
Temple University, Philadelphia PA, USA.