Calculating rate constants (GT and fastest path)
So, what do you do with the database of stationary points when you've finished the sampling part of the process? Usually, calculate rate constants (of various flavours) and examine mechanistically discrete paths with significant rate constants (i.e., make nice movies!).
Additionally, it's often helpful as you're building up the database to monitor the various rate constants that pathsample produces, by running separate one-processor pathsample jobs as described below. For safety, there should only be one running pathsample job in a directory at a time, though, so do this on a copy or between restarts of the sampling jobs.
It's also possible to ask pathsample to calculate the GT rate constant every n cycles of a sampling run (put n as the second integer on the GT keyword line in the pathdata input file), but this can be prohibitively expensive if n is small and the database is large.
pathdata should contain the usual keywords that set up information about the system: NATOMS, PERMDIST, etc., as used in the database-generating runs. However, remember to remove or comment out the sampling commands: SHORTCUT, FREEPAIRS, etc. A common gotcha: DIRECTION AB means we're interested in paths from the B state to the A state. The same is true with the subscripts on rate constants: i.e., we're using the spectroscopists' convention. The direction can also be changed, if desired, from the sampling runs.
The following keywords must be present:
CYCLES 0 TEMPERATURE x
where x is the reduced temperature, i.e. for CHARMM, it's the energy in the prevailing units of kcal/mol corresponding to the desired temperature: 0.592 kcal/mol for 298 K. The CYCLES 0 ensures that no further sampling is performed. The temperature need not be the same as that used for the sampling runs; MAXTSENERGY y can also changed (or included) in the rate constant calculations.
The different flavours of rate constant calculation...
Find the 'fastest' path: the single discrete path that makes the largest individual contribution to the phenomenological DPS steady-state rate constant (equation ?? in Wales02) in the desired direction, by including the following keyword:
DIJKSTRA i
where minima with integer i direct connections or fewer are removed from the analysis. Dijkstra's shortest-path algorithm is used with an edge-weight formulation based on log-weighted adjacency matrix elements, as described in EvansW04. It's slightly adapted from the original to ensure that the edge weights are non-negative (see an appendix of Semen Trygubenko's Ph.D. thesis). It's important to note that this analysis only compares the contribution to the rate constant from paths that don't revisit any of the minima: i.e., it doesn't include the effect of full recrossings of TS's, whereby a particular set of stationary points gives rise to an infinite number of paths.
This use of Dijkstra's algorithm is also distinct from that used in the CONNECT algorithm (CarrTW05), with a distance-based edge-weight formulation.
Key output:
standard output file: search for '??' for the contribution to the DPS steady-state (SS) rate constant from this 'fastest' path. The mean first-passage time (= the reciprocal of the rate constant) from a GT calculation on just the stationary points on this path is also quoted on the '??' line. This latter time includes the effect of recrossings (which acts to increase the resulting rate constant from the SS contribution) but does not take account of off-path transitions (which also acts to increase the resulting rate, as the presence of off-path transitions reduces a rate constant).
File EofS: energies of the stationary points on the fastest path as a function of the integrated Cartesian path length. Feed into gnuplot or similar graph-plotting program!
File redopoints: contains the coordinates of the stationary points in bare xyz format. Give as input to OPTIM to re-create the approximate steepest-descent paths, writing more frames of coordinates, to make smoother [movies] :-)
For CHARMM: file stationary.points.pdb contains the coordinates for the sequence of stationary points in protein data bank format. Feed into a protein structure viewer as a trajectory file.
More to come...
--jmc49 21:45, 2 July 2008 (BST)