Calculating Rg using VMD

I am working on a system with 60 copies of a molecule which has 1287 atoms in it. The script below loads the trajectory and calculates the radius of gyration of each copy.

mol new system.gro

for {set num 0} {$num < 201} {incr num} {
mol addfile cen_md_$num.xtc waitfor all
}

for {set i 1} {$i <61} {incr i} {
set start [expr 1 + 1287 * ($i-1) ]
set end [expr 1287 * $i]
puts "$start $end"
set sel [atomselect top "serial $start to $end"]
set filename [open "rg-dextran-number-$i.txt" w]
puts "Number of frames is [molinfo top get numframes]";
for {set frame 0} {$frame < [molinfo top get numframes]} {incr frame} {
animate goto $frame
$sel update
puts $filename [format "%8.3f" [measure rgyr $sel]]
}
close $filename
}

quit

VMD 'measure hbonds' syntax

The figure below makes it easy to remember the syntax of VMD's measure hbonds command.

Contact Maps using VMD

Contact maps are a quick way to identify residues of a protein (or monomer) that interact with residue from another protein, ligand or monomer. The script below helps obtain a distance list:

set seg1 [atomselect top "segname SG03 and name CA"]

set seg2 [atomselect top "segname SG04 and name CA"]

set file [open "Contact_map.dat" w]

set list1 [$seg1 get index]

set list2 [$seg2 get index]

foreach atom1 $list1 {

        foreach atom2 $list2 {

                set index1 [atomselect top "index $atom1"]

                set index2 [atomselect top "index $atom2"]

                set resid1 [[atomselect top "index $atom1"] get resid]

                set resid2 [[atomselect top "index $atom2"] get resid]

                set resnm1 [[atomselect top "index $atom1"] get resname]

                set resnm2 [[atomselect top "index $atom2"] get resname]

     puts $file "$resnm1 $resid1 $resnm2 $resid2 [veclength [vecsub [measure center $index1] [measure center $index2]]]"

                $index1 delete

                $index2 delete

        }

puts $file " "   # I include this line to make it easier to make contour plots using gnuplot.

}

close $file

Aligning hexameric proteins along an axis

This tcl script aligns a hexamer such that its central axis is aligned along the y-axis.


set all [atomselect top protein]
$all moveby [vecscale [measure center $all] -1]

set res10 [atomselect top "resid 10"]
set res40 [atomselect top "resid 40"]

set cenres10 [measure center $res10]
set cenres40 [measure center $res40]

set vector [vecsub $cenres10 $cenres40]

$all move [transvecinv $vector]
$all move [transaxis z 90]

foreach i {A B C D E F} {
        set chain [atomselect top "chain $i and protein and resid 4 to 105"]
        $chain writepdb chain$i.pdb
}

This script is also applicable to trimers, tetramers and other aggregates that have a well defined central axis.

Manipulating dihedral angles in VMD

VMD's scripting interface is quite powerful. It has a trans bond command which gives the transformation matrix required to achieve desired dihedral rotation. If you want to edit the dihedral 1-2-3-4 (i.e. rotate around bond 2-3), then you first select the set of atoms on which the transformation is to be applied; 1,2,5,6,7,8,9 in this case. To this selection you apply the transformation;

$moveselection move [trans bond [lindex [$a1 get {x y z}] 0] [lindex [$a2 get {x y z}] 0] 10 deg]

trans bond expects only the atom indices. The value 10 indicates the magnitude and sign of change to the dihedral 1-2-3-4. deg says the change should be in degrees and not radians. To get a nice picture of rotation around a bond use this script:

set movesel [atomselect top "index 1 2 5 6 7 8 9"]
for {set i 0} {$i <360 } {incr i 10}
{
$movesel move [trans bond [lindex [$a1 get {x y z}] 0] [lindex [$a2 get {x y z}] 0] 10 deg]
render TachyonInternal [format "rotate%03d.tga" $i]
}

Making 2D heat-maps of protein secondary structure

