Water Sphere tcl Script

# finds a center of mass of the molecure (ubquitin), place a shpere of water 
# around it.

# To run execute: vmd -dispdev text  -e waterSphere.tcl


proc center_of_mass {selection} {

        # some error checking
        if {[$selection num] <= 0} {
                error "center_of_mass: needs a selection with atoms"
        }
        # set the center of mass to 0
        set com [veczero]
        # set the total mass to 0
        set mass 0
        # [$selection get {x y z}] returns the coordinates {x y z} 
        # [$selection get {mass}] returns the masses
        # so the following says "for each pair of {coordinates} and masses,
        #  do the computation ..."
        foreach coord [$selection get {x y z}] m [$selection get mass] {
           # sum of the masses
           set mass [expr $mass + $m]
           # sum up the product of mass and coordinate
           set com [vecadd $com [vecscale $m $coord]]
        }
        # and scale by the inverse of the number of atoms
        if {$mass == 0} {
                error "center_of_mass: total mass is zero"
        }
        # The "1.0" can't be "1", since otherwise integer division is done
        return [vecscale [expr 1.0/$mass] $com]
}

################################################################
#                 MAIN PART STARTS HERE
################################################################


set psf         ubq.psf
set pdb         ubq.pdb
set box         ubq_box
set psfDrop     ubq_ws.psf
set pdbDrop     ubq_ws.pdb

package require psfgen

resetpsf

mol load psf $psf pdb $pdb
set sel [atomselect top all]
# find mass center
set center [center_of_mass $sel]
puts "center of mas is at $center"


foreach {xmass ymass zmass} $center { break }


set num0 9999
set Rmin 0.0

while {$num0 != 0} {
    set Rmin [expr $Rmin +1.0]
    set probSel [atomselect top "not (sqr(x-$xmass) + sqr(y-$ymass) +
sqr(z-$zmass) <= sqr($Rmin))"]
    set num0 [$probSel num]
    puts "$num0 $Rmin"
}


package require solvate

set xmin [expr $xmass -$Rmin]
set xmax [expr $xmass +$Rmin]

set ymin [expr $ymass -$Rmin]
set ymax [expr $ymass +$Rmin]

set zmin [expr $zmass -$Rmin]
set zmax [expr $zmass +$Rmin]

puts " $xmin $ymin $zmin $xmax $ymax $zmax"
#solvate $psf $pdb -o $box -minmax {{15.4399995804 12.8879995346 
-0.365999996662} {46.125 45.2680015564 36.2509994507} }

set min "$xmin $ymin $zmin"
set max "$xmax $ymax $zmax"
set minmax [list $min $max]

solvate $psf $pdb -o $box -minmax $minmax

mol delete top

resetpsf

mol load psf ${box}.psf pdb ${box}.pdb

readpsf  ${box}.psf
coordpdb ${box}.pdb

set selDel [atomselect top "not  (sqr(x-$xmass) + sqr(y-$ymass) + sqr(z-$zmass)
<= sqr($Rmin))" ]
puts  " not within [$selDel num]"

set testSel [atomselect top "not (sqr(x-$xmass) + sqr(y-$ymass) + 
sqr(z-$zmass)  <= sqr($Rmin)) and (not water)"]
puts " not within and not water [$testSel num]"

if { [$testSel num] != 0} { 
    puts "ERROR: there are [$testSel num] non water molecules outside the shell"
    puts "EXIT"
    exit
}



set delList [$selDel get {segid resid}]

set delList [lsort -unique $delList]

foreach record $delList {
    foreach {segid resid} $record { break }
    delatom $segid $resid
}

writepsf $psfDrop
writepdb $pdbDrop

# remove temprorary files generated by the script 
file delete ${box}.psf ${box}.pdb combine.pdb combine.psf

puts "CENTER OF MASS IS AT: $center"
puts "SPHERE RADIUS:        $Rmin"

exit