Analysis scripts

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

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
# (this assumes the url_get program was set up correctly (see ...)
vmd > mol pdbload 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'.