Example Protein in minimum water sphere

The script(min_wat.tcl) for solvating protein is given below:

# finds a center of mass of the molecure place a shpere of water around it.
# from NAMD 2003 summer tutorials

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         myo.psf         
set pdb         myo.pdb
set box         myo_box
set psfDrop     myo_ws.psf   ;# psf file of solvated protein
set pdbDrop     myo_ws.pdb   ;# pdb file of solvated protein

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 mass 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"

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"
mol delete top
mol load psf ${psfDrop} pdb ${pdbDrop}

First we load the psf and pdb files to vmd, then source the tcl script.

Load the psf file to vmd using the selection File> New Molecule, then pdb file using File $>$ Load Data Into Molecule in VMD Main window.Source the script by using selection Extensions $>$ tkcon inVMD Main window

source min_wat.tcl

The above script generates the water sphere of radius of 27.0$\AA$ centered at center of mass(2.824,3.829,10.110)of the protein.

Here is the snap of solvated apomyoglobin.