set sel [atomselect top "water and same residue as (within 2 of protein)"] set n [molinfo top get numframes] for { set i 0 } { $i < $n } { incr i } { $sel frame $i $sel update $sel writepdb water_$i.pdb }
The frame option sets the frame of the selection, update tells the atom selection to recompute which waters are near the protein, and writepdb writes the selected waters to a file.
proc total_mass {selection} { set sum 0 foreach mass [$selection get mass] { set sum [expr $sum + $mass] } return $sum }Here's another (slower) way to do the same thing. This works because the mass returned from the selection is a list of lists. Putting it inside the quotes of the eval makes it a sequence of vectors, so the vecadd command will work on it.
proc total_mass1 {selection} { set mass [$selection get mass] eval "vecadd $mass" }
proc minmax {molid} { set sel [atomselect $molid all] set coords [$sel get {x y z}] set coord [lvarpop coords] lassign $coord minx miny minz lassign $coord maxx maxy maxz foreach coord $coords { lassign $coord x y z 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}} } return [list [list $minx $miny $minz] [list $maxx $maxy $maxz]] }
proc gyr_radius {sel} { # make sure this is a proper selection and has atoms if {[$sel num] <= 0} { error "gyr_radius: must have at least one atom in selection" } # gyration is sqrt( sum((r(i) - r(center_of_mass))^2) / N) set com [center_of_mass $sel] set sum 0 foreach coord [$sel get {x y z}] { set sum [vecadd $sum [veclength2 [vecsub $coord $com]]] } return [expr sqrt($sum / ([$sel num] + 0.0))] }
Applying this to the alanin.pdb coordinate file
vmd > mol load pdb alanin.pdb vmd > set sel [atomselect top all] vmd > gyr_radius $sel Info) 5.45443
proc frame_rmsd {selection frame1 frame2} { set mol [$selection molindex] # check the range set num [molinfo $mol get numframes] if {$frame1 < 0 || $frame1 >= $num || $frame2 < 0 || $frame2 >= $num} { error "frame_rmsd: frame number out of range" } # get the first coordinate set set sel1 [atomselect $mol [$selection text] frame $frame1] set coords1 [$sel1 get {x y z}] # get the second coordinate set set sel2 [atomselect $mol [$selection text] frame $frame2] set coords2 [$sel2 get {x y z}] # and compute the rmsd values set rmsd 0 foreach coord1 $coords1 coord2 $coords2 { set rmsd [expr $rmsd + [veclength2 [vecsub $coord2 $coord1]]] } # divide by the number of atoms and return the result return [expr $rmsd / ([$selection num] + 0.0)] }
The following uses the frame_rmsd function to list the rmsd of the molecule over the whole trajectory, as compared to the first frame.
vmd > mol load psf alanin.psf dcd alanin.dcd vmd > set sel [atomselect top all] vmd > for {set i 0} {$i < [molinfo top get numframes]} {incr i} { ? puts [list $i [frame_rmsd $sel $i 0]] ? } 0 0.0 1 0.100078 2 0.291405 3 0.523673 .... 97 20.0095 98 21.0495 99 21.5747
The last example shows how to set the beta field. This is useful because one of the coloring methods is 'Beta', which uses the beta values to color the molecule according to the current color scale. (This can also be done with the occupancy field.) Thus redefining the beta values allows you to color the molecules based on your own definition. One useful example is to color the molecule based on the distance from a specific point (for this case, coloring a poliovirus protomer based on its distance to the center of the virus (0, 0, 0) helps bring out the surface features).
proc betacolor_distance {sel point} { # get the coordinates foreach coord [$sel get {x y z}] { # get the distance and put it in the "newbeta" list set dist [veclength2 [vecsub $coord $point]] lappend newbeta $dist } # set the beta term $sel set beta $newbeta }And here's one way to use it:
# load pdb2plv.ent using anonymous ftp to the PDB vmd > mol load webpdb 2plv vmd > set sel [atomselect top all] vmd > betacolor_distance $sel {0 0 0}Then go to the graphics menu and set the 'Coloring Method' to 'Beta'.