Creating movies (.mpg) of paths using OPTIM

The aim is to take a discrete path -- just the sequence of stationary points -- and use OPTIM to fill in the gaps between the stationary points by writing extra frames of coordinates from the two approximate steepest-descent paths that set off from each transition state. Otherwise, movies made from just the stationary points would be very jerky!

N.B. what's described below works for charmm and amber9/NAB...can't say about other systems.

You need the key ingredient to be present in the same directory: a file called redopoints, which contains just the coordinates of the stationary points (minima and TSs) on the path you're interested in, in bare xyz format (like points). redopoints contains the sequence of stationary points encountered along the path, starting from the end minimum and travelling backwards (!)

Generate a redopoints file by running a DIJKSTRA <n> calculation in PATHSAMPLE (or a k-shortest-path calculation). Make sure that your two end points of interest are the minima indexed in min.A and min.B for the PATHSAMPLE calculation.

Next, set up an OPTIM NEWCONNECT job in the usual way for your system. Include a NEWCONNECT line with at least as many cycles allowed as there are transition states in the path. More cycles is a good idea, in case one of the connections gives a different minimum! Make sure files odata, input.crd / start, finish and perm.allow are present as applicable for your system.

Add the following lines to odata:

ALLPOINTS
REDOPATH
REOPTIMISEENDPOINTS

where REDOPATH tells OPTIM to read in a redopoints file (see below), and ALLPOINTS tells OPTIM to write lots of relevant output files. Also the following may be commented out or deleted, because the Hessian calculated at the end of the run is no longer necessary.

COMMENT ENDHESS
COMMENT ENDNUMHESS
COMMENT DUMPALLPATHS

Make sure there's a PATH line something like this one:

PATH 100 0.001

with a 'large' value for the first argument: the number of coordinate frames to save on either side of the TS. The second number is the threshold on the energy difference between two configurations, below which they're considered the same for the purposes of writing extra frames. This is to improve the distribution of frames along the path, so that we get more frames where the energy is changing most rapidly.

Ready to go... set off the OPTIM job.

It should be reasonably quick: no actual NEB iterations should be necessary. Having said that... check the main OPTIM output file: sometimes different connected minima are found by stepping off a particular TS than you were expecting. If this happens, try rerunning with a smaller value of PUSHOFF. If this doesn't help, then the NEWCONNECT commands in odata will tell OPTIM to try to bridge the gap by finding new TSs in the usual way to form a connected path. If OPTIM struggles to reconverge the TSs (for instance, you see 'bfgsts> WARNING *** initial eigenvalue is positive - increase NEBK?' in the output), try some or all of the following: increase the number of number of iterations allowed in the calculation of the smallest Hessian eigenvalue (1st argument on BFGSTS line, check also that the STEPS value is consistent); decrease the second argument to MAXERISE if you're also using NOHESS; decrease the convergence criterion for the calculation of the smallest non-zero Hessian eigenvalue (4th argument to BFGSTS). See also BFGSTSTOL and BFGSTSPC.

Interesting output: EofS, which you can look at with gnuplot or similar, and path.xyz.

For amber9, you can simply load path.xyz and coords.prmtop directly into VMD with

vmd -parm7 coords.prmtop -xyz path.xyz

and the atom types, connectivities etc should be picked up from the .prmtop file.

For CHARMM, path.xyz may need to be edited to give input appropriate for VMD or other protein-y viewer. The number of atoms and comment line at the top of each frame of coordinates need to be removed, as do the atom labels at the start of each line of coordinates (these can be replaced by spaces).

There are, of course, many ways of doing this, but here's my favourite :-)

sed -e '/Energy/d' -e '/ 160$/d' path.xyz > foo

Note the number of atoms in the second regular expression (after the space)... it'd probably be best to replace this with an appropriate value for your system ;-)

and, using vi in command mode, do a

:%s/^[A-Z]/ /g

on the opened file foo. You'll need to edit this if your atom labels are 2+ characters long...

Alternatively, because this can take a looong time for large files, use cut to strip the atom labels:

cut --complement -c 1-2 foo > newfilename
mv newfilename foo

Then, to convert this points format to pdb format, run

~/svn/SCRIPTS/CHARMM/points2pdb_combined.py input.crd foo

where input.crd can be any old coordinate file for this system (molecule and force field), in .crd format. It's just used as a template for the residue names and atom types.

Output from this is a file foo.pdb: coordinate frames in PDB format, separated by the END line used in PDB files.

Getting VMD to make the movie

Fire up VMD, load the coordinates file and make your protein look pretty...

In the VMD Main console, follow the menus Extensions -> Visualization -> Movie Maker

In the Movie Generator, make sure that Movie Settings -> Trajectory is selected, and that the name of the working directory and movie title are what you want. Play around with the other settings as you wish. Making a .mpg is recommended (?) - i.e. Format -> MPEG-1 (ppmtompeg) is selected, which is usually the default.

Click Make Movie, then go have a cup of tea :-) One gotcha here that caused me some consternation (!) - if you have Renderer -> Snapshot (Screen Capture) selected then you should literally go away to make tea. If you change the screen focus/change window etc then it may mess up the movie (you would see all or some parts of the protein "frozen" from the point where the focus changed.

Then, sit back and enjoy your masterpiece! I have used Gnome Mplayer to view the .mpg on my workstation (ubuntu)

--jmc49 22:52, 28 July 2008 (BST)