κ-dynamics is an accelerated molecular dynamics method for systems with slow transitions between stable states. Short trajectories are integrated from a transition state separating a reactant state from products. The first trajectory found that leads directly to a product without recrossing the transition state and starts in the reactant state is followed. The transition time is drawn from a distribution given by the transition state theory rate and the number of attempted trajectories. Repeating this procedure from each state visited gives a statistically exact state-to-state trajectory.

A great challenge of modeling atomic motion is bridging the gap between the atomic vibration on the femtosecond timescale and the interesting rare events which can occur on the millisecond to second timescale. Several computational methods have been developed for accelerating molecular dynamics1 using approaches such as increasing the temperature,2 applying a bias potential to destabilize minima,3 and integrating parallel replicas,4 to increase the frequency of rare events. The primary assumption for these methods is that the interesting dynamics of the system can be described by rare events of short duration which take the system between stable states. Then the fast dynamics within each state can be modeled statistically in an average or approximate way while the slow dynamics between states is modeled correctly.

Specifically, the probability Pi to find the system in a stable state i satisfies the master equation,

\begin{equation}\frac{dP_i}{dt} = \sum _j - {\rm k}_{i \rightarrow j} \,P_i + {\rm k}_{j \rightarrow i} \,P_j,\end{equation}
dPidt=jkijPi+kjiPj,
(1)

where kij describes the rate (constant) to make a transition from state i to j. If kij is known for all pairs of states, Eq. (1) is a linear system of equations that can, in principle, be solved exactly. However, if the number of states is very large or infinite, one often has to resort to the kinetic Monte Carlo (KMC) algorithm,5,6 a stochastic procedure that is statistically equivalent to solving Eq. (1), but does not require calculating all possible rates in advance. KMC generates a stochastic sequence of states visited by the system and the times at which transitions between those states take place. If, at some point, the system is found in state i, the probability Pi(t) to escape this state in time t is given by an exponential distribution of the form

\begin{equation}P_{i \rightarrow }(t)={\rm k}_{i \rightarrow } \exp (-{\rm k}_{i \rightarrow \,} t),\end{equation}
Pi(t)=kiexp(kit),
(2)

where ki = ∑jkij is the total rate of escaping i. A way to generate t satisfying this distribution is to use

\begin{equation}t = \frac{\ln (1/\mu)}{{\rm k}_{i \rightarrow }},\end{equation}
t=ln(1/μ)ki,
(3)

where μ is a uniform random number on (0, 1]. Once t is determined, the product state j is chosen stochastically with a probability equal to the branching ratio

\begin{equation}P_{i \rightarrow j}=\frac{{\rm k}_{i \rightarrow j}}{{\rm k}_{i \rightarrow }}.\end{equation}
Pij=kijki.
(4)

There are two significant problems with using the above approach in molecular simulations. First, from each state i, all product states j need to be found. In a high dimensional system with hundreds of atoms, there can be a very large number of such states. Second, calculating the rates kij between these states is expensive.

One way to avoid the high computational cost of calculating kij is to use the transition state theory (TST) approximation.7,8 The TST approach requires the identification of a transition state (TS), which is a dividing hypersurface that separates reactants and products. If the TS contains the reaction bottleneck, the true rate can be approximated as the equilibrium flux through the surface in one direction. In most cases, however, it is difficult to identify a good TS a priori. Trajectories which cross the TS at short time can recross and reveal themselves as unreactive at longer times. Since TST is a variational theory that gives an upper bound on the true rate,9 the TS can be optimized to minimize the recrossings and the flux through the surface. The overestimation of the TST rate can also be calculated from the recrossings of trajectories initiated at the TS.10,11 This so-called dynamical correction factor, or transmission coefficient, κ ∈ [0, 1], is the ratio of the true rate to the TST rate,

\begin{equation}{\rm k}_{i \rightarrow } = \kappa _{i \rightarrow } {\rm k}^{\rm TST}_{i \rightarrow }.\end{equation}
ki=κikiTST.
(5)

