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!

What's described below works for charmm, and contains charmm-specific setup information. Don't rely on this necessarily working without modification for other systems...

Set up an OPTIM NEWCONNECT job in the usual way. 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, finish and perm.allow (if applicable) are present. input.crd and finish must contain the start and end minima, in CHARMM's crd format.

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 maybe added, because the Hessian calculated at the end of the run is no longer necessary.

NOHESS
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.

You also 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 (!) redopoints* files are generated in Dijkstra and k-shortest-path rate constant calculations, for instance.

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.

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

path.xyz needs 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...

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

python ~wales/bin/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.

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