Analysis scripts

Following are some more examples of routines that could be used for analysing molecules. These are not the best routines to use since many of these can be implemented with the measure command, which calls a much faster built-in function.

Finding waters near a protein

This example finds the waters near the protein for each frame of a trajectory and writes out a PDB file containing those waters:

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.

Total mass of a selection

proc total_mass {selection} {
  set sum 0
  foreach mass [$selection get mass] {
    set sum [expr {$sum + $mass}]
  }
  return $sum
}

Note the curly braces after the expr command in the above example. Omitting those braces causes this script to run about three times slower! The moral of the story is: always put curly braces around the expression that you pass to expr.

Here's another (slightly 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"
}

Coordinate min and max

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 top all]
  set sx [$sel get x]
  set sy [$sel get y]
  set sz [$sel get z]
  set minx [lindex $sx 0]
  set miny [lindex $sy 0]
  set minz [lindex $sz 0]
  set maxx $minx
  set maxy $miny
  set maxz $minz
  foreach x $sx y $sy z $sz {
    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]]
}

Radius of gyration

Compute the radius of gyration for a selection (see also measure rgyr). The square of the radius of gyration is defined as $\sum_i m_i (\vec r_i - \vec r_c)^2/\sum_i m_i$. This uses the center_of_mass function defined earlier in this chapter; a faster version would replace that with measure center. Note that the measure rgyr command does the same thing as this script, only much much faster.

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

Root mean square deviation

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} {
  # 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'.