Typically, calculation of κ requires following trajectories for short time as they escape the TS region. Voter and Doll have extended the concept of the transmission coefficient to multiple product states,12 

\begin{equation}{\rm k}_{i \rightarrow j} = \kappa _{i \rightarrow j} {\rm k}^{\rm TST}_{i \rightarrow },\end{equation}
kij=κijkiTST,
(6)

so that the true rate to each product state j can be calculated by considering the fate of trajectories that are launched from the same TS. The branching ratio can also be expressed in terms of the transmission coefficients

\begin{equation}P_{i \rightarrow j}=\frac{\kappa _{i \rightarrow j}}{\kappa _{i \rightarrow }},\end{equation}
Pij=κijκi,
(7)

where κi = ∑jκij.

The importance of κ is strongly dependent upon the dynamics of the system and the ability to choose a good TS. In solid systems, the harmonic approximation to TST, in which the TS is taken to be a plane perpendicular to the negative mode at the saddle point for a reaction, can be very accurate.13 Then, if the saddle points for all accessible reaction pathways can be found, the KMC algorithm can be used to model the long time dynamics of the system.14 However, in the general case where κ < 1, either because a good TS cannot be identified or dynamically correlated events are important, a better approach is needed to model the dynamics.

The κ-dynamics method presented here gives the exact classical rare event kinetics leaving a state i without having to identify any product state ja priori or to calculate the true rates ki →  or kij. Exact, in this context, is a trajectory that is statistically equivalent to one calculated with KMC using the true rates kij.

The κ-dynamics approach is related to other path sampling,15 interface sampling,16 and forward flux sampling methods17 but is different in that a single rare event trajectory is calculated through configuration space (wherever it goes), instead of an ensemble of trajectories to a predefined final state.

A cartoon of the κ-dynamics method is shown in Fig. 1. The solid circles represent the stable regions around the initial state i and two product states j and l. A requirement of the method is the definition of a reaction coordinate s(x), where x represents the configuration space of the system. Then a particular value s‡ defines the TS dividing surface s(x) = s‡, which separates the reactant state i from all product states.

FIG. 1.

κ-dynamics trajectories are sampled from a TS surface that divides initial state i from possible final states j and l. Trajectories are followed until one is found which starts in i and goes directly to a product state without recrossing the TS. In this cartoon, N = 6 trajectories were required.

FIG. 1.

κ-dynamics trajectories are sampled from a TS surface that divides initial state i from possible final states j and l. Trajectories are followed until one is found which starts in i and goes directly to a product state without recrossing the TS. In this cartoon, N = 6 trajectories were required.

Close modal

For the method to be exact, we need to make a transition to a product state j with the correct branching probability Pij [Eq. (4)] in a time t drawn from the correct distribution Pi(t) [Eq. (2)]. Neither of these quantities depend upon a reaction coordinate or TS, but we use a TS surface in the κ-dynamics method to find a statistically correct product state and transition time.

In the following, we first describe the method and show that it works and then give a numerical example and a discussion of the implementation details. The κ-dynamics method has the following steps:

  1. Sample a Boltzmann distribution in a TS defined by s(x) = s‡ enclosing an initial state.

  2. Select uncorrelated forward flux weighted initial points on the TS with umbrella sampling.

  3. Integrate N trajectories (forward and backward) until a “successful” one is found that goes directly to a product state without recrossing the TS and originates in the initial state.

  4. If no such trajectory is found after Nmax attempts, choose a new value of s‡ which increases the free energy of s(x) = s‡ and return to step 1;

  5. Increase the simulation time by an amount
    \begin{equation}t=\sum _{n=1}^N \frac{\ln (1/\mu _n)}{{\rm k}_{i \rightarrow }^{\rm TST}},\end{equation}
    t=n=1Nln(1/μn)kiTST,
    (8)
    where μn are independent uniform random numbers on (0, 1].
  6. Repeat from step 1 where the product state is the new initial state.

