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 8 of 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 Appendix E 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 involving only those stationary points.
This use of Dijkstra's algorithm is also distinct from that used in the CONNECT algorithm (CarrTW05), with a distance-based edge-weight formulation.
For the case where the end point states contain more than one minimum, the Dijkstra analysis is performed separately for each start min in turn, to the rest of the minima in the database, and the path between the pair with the largest contribution to the SS rate constant is taken as the overall 'fastest'. For a given pair of end-point minima, the remaining ones in the A and B states are treated as intervening.
Key output:
standard output file: search for 'Largest contribution to SS rate constant' to find the contribution to the DPS steady-state (SS) rate constant from this 'fastest' path between any start and any end minimum, and the number of TS's on it. The mean first-passage times (= the reciprocal of the rate constants) in the forwards and backwards directions from a GT calculation on just the stationary points on this path are also quoted on the 'MFPT for fastest path' line. This latter time does not impose the condition of steady state on the intervening minima. It also 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). It's the rate constants corresponding to these GT times that are given on the 'rates without conditional probability' and 'rates with ' lines. The conditional probability bit refers to the factor p_x/p_X, where x=[a,b] and X=[A,B]: the probability of being in minimum x given that we're in state X.
On the 'detailed balance, ratio should' line the forwards and backwards GT rate constants are used to see how closely (or not!) the detailed balance relation (k_AB.p_B = k_BA.p_A) is obeyed, to assess the validity of the steady-state approximation used in the determination of the fastest path. The value quoted is k_AB.p_B/(k_BA.p_A) if DIRECTION BA is set (and vice-versa otherwise).
Sample chunk of output from a run with DIJKSTRA 1:
Dijkstra> 47099 minima with 1 connections or fewer will not be considered
Dijkstra> steps to A region converged in 47 cycles; maximum= 97 disconnected= 48089
Dijkstra> steps to B region converged in 38 cycles; maximum= 78 disconnected= 48089
Dijkstra> Best path for min 1 and any B minimum, k^SS B<-A 0.2648920258E-26
82 27545 78123 76027 93353 93573 103810 104873 104876 108780 113457 113455 113816 111592 111596 111595 111599 113595 113540 111670 108193 107117 114410 113282 106425 106426 109207 113864 107981 92364 92361 92362 29005 81267 49822 13885 13883 13880 82260 82261 26172 1
Dijkstra> Best path between any A minimum and any B minimum:
82 27545 78123 76027 93353 93573 103810 104873 104876 108780 113457 113455 113816 111592 111596 111595 111599 113595 113540 111670 108193 107117 114410 113282 106425 106426 109207 113864 107981 92364 92361 92362 29005 81267 49822 13885 13883 13880 82260 82261 26172 1
Dijkstra> Largest contribution to SS rate constant B<-A for any A and B is 0.2648920258E-26 for 41 transition states:
Dijkstra> Note that path is printed backwards starting with B, ending with A
E+ Ets E-
82 -555.5906153117 36918 -549.1427629896 27545 -550.7885223554
27545 -550.7885223554 105580 -543.5187142799 78123 -549.6194842919
78123 -549.6194842919 106531 -548.9845557840 76027 -551.8813095578
76027 -551.8813095578 130259 -540.4728394131 93353 -542.4178176348
93353 -542.4178176348 130266 -542.3606545656 93573 -542.3609660040
93573 -542.3609660040 145394 -534.7644579433 103810 -535.6463272297
103810 -535.6463272297 159813 -534.0226300928 104873 -534.7646896816
104873 -534.7646896816 159646 -525.0030493670 104876 -526.2912996602
104876 -526.2912996602 159359 -523.8936274377 108780 -530.2551921098
108780 -530.2551921098 159972 -524.2397843902 113457 -527.3469929975
113457 -527.3469929975 159495 -522.7281206378 113455 -522.8085028867
113455 -522.8085028867 162260 -518.9479471804 113816 -519.9294092195
113816 -519.9294092195 160045 -518.4098972167 111592 -520.7709917382
111592 -520.7709917382 156953 -520.6519499996 111596 -520.6537222572
111596 -520.6537222572 156948 -509.2607120466 111595 -510.7824074831
111595 -510.7824074831 156952 -510.7615417725 111599 -514.0064177058
111599 -514.0064177058 159697 -513.9535007538 113595 -513.9756163730
113595 -513.9756163730 159695 -513.2434291865 113540 -515.2702307728
113540 -515.2702307728 159696 -514.0369077195 111670 -516.4878259009
111670 -516.4878259009 157039 -513.3376855362 108193 -513.4635144282
108193 -513.4635144282 152012 -511.8132316150 107117 -513.5373413547
107117 -513.5373413547 161138 -510.2038266180 114410 -512.1777532197
114410 -512.1777532197 161142 -512.1643025601 113282 -512.2181171004
113282 -512.2181171004 161140 -503.8542455138 106425 -508.1494849751
106425 -508.1494849751 149226 -503.2919483110 106426 -506.8230208740
106426 -506.8230208740 159225 -506.6494193252 109207 -507.4352131723
109207 -507.4352131723 160121 -496.8941335219 113864 -497.9365618057
113864 -497.9365618057 160122 -497.2537286432 107981 -500.2176671595
107981 -500.2176671595 151632 -496.4669107556 92364 -497.0219772986
92364 -497.0219772986 128174 -496.3379764022 92361 -497.0175900124
92361 -497.0175900124 132028 -495.5050402587 92362 -495.9664182959
92362 -495.9664182959 130702 -495.2502792770 29005 -496.1568838964
29005 -496.1568838964 110518 -485.9257443857 81267 -490.5071327086
81267 -490.5071327086 110523 -490.4186741312 49822 -491.3243280411
49822 -491.3243280411 63967 -491.3031137536 13885 -491.3671659476
13885 -491.3671659476 110345 -486.3772296810 13883 -487.1372548014
13883 -487.1372548014 22162 -481.8341455402 13880 -483.2413610903
13880 -483.2413610903 112061 -476.0721186686 82260 -478.9507687081
82260 -478.9507687081 112062 -478.9115947145 82261 -479.8779631229
82261 -479.8779631229 112063 -476.4618088629 26172 -477.8753536184
26172 -477.8753536184 160106 -476.4919243777 1 -476.9018216251
Dijkstra> dumping minima and transition state coordinates to file redopoints, steps= 41
Dijkstra> detailed balance, ratio should be one if SS applies: 136.2100556
Dijkstra> MFPT for fastest path: A<-B value= 0.7030943368E+25 B<-A value= 0.2425804971E-08
Dijkstra> rates without conditional probability: A<-B value= 0.1422284248E-24 B<-A value= 412234294.2
Dijkstra> rates with conditional probability: A<-B value= 0.1422284248E-24 B<-A value= 412234294.2
main> CPU time spent in Dijkstra = 51.940 s
main> total CPU time spent in PATHSAMPLE = 101.49 s
Files min.[A,B,data].fastest, ts.data.fastest,points.[min,ts].fastest: a new database containing only the stationary points on the fastest path.
File Epath: numbered list of the energies of the stationary points on the fastest path, in reverse order from end to start. Feed into gnuplot or similar graph-plotting program! Note that the first column of data just contains the line number in this file; the second column holds the energies.
File redopoints: coordinates of the stationary points n reverse order along the path from end to start, 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, from end to start, in protein data bank format. Feed into a protein structure viewer as a trajectory file.
Another important type of calculation uses the graph-transformation (GT) algorithm to obtain phenomenological rate constants (i.e., ones that can be compared to experiment :-) in both directions (no matter which DIRECTION is specified in pathdata). The reference is TrygubenkoW06. GT rate constants are different from the ones from the above Dijkstra analysis in a number of ways: 1. the steady state approximation is not used in GT 2. all possible discrete paths through a database are accounted for in the GT rate 3. as a consequence of 2., the effect of full recrossings of TS's is included in GT. So, GT rates are for a particular whole database, whereas Dijkstra rates are for a single path in a particular database.
There are 4 variations on the GT theme coded in PATHSAMPLE, labelled GT, SGT, DGT and SDGT according to their keywords. I don't have much to add to the documentation about their relative merits; I've only used GT, which is what's described below.
Set up pathdata as for a Dijkstra calculation, but instead of DIJKSTRA i put
GT i
where i has the same meaning.
Here's a sample chunk from an output file:
... GT> Peak array size for branching ratios= 0.03897 Gb GT> For direction A<-B: k(A<-B)= 0.5156436361E-21 k(B<-A) (from detailed balance)= 742707728.3 total= 742707728.3 GT> For direction B<-A: k(A<-B) (from detailed balance)= 0.1105416573E-27 k(B<-A)= 159.2187655 total= 159.2187655 GT> detailed balance, ratio should be one if SS applies: 4664699.703 GT> non-DB rates: k(A<-B)= 0.5156436361E-21 k(B<-A) = 159.2187655 total= 159.2187655 main> CPU time spent in GT = 95.600 s main> total CPU time spent in PATHSAMPLE = 99.250 s ...
Look for the 'non-DB rates:' line for the GT-calculated rate constants in the two possible directions.
For each direction, the GT rate constant is also used to give a rate in the other direction via the detailed-balance relationship. As in the Dijkstra calculation, the value quoted on the 'detailed balance, ratio should' line is k_AB.p_B/(k_BA.p_A), but using the full GT rate constants here. This ratio is always this way up, irrespective of what's set by the DIRECTION keyword.
No other relevant files are produced by GT.
Running rate-constant calculations on a grouped database
All of the above methods can also be applied to a database that has been grouped into sets of minima and the ensembles of transition states that connect them. Just add a line for the grouping method you want (e.g. REGROUPFREE, REGROUPPE, REGROUPRATE) to a pathdata file that contains appropriate rate-constant instructions as described above. If the value of Planck's constant, in the prevailing energy units times seconds, is not 1.0 (the default) then don't forget to set it appropriately using a PLANCK value line. For CHARMM, in units of (kcal/mol) s, use PLANCK 9.546D-14
For example:
... CYCLES 0 PLANCK 9.546D-14 REGROUPFREE 5.0 DIJKSTRA 0 ...
Any fastest paths printed will now comprise sequences of groups of individual minima and the sets of TS's that connect them, and all the individual stationary points in those groups will be written to new *.fastest database files.
Care should be taken that the parameters chosen in the grouping scheme give rise to groups where, in broad terms, the inter-group rate constants are at least an order of magnitude smaller than intra-group rates. Otherwise, equilibration within a group probably can't occur before a transition out of the group happens, and the approximation doesn't hold that the contribution to the overall rate constant from the intra-group transitions is negligible compared to the inter-group ones.
--jmc49 22:39, 9 July 2008 (BST)