Analysis scripts

Get the total mass given a selection.

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

Find the min and max coordinate values of a given molecule in the x,
y, and z directions (see also the `measure` command 'minmax').
The
function takes the molecule id and returns two vectors; the first
contains the min values and the second contains the max.

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

Compute the radius of gyration for a selection (see also `measure rgyr`). The square of the radius of gyration is
defined as
.
This uses the center_of_mass function defined earlier in
this chapter; a faster version would replace that with `measure
center`.

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

Compute the rms difference of a selection between two frames of a trajectory. This takes a selection and the values of the two frames to compare.

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} { set newbeta {} # 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]] lvarpush 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'.