To show that this procedure gives the correct dynamics, we start with the definition of the TST rate;

\begin{equation}{\rm k}_{i \rightarrow }^{\rm TST} = \frac{1}{Q_i} \int e^{-\beta H} {\rm d}x {\rm d}p\, \dot{s}\,\delta (s-s^\ddag) \Theta (\dot{s}),\end{equation}
kiTST=1QieβHdxdpṡδ(ss)Θ(ṡ),
(9)

where Qi is the partition function of the initial state i and

$\dot{s}=\dot{x}\cdot \nabla s$
ṡ=ẋ·s is the velocity measured along the reaction coordinate. TST assumes that all trajectories with a positive velocity at the TS are reactive, which is enforced by the Heaviside step function
$\Theta ({\dot{s}})$
Θ(ṡ)
in the above definition. The TST assumption, however, is generally incorrect. To select only the reactive trajectories, we use the indicator function Γi, j(x, p) which is 1 if a trajectory from the phase space point (x, p) reaches i before j, and 0 otherwise. Each reactive trajectory is counted only once to give the true rate,11,16

\begin{eqnarray}{\rm k}_{i \rightarrow} &=& \frac{1}{Q_i}\!\int \!\! e^{-\beta H} {\rm d}x {\rm d}p\, \dot{s}\delta (s-s^\ddag) \Theta (\dot{s})\nonumber\\ [4pt]&& \times\, \Gamma _{\bar{i},{\rm TS}}(x,p) \Gamma _{i,\bar{i}}(x,-p),\end{eqnarray}
ki=1QieβHdxdpṡδ(ss)Θ(ṡ)×Γi¯,TS(x,p)Γi,i¯(x,p),
(10)

where we have chosen the crossing point that goes directly to

$\bar{i}$
i¯ (any stable state that is not i) before recrossing the TS, and when one reverses the momentum, the trajectory must start in i before any other state
$\bar{i}$
i¯
.

In the κ-dynamics method, crossing points of the TS are sampled with a forward flux weighting. An average over these crossing points can be written as

\begin{equation}\langle \cdots \rangle _{\rm FF} = \frac{\int e^{-\beta H} {\rm d}x {\rm d}p\, \dot{s}\, \delta (s-s^\ddag) \Theta (\dot{s}) (\cdots) }{\int e^{-\beta H} {\rm d}x {\rm d}p\, \dot{s}\, \delta (s-s^\ddag) \Theta (\dot{s})}.\end{equation}
FF=eβHdxdpṡδ(ss)Θ(ṡ)()eβHdxdpṡδ(ss)Θ(ṡ).
(11)

Using this notation, the dynamical correction factor is

\begin{equation}\kappa _{i \rightarrow } = \frac{{\rm k}_{i \rightarrow }}{{\rm k}^{\rm TST}_{i \rightarrow }} = \langle \Gamma _{\bar{i},{\rm TS}}(x,p) \Gamma _{i,\bar{i}}(x,-p) \rangle _{\rm FF},\end{equation}
κi=kikiTST=Γi¯,TS(x,p)Γi,i¯(x,p)FF,
(12)

which is an average over crossing points that contribute 1 if reactive and 0 if unreactive. Extending this to multiple product states is straightforward,

\begin{equation}\kappa _{i \rightarrow j} = \frac{{\rm k}_{i \rightarrow j}}{{\rm k}^{\rm TST}_{i \rightarrow }} = \langle \Gamma _{j,{\rm TS}}(x,p) \Gamma _{i,\bar{i}}(x,-p) \rangle _{\rm FF}.\end{equation}
κij=kijkiTST=Γj,TS(x,p)Γi,i¯(x,p)FF.
(13)

Then the probability of a successful trajectory reaching the final state j is

