VMD-L Mailing List

From: Mustafa Tekpinar (tekpinar_at_buffalo.edu)
Date: Fri Jan 27 2012 - 16:26:29 CST

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).

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

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