I needed to make 2D heat-maps of protein structure with time on x-axis and protein residues on y-axis. This helps to see the time evolution of structure quite nicely. This is easily accomplished by loading the trajectory in VMD and sourcing the file ss_traj.tcl and going through your entire trajectory once. This creates a buffer with all the structure information. ss_traj.tcl creates a ss_traj.dat file which has the structure information. This ss_traj.dat file is then used by ss_plot.c, after compilation of course, to generate a .ps file with the desired output. All this is explained quite well here. I made a personal copy of the files here:
ss_traj.tcl
ss_plot.c
Sample ss_traj.dat file.

Once you have the output in .ps format you can do interesting things to change the file. The 7th line has the scaling information. The width and height of the output were not to my satisfaction; so I changed the scaling factors to adjust the display. Also the executable could crash when running due to memory requirements. Simply change the dimensions of the array in ss_plot.c according to your needs.

Using psfgen to delete unwanted atoms in simulations

psfgen's user manual can be found here.

Extracting coordinates of protein from a simulation

set dir [pwd]
mol new     ${dir}/ache.psf       type {psf}
mol addfile ${dir}/ache.pdb       type {pdb} waitfor all
mol addfile ${dir}/acheminnpt.dcd type {dcd} waitfor all

set molnum 0
set skip 1
set backframe0 [atomselect $molnum backbone frame 0]

for {set i 0} {$i < [molinfo $molnum get numframes]} {incr i $skip} {
set pro  [atomselect $molnum protein  frame $i]

# Change 'protein' to "resid N to Z" to suit your needs.
set back [atomselect $molnum backbone frame $i]
set trans_mat [measure fit $back $backframe0]
$pro move $trans_mat
set k [expr $i/$skip]
$pro writedcd ${dir}/dcd/project-pdb-$k.dcd
}

Now generate a protein-only psf file. At this stage you should have n number of dcd files. To join them all in to one dcd file,

mol new    ${dir}/ache_protein.psf
set pro_frames [expr [molinfo $molnum get numframes]/$skip ]
for {set j 0} {$j < $pro_frames} {incr j 1} {
mol addfile ${dir}/dcd/project-pdb-$j.dcd
}
set top [molinfo top]
animate write dcd {${dir}/dcd/project_protein.dcd} beg 0 end [expr $pro_frames - 1] skip 1 [molinfo top]

The last step is where the individual files are stitched together.

Making movies with VMD/PyMol

These two scripts help greatly in making and handling movies. I need scripts instead of the GUI because most of my data is on machines on the other coast and it gets rather cumbersome to pull a gigantic trajectory file just to see a quick snapshot of the system. I could make PDB (or better dcd) files of parts of trajectory I am interested in, but I will need to make movies at some point.

view_change_render.tcl
trajectory_movie.tcl
http://www.myofilament.org/documents/howto/modeling/VMD_Tcl.htm

PyMol uses get_view and "set view" commands to save and retrieve views of molecules. VMD handles things a little differently. It needs 4 parameters per view. They are [rotate/center/scale/global]_matrix. Saving and retrieving are done by processes below. They have been taken from view_change_render.tcl
I added the last part get_vp to see what the view point actually looks like (similar to get_view in PyMol).

proc save_vp {view_num} {
  global viewpoints
  if [info exists viewpoints($view_num)] {unset viewpoints($view_num)}
  # get the current matricies
  foreach mol [molinfo list] {
    set viewpoints($view_num,$mol,0) [molinfo $mol get rotate_matrix]
    set viewpoints($view_num,$mol,1) [molinfo $mol get center_matrix]
    set viewpoints($view_num,$mol,2) [molinfo $mol get scale_matrix]
    set viewpoints($view_num,$mol,3) [molinfo $mol get global_matrix]
  }
} 

 proc retrieve_vp {view_num} {
  global viewpoints
  foreach mol [molinfo list] {
    if [info exists viewpoints($view_num,$mol,0)] {
      molinfo $mol set rotate_matrix   $viewpoints($view_num,$mol,0)
      molinfo $mol set center_matrix   $viewpoints($view_num,$mol,1)
      molinfo $mol set scale_matrix   $viewpoints($view_num,$mol,2)
      molinfo $mol set global_matrix   $viewpoints($view_num,$mol,3)
    } else {
      puts "View $view_num was not saved"}
  }
}

 proc get_vp {view_num} {
  global viewpoints
  foreach mol [molinfo list] {
    if [info exists viewpoints($view_num,$mol,0)] {
puts "global viewpoints"
puts "molinfo $mol set rotate_matrix     \"$viewpoints($view_num,$mol,0)\""
puts "molinfo $mol set center_matrix     \"$viewpoints($view_num,$mol,1)\""
puts "molinfo $mol set scale_matrix      \"$viewpoints($view_num,$mol,2)\""
puts "molinfo $mol set global_matrix     \"$viewpoints($view_num,$mol,3)\""
    } else {
      puts "View $view_num was not saved"}
  }
}

Once you save this output (as a text file) it becomes very easy to retrieve the view even 

after closing down VMD. Simply copy and paste the output of get_vp.

Loops in Tcl/Tk

If loop

if { test1 } {
body1 
} elseif { test2 } {
body2 
} else {
bodyn
}

For loop

for { set i 1 } { $i <= 10 } { incr i } { # This will append all multiples of all numbers for the number n # into a variable called table set table "$table $i x $n = [expr $i \* $n]\n" }

 #Make the list
set list [list "Item 1" {Item 2} Last]
set text ""
foreach item $list {
 #Append the every item to $text
 set text "$text\n$item"
}

Foreach loop

#Make the list 
set list [list "Item 1" {Item 2} Last] 
set text "" 

foreach item $list { 
#Append the every item to $text 
 set text "$text\n$item"
 }

Making VMD detect improperly named atoms

mol reanalyze does the job.
More here; http://www.ks.uiuc.edu/Research/vmd/mailing_list/vmd-l/19657.html

Solvating molecules with VMD

To perform MD simulations with NAMD with Charmm ff, I typically use Charmm to build the initial psf/pdb of the protein and then use VMD to solvate it. Now VMD does not like any psf file generated by Charmm. It likes to its own version of psf file. For this one has to first read in a pdb file and read out a psf file. For this the following script is helpful:


package require psfgen
topology ./top_all27_prot_na.inp
segment PROT { pdb xp5helixonXaxis.out.pdb }
coordpdb xp5helixonXaxis.out.pdb PROT
guesscoord
writepsf vmd.psf

package require solvate
solvate vmd.psf xp5helixonXaxis.out.pdb  -o solvate  -s WT -t 15

Useful resource:
http://www.ks.uiuc.edu/~villa/dna/

Fitting Protein Structures in VMD

# This script can be used to first load a reference structure and
# fit the second structure/trajectory to the reference.

mol new reference.pdb

mol new enzfibsol.psf
mol addfile md2.dcd type dcd first 0 last -1 step 1 waitfor all
mol addfile md3.dcd type dcd first 0 last -1 step 1 waitfor all
mol addfile md4.dcd type dcd first 0 last -1 step 1 waitfor all
mol addfile md5.dcd type dcd first 0 last -1 step 1 waitfor all

set backbone [atomselect 0 backbone]

for {set i 0} {$i < [molinfo 1 get numframes]} {incr i 1} {

set bk [atomselect 1 backbone frame $i]
set all [atomselect 1 all frame $i]
$all move [measure fit $bk $backbone]

}

VMD script to pick a box within a box

set cellxmin [lindex [lindex [measure minmax $all] 0] 0]
set cellxmax [lindex [lindex [measure minmax $all] 1] 0]
set cellymin [lindex [lindex [measure minmax $all] 0] 1]
set cellymax [lindex [lindex [measure minmax $all] 1] 1]
set cellzmin [lindex [lindex [measure minmax $all] 0] 2]
set cellzmax [lindex [lindex [measure minmax $all] 1] 2]

for {set i 0} {$i < 100} {incr i} {
# This lets you pick a random box of size 10x10x10 in the simulation box
set xmin [expr $cellxmin + ( ($cellxmax-$cellxmin-$boxsize) * rand())]
set ymin [expr $cellymin + ( ($cellymax-$cellymin-$boxsize) * rand())]
set zmin [expr $cellzmin + ( ($cellzmax-$cellzmin-$boxsize) * rand())]
set xmax [expr $xmin + $boxsize]
set ymax [expr $ymin + $boxsize]
set zmax [expr $zmin + $boxsize]

set selection [atomselect top "(x>$xmin and x<$xmax and y>$ymin and y<$ymax and z>$zmin and z<$zmax)"]
# This gives the number of atom in that selection
set n [$selection num];
#This gives the number of non-bonded interactions in the 10x10x10 box selected.
puts $file [expr $n * ($n -1)/2 ]
}

NAMD inputs from VMD

set all [atomselect top all]
measure minmax $all
set cellxmin [lindex [lindex [measure minmax $all] 0] 0]
set cellxmax [lindex [lindex [measure minmax $all] 1] 0]
set cellymin [lindex [lindex [measure minmax $all] 0] 1]
set cellymax [lindex [lindex [measure minmax $all] 1] 1]
set cellzmin [lindex [lindex [measure minmax $all] 0] 2]
set cellzmax [lindex [lindex [measure minmax $all] 1] 2]

puts "cellOrigin              [format "%7.3f %7.3f %7.3f" [lindex [measure center $all] 0] [lindex [measure center $all] 1] [lindex [measure center $all] 2]]"
puts "cellBasisVector1        [format "%7.3f %7.3f %7.3f"      [vecsub $cellxmax $cellxmin] 0 0]"
puts "cellBasisVector2        [format "%7.3f %7.3f %7.3f"  0   [vecsub $cellymax $cellymin] 0]"
puts "cellBasisVector3        [format "%7.3f %7.3f %7.3f"  0 0 [vecsub $cellzmax $cellzmin] ]"



Pretty primitive, but it works.

Common VMD Tcl/Tk commands

mol modselect 0 0 resname BGLC and noh
mol modselect 0 1 resname BGLC and noh
mol modselect 0 2 resname BGLC and noh
mol modselect 0 3 resname BGLC and noh

mol smoothrep 0 0 5
mol smoothrep 1 0 5
mol smoothrep 2 0 5
mol smoothrep 3 0 5


mol modstyle 0 0 Licorice 0.300000 10.000000 10.000000

mol addrep 0
mol modselect 1 0 protein
mol modstyle  1 0 Licorice 0.300000 10.000000 10.000000
mol modstyle  1 0 vdw
mol smoothrep 1 1 5

mol off 0
mol off 1
mol off 2
mol off 3

mol showrep 0 0 0
mol showrep 0 0 1

color Display Background white

Calculating phi/psi values in VMD

First source the following script:
--- start ---
proc all_dihed_angle { a1 a2 a3 a4 } {
  # Delete all existing Dihdral labels so that the one we add has index 0.
  label delete Dihedrals

  # Use the top molecule
  set molid [molinfo top]

  # Add the dihedral monitor
  label add Dihedrals $molid/$a1 $molid/$a2 $molid/$a3 $molid/$a4

  # Get all values
  return [label graph Dihedrals 0]
}
--- end ---

Then run this script:

mol new enzfibsol.psf
mol addfile md7.dcd first 0 last 10 waitfor all
set phifile [open "cellulose_dihed_phi.tcl" w]
set psifile [open "cellulose_dihed_psi.tcl" w]

set oxygen [atomselect top "resname BGLC and resid 436 and name O1"]
set ox [$oxygen get index]
set h1 [expr {$ox-1}]
set c1 [expr {$ox-2}]
set c4 [expr {$ox-9}]
set h4 [expr {$ox-8}]

puts $phifile [all_dihed_angle $h1 $c1 $ox $c4]
puts $psifile [all_dihed_angle        $c1 $ox $c4 $h4]

close $phifile
close $psifile
quit

Flush VMD

While running VMD analysis scripts over large number of trajectories, it is a good idea to free up memory after every step like so:
for {set j 1} {$j < 101} {incr j}
{
mol addfile /.$j.dcd first 0 last -1 waitfor all
set ca [atomselect top "name CA"]
set cb [atomselect top "name CB"]
....
$ca delete
$cb delete
animate delete all
}
The last three lines will make life much better.

Making figures using VMD

I have lately been playing around with VMD to make pictures. After having used their "screenshot" feature to make pictures all this time I decided to try out rendering images using Tachyon and boy! what a difference. A simple example is attached. I'll let you guess which of the two was made using the Tachyon ray tracer. More info here.