\begin{equation}\frac{\langle \Gamma _{j,{\rm TS}}(x,p) \Gamma _{i,\bar{i}}(x,-p) \rangle _{\rm FF}}{\langle \Gamma _{\bar{i},{\rm TS}}(x,p) \Gamma _{i,\bar{i}}(x,-p) \rangle _{\rm FF}} = \frac{\kappa _{i \rightarrow j}}{\kappa _{i \rightarrow }} = P_{i \rightarrow j},\end{equation}
Γj,TS(x,p)Γi,i¯(x,p)FFΓi¯,TS(x,p)Γi,i¯(x,p)FF=κijκi=Pij,
(14)

which is exactly that is required from Eq. (7).

In κ-dynamics we obtain the number of trial points, N, required to find a successful trajectory. If this procedure was repeated from the same initial state i, the average 〈N〉 could be calculated, as well as the dynamical correction factor κi = 1/〈N〉 and the escape time,

\begin{equation}t=\frac{\ln (1/\mu)}{k_{i \rightarrow }} = \frac{\ln (1/\mu)}{\kappa _{i \rightarrow } {\rm k}_{i \rightarrow }^{\rm TST}} = \langle N \rangle \frac{\ln (1/\mu)}{{\rm k}_{i \rightarrow }^{\rm TST}}.\end{equation}
t=ln(1/μ)ki=ln(1/μ)κikiTST=Nln(1/μ)kiTST.
(15)

But, as shown in the following, we do not not need to calculate κi or 〈N〉 to draw a statistically correct time of escape for the first successful trajectory leaving i.

The probability that a single trajectory from the TS is successful is κi. The probability that N attempted trajectories will be required to find a successful one is given by the geometric distribution

\begin{equation}P(N)= (\kappa _{i \rightarrow }) (1-\kappa _{i \rightarrow })^{N-1},\end{equation}
P(N)=(κi)(1κi)N1,
(16)

which specifies N − 1 unsuccessful trajectories before a successful one. Given that a successful trajectory is found in N attempts, the escape time [Eq. (8)] is the sum of N times drawn from the TST exponential distribution

\begin{equation}P^{\rm TST}_{i \rightarrow }(t)={\rm k}^{\rm TST}_{i \rightarrow } \exp \bigl(-{\rm k}^{\rm TST}_{i \rightarrow } t\big).\end{equation}
PiTST(t)=kiTSTexpkiTSTt.
(17)

The sum of N exponentially distributed random numbers gives a distribution that was derived by the Danish engineer, A. K. Erlang, when studying the network traffic in a village telephone exchange.18 For the TST times of Eq. (17), the Erlang distribution for N events is

\begin{equation}P_{i \rightarrow }(t; N)=\frac{\big({\rm k}^{\rm TST}_{i \rightarrow }\big)^N t^{N-1} \exp \bigl(-{\rm k}^{\rm TST}_{i \rightarrow } t\big)}{(N-1)!}.\end{equation}
Pi(t;N)=kiTSTNtN1expkiTSTt(N1)!.
(18)

The distribution of κ-dynamics escape times is obtained by combining the geometric probability distribution of N of Eq. (16), with the Erlang N-distribution of Eq. (18),

\begin{eqnarray}P_{i \rightarrow }(t) &=& \sum _{N=1}^\infty P(N) P_{i \rightarrow }(t; N) \nonumber \\&=& \kappa _{i \rightarrow } {\rm k}^{\rm TST}_{i \rightarrow } \exp \bigl(-\kappa _{i \rightarrow } {\rm k}^{\rm TST}_{i \rightarrow } t\big) \nonumber \\&=& {\rm k}_{i \rightarrow } \exp (-{\rm k}_{i \rightarrow } t).\end{eqnarray}
Pi(t)=N=1P(N)Pi(t;N)=κikiTSTexpκikiTSTt=kiexp(kit).
(19)

This derivation requires the expansion of the exponential in Eq. (18) and the definition of κi from Eq. (5).

We have shown that κ-dynamics trajectories have the correct product state distribution [Eq. (14)] and transition times drawn from the correct exponential distribution according to the true rate [Eq. (19)], so the method gives a correct state-to-state dynamical trajectory.

Two numerical examples are provided to demonstrate the κ-dynamics method. Both examples involve the diffusion of Al atoms on the Al(100) surface, using an embedded atom interatomic potential.19 First, the mechanisms of adatom diffusion are calculated with κ-dynamics. The TS surface was chosen to be of the bond-boost form introduced by Miron and Fichthorn,20 

\begin{equation}s(x) = \max \left(\frac{\left&#x007c;x_a-x_b\right&#x007c; - \langle \left&#x007c;x_a-x_b \right&#x007c;\rangle }{\langle \left&#x007c;x_a-x_b\right&#x007c;\rangle } \right) = s^\ddag,\end{equation}
s(x)=maxxaxbxaxbxaxb=s,
(20)

where xa are the coordinates of atom a, and the average indicates the equilibrium bond length between atoms a and b. The value s‡ is the maximum fractional stretch of the bond between any atomic pair a and b. A value of s‡ = 0.475 was used so that stable states were distinguished by a bond length changing by more than 47.5%.

The eight most probable reaction mechanisms for adatom diffusion are shown in Fig. 2(c). The exchange (2) dominates at 200 K. At 400 K, there are significant contributions from other processes, particularly (5) which is a concerted double-exchange. The branching ratio calculated from κ-dynamics is within statistical uncertainty of direct MD (see log-scale inset). The harmonic TST branching ratios are significantly different because the concerted exchange processes (3, 5, 6, and 7) have no corresponding saddle point; they are only possible because of anharmonicity in the potential energy surface.

FIG. 2.

The branching ratio for Al adatom diffusion on Al(100) using κ-dynamics is shown to be consistent with MD at (a) 200 K and (b) 400 K for the dominant mechanisms (c).

FIG. 2.

The branching ratio for Al adatom diffusion on Al(100) using κ-dynamics is shown to be consistent with MD at (a) 200 K and (b) 400 K for the dominant mechanisms (c).

Close modal

In the second example, shown in Fig. 3, a metastable square pyramidal cluster is placed upon the Al(100) surface. An initial bond-boost TS surface was chosen with a value of s‡ = 0.3. The surface was sampled with replica-exchange using a harmonic constraint.21 Up to Nmax = 100 independent trajectories were integrated to find a product state. If no successful process was found, the TS surface was pushed toward products by an amount Δs‡ = 0.05. A maximum of s‡ = 0.5 was required to find a successful trajectory in 100 attempts, corresponding to values of κi greater than 0.01. TST rates were found using umbrella sampling of the s(x) = s‡ surfaces and the weighted histogram method to find the free energy of the TS with respect to the reactants.22 Twenty-six transitions were observed before the cluster flattened onto the terrace, in a timescale of 0.9 s at 100 K. The transitions in which the atoms descend onto the surface involve concerted exchange events. In the final event three atoms simultaneously descend onto the terrace. The number of force evaluations in the κ-dynamics simulation corresponds to a MD time of 100 μs, giving an acceleration factor of 104.

FIG. 3.

κ-dynamics states (numbered) from a trajectory of a metastable Al pyramid collapse on Al(100) at 100 K.

FIG. 3.

κ-dynamics states (numbered) from a trajectory of a metastable Al pyramid collapse on Al(100) at 100 K.

Close modal

The primary limitation of the κ-dynamics method is the need to identify a suitable reaction coordinate that separates reactants from products. In the cases where the reactive events involve bond breaking, the bond stretch reaction coordinate of Eq. (20) makes sense, but, for example, in conformational rearrangements of biomolecules, the choice of a good reaction coordinate is far less intuitive. Fortunately, the assumed reaction coordinate does not need to be optimal. The κ-dynamics trajectories from the TS are typically very short for systems with low values of κ. For example, if 106 trajectories can be integrated, a TS with κ = 10−6 is acceptable. Furthermore, a flexible TS can be used to maximize κ.

Two other computational details are worth noting. First, each part of the algorithm is implemented in parallel, including sampling the TS to find equilibrium crossing points; the use of multiple replicas and replica exchange to enhance sampling; and following trajectories from the TS. Nmax is limited only by the number of processors available for the calculation. Second, the state-to-state trajectory, which is typically the most interesting aspect of the calculation, can be found at a fraction of the total computational work. Generating a trajectory requires sampling the TS and finding a single successful short-time trajectory out of each state visited. Then, as much computational effort as desired can be used to determine kTST and the time for each transition in the trajectory. An order of magnitude estimate requires less computational work than a precise time history.

This work was funded by NSF Grant Nos. CHE-0645497 and CHE-0848571 and the Robert A. Welch foundation under Grant Nos. F-1601 and F-1514. Computer resources at the Texas Advanced Computing Center were used.

1.
A. F.
Voter
and
M. R.
Sørensen
,
Mater. Res. Soc. Symp. Proc.
591
,
427
(
2000
).
2.
M. R.
Sørensen
and
A. F.
Voter
,
J. Chem. Phys.
112
,
9599
(
2000
).
3.
A. F.
Voter
,
Phys. Rev. Lett.
78
,
3908
(
1997
).
4.
A. F.
Voter
,
Phys. Rev. B
57
,
R13985
(
1998
).
5.
A. B.
Bortz
,
M. H.
Kalos
, and
J. L.
Lebowitz
,
J. Comput. Phys.
17
,
10
(
1975
).
6.
D. T.
Gillespie
,
J. Comput. Phys.
22
,
403
(
1976
).
7.
H.
Eyring
,
J. Chem. Phys.
3
,
107
(
1935
).
8.
E.
Wigner
,
Trans. Faraday Soc.
34
,
29
(
1938
).
9.
J. C.
Keck
,
J. Chem. Phys.
32
,
1035
(
1960
).
10.
J. C.
Keck
,
Dicuss. Faraday Soc.
33
,
173
(
1962
).
11.
J. B.
Anderson
,
J. Chem. Phys.
58
,
4684
(
1973
).
12.
A. F.
Voter
and
D.
Doll
,
J. Chem. Phys.
82
,
80
(
1985
).
13.
A. F.
Voter
and
D.
Doll
,
J. Chem. Phys.
80
,
5832
(
1984
).
14.
G.
Henkelman
and
H.
Jónsson
,
J. Chem. Phys.
115
,
9657
(
2001
).
15.
P. G.
Bolhuis
,
D.
Chandler
,
C.
Dellago
, and
P. L.
Geissler
,
Annu. Rev. Phys. Chem.
53
,
291
(
2002
).
16.
T. S.
van Erp
,
D.
Moroni
, and
P. G.
Bolhuis
,
J. Chem. Phys.
118
,
7762
(
2003
).
17.
R. J.
Allen
,
D.
Frenkel
, and
P. R. ten
Wolde
,
J. Chem. Phys
124
,
024102
(
2006
).
18.
A. K.
Erlang
,
Elektrotkeknikeren
13
,
5
(
1917
).
19.
A. F.
Voter
and
S. P.
Chen
,
Mat. Res. Soc. Symp. Proc.
82
,
175
(
1987
).
20.
R. A.
Miron
and
K. A.
Fichthorn
,
J. Chem. Phys.
119
,
6210
(
2003
).
21.
R. H.
Swendsen
and
J.-S.
Wang
,
Phys. Rev. Lett.
57
,
2607
(
1986
).
22.
S.
Kumar
,
D.
Bouzida
,
R. H.
Swendsen
,
P. A.
Kollman
, and
J. M.
Rosenberg
,
J. Comput. Chem.
13
,
1011
(
1992
).