Efficient dynamical correction of the transition state theory rate estimate for a flat
energy barrier
Harri Mökkönen, Tapio Ala-Nissila, and Hannes Jónsson
Citation: The Journal of Chemical Physics 145, 094901 (2016);
View online: https://doi.org/10.1063/1.4962167
View Table of Contents: http://aip.scitation.org/toc/jcp/145/9
Published by the American Institute of Physics
Articles you may be interested in
Transition state theory approach to polymer escape from a one dimensional potential well
The Journal of Chemical Physics 142, 224906 (2015); 10.1063/1.4921959
A climbing image nudged elastic band method for finding saddle points and minimum energy paths
The Journal of Chemical Physics 113, 9901 (2000); 10.1063/1.1329672
Polymer escape from a metastable Kramers potential: Path integral hyperdynamics study
The Journal of Chemical Physics 133, 184902 (2010); 10.1063/1.3493292
From transition paths to transition states and rate coefficients
The Journal of Chemical Physics 120, 516 (2003); 10.1063/1.1630572
Hybrid functionals based on a screened Coulomb potential
The Journal of Chemical Physics 118, 8207 (2003); 10.1063/1.1564060
An automated nudged elastic band method
The Journal of Chemical Physics 145, 094107 (2016); 10.1063/1.4961868
THE JOURNAL OF CHEMICAL PHYSICS 145, 094901 (2016)
Efficient dynamical correction of the transition state theory rate estimate
for a flat energy barrier
Harri Mökkönen,
1,2,3
Tapio Ala-Nissila,
1,2,4
and Hannes Jónsson
1,3,a)
1
Department of Applied Physics, Aalto University School of Science, P.O. Box 1100, FIN-00076 Aalto,
Espoo, Finland
2
COMP CoE, Aalto University School of Science, P.O. Box 1100, FIN-00076 Aalto, Espoo, Finland
3
Faculty of Physical Sciences and Science Institute, University of Iceland, VR-III, 107 Reykjavík, Iceland
4
Department of Physics, Brown University, Providence Rhode Island 02912-1843, USA
(Received 13 May 2016; accepted 22 August 2016; published online 6 September 2016)
The recrossing correction to the transition state theory estimate of a thermal rate can be difficult
to calculate when the energy barrier is flat. This problem arises, for example, in polymer escape
if the polymer is long enough to stretch between the initial and final state energy wells while the
polymer beads undergo diffusive motion back and forth over the barrier. We present an efficient
method for evaluating the correction factor by constructing a sequence of hyperplanes starting at
the transition state and calculating the probability that the system advances from one hyperplane
to another towards the product. This is analogous to what is done in forward flux sampling except
that there the hyperplane sequence starts at the initial state. The method is applied to the escape of
polymers with up to 64 beads from a potential well. For high temperature, the results are compared
with direct Langevin dynamics simulations as well as forward flux sampling and excellent agreement
between the three rate estimates is found. The use of a sequence of hyperplanes in the evaluation
of the recrossing correction speeds up the calculation by an order of magnitude as compared with
the traditional approach. As the temperature is lowered, the direct Langevin dynamics simulations
as well as the forward flux simulations become computationally too demanding, while the harmonic
transition state theory estimate corrected for recrossings can be calculated without significant increase
in the computational effort. Published by AIP Publishing. [http://dx.doi.org/10.1063/1.4962167]
I. INTRODUCTION
Transitions inducted by thermal fluctuations in atomic
systems such as chemical reactions, diffusion events, and
configurational changes are often much less frequent than
atomic vibrations. In order to estimate the rate of such
rare events, a direct dynamical simulation requires much
too long a simulation time while a statistical approach can
be applied because of the separation of time scales. The rate
theory developed by Eyring and Polanyi
1
and, in a different
form, by Pelzer and Wigner
2
for chemical reactions can be
applied in many different contexts. We will refer to this
as transition state theory (TST). It provides a method for
estimating the rate of rare events by performing a statistical
average over the fast, oscillatory motion and focuses on the
probability of significant transitions. A key concept there
is the transition state, the region of configurational space
representing a bottleneck for the transition and corresponding
to a free energy barrier. One of the key assumptions of
TST is that the transition state separating the reactant and
product regions is only crossed once. If the system makes it
to the transition state and is heading away from the initial
state, it is assumed that the trajectory ends up in a product
state for an extended period of time, before a possible back-
reaction can occur. This approximation can then be checked
a)
hj@hi.is
and corrected by calculating short time trajectories started
at the transition state to obtain the so-called recrossing (or
“dynamical”) correction which is based on the reactive flux
formulation.
3–10
It turns out that TST gives an overestimate
of the transition rate and the correction factor is κ ≤ 1. If the
transition state is variationally optimised,
3,11
the correction
factor is as close to unity as possible (for a recent review
of variational transition state theory, see Ref. 11). This two-
step procedure for obtaining the value of the rate offers a
great advantage over the direct calculation of the rate from
trajectories starting at the initial state. It can take an impossibly
long time to simulate even one such reactive trajectory, while
the trajectories needed for the correction to the TST estimate
are short since they start at the transition state. We will refer
to this approach as the two step Wigner-Keck-Eyring (WKE)
procedure.
12
The TST estimate of a transition rate requires, in general,
the evaluation of the free energy of the system using some
thermal sampling. But, a simpler approach is to apply a
harmonic approximation where the rate is estimated by
identifying the maximum energy along the minimum energy
path (MEP) of the transition. This point corresponds to a
first order saddle point on the energy surface and gives
the activation energy. The entropic prefactor is obtained by
evaluating the vibrational modes at the saddle point and at
the initial state minimum. This simplified version of TST is
referred to as harmonic transition state theory (HTST).
13,14
0021-9606/2016/145(9)/094901/7/$30.00 145, 094901-1 Published by AIP Publishing.
094901-2 Mökkönen, Ala-Nissila, and Jónsson J. Chem. Phys. 145, 094901 (2016)
A recrossing of the transition state can occur for two
different reasons. One is the fluctuating force acting on the
system due to the heat bath. If such a force is large enough
and acts in the direction opposite to the velocity of the system
soon after it has crossed the transition state, the system can go
back through the transition state to the part of configuration
space corresponding to the initial state. The other reason for a
recrossing of the transition state is related to the shape of the
potential energy surface. If the reaction path is curved near the
transition state, the system can enter a repulsive region that
creates a force on the system that sends it back to the initial
state. The more recrossings that occur, the smaller κ becomes.
The recrossings are particularly important if the energy along
the reaction path in the region near the transition state is
relatively constant, i.e., the energy barrier is flat.
A different approach to the estimation of thermal transi-
tion rates was developed by Kramers
15
and later generalised
to multidimensional systems by Langer
16
(for reviews, see
Refs. 17 and 18). There, a harmonic approximation to the
energy surface is used and a statistical estimate is made for
the recrossings due to the fluctuating force from the heat bath.
The advantage of this approach is that some of the recrossings
are taken into account in the rate estimate without requiring a
dynamical recrossing correction. A disadvantage as compared
to the two step WKE procedure is that some of the recrossings
are not included, in particular those resulting from the shape
of the potential energy surface far from the saddle point (see,
e.g., Ref. 19). Another disadvantage of the Kramers/Langer
approach is a harmonic approximation of the energy surface
in the direction of the reaction path at the transition state. As
a result, the rate is estimated to vanish if the energy barrier
is flat, i.e., if the second derivative of the energy along the
reaction path is zero. Such flat top energy barriers can occur
in various applications.
One example of a flat barrier problem is the transition
of a polymer from a potential well, the so-called polymer
escape problem.
20–27
Experimental examples of systems of
this sort include polymer translocation,
28,29
where a polymer
is crossing a membrane through a pore
30
or narrow µm-scale
channels with traps.
31
Recent experiments by Liu et al. involve
the escape of a DNA molecule from an entropic cage.
32
Similar
translocation and escape processes are common in cell biology
and have possible bioengineering applications, such as DNA
sequencing
33
and biopolymer filtration.
34
In a recent study of a model polymer escape problem,
the application of HTST followed by dynamical corrections
was shown to give accurate results as compared to direct
dynamical simulations, while the Kramers-Langer approach
gave a significant underestimate of the transition rate for the
longer polymers.
35,36
Flat energy barriers are also common
in magnetic transitions involving the temporary domain wall
mechanism.
37,38
The evaluation of the recrossing correction in the WKE
procedure can involve large computational effort for extended,
flat energy barriers. This occurs, for example, for the escape
of long polymers that are long enough to stretch between
the initial and final state energy wells while the polymer
beads undergo diffusive motion back and forth over the
barrier. Well established procedures exist for the evaluation of
the recrossing correction from trajectories that start and the
transition state and eventually make it to either the final state or
back to the initial state (see, for example, Ref. 39). However,
the diffusive motion along the flat energy barrier can make
such trajectories long and the calculation computationally
demanding.
We present here an approach for calculating the
recrossing correction in such challenging problems by using
a procedure that is similar to so-called forward flux sampling
(FFS)
40–42
where a sequence of hyperplanes is constructed
and trajectories are calculated to estimate the probability that
the system advances from one hyperplane to the next. This
method is similar to the “milestoning” approach developed
by Faradjian and Elber.
43
While the FFS method has been
proposed as a way to calculate transition rates starting from
the initial state and ending at a final state, we start the
hyperplane sequence at the transition state and evaluate the
probability that the system makes it all the way to the final
state given that it starts at the transition state. This turns out
to be a more efficient procedure than the standard method for
evaluating the recrossing correction by calculating individual
trajectories that make it all the way from the transition state
to the final state.
39
We present results on the computational
efficiency of these different methods for estimating the escape
rates of polymers with up to 64 beads, where there is a
pronounced flat barrier.
The article is organised as follows: In Sec. II, the test
problem, parameters, and numerical methods are described.
The results of the various calculations are presented and the
efficiency is compared in Sec. III. The article concludes with
a summary and discussion in Sec. IV.
II. SYSTEM AND METHODS USED
A. Description of the system
The polymer is described by a set of N identical beads
that are connected to two neighboring beads except for the
end points. The configuration of the polymer is described by
the coordinates of the beads r B {r
n
}
N
n=1
. The centre of mass
of the polymer is X
0
=
1
N
N
n=1
r
n
. The dynamics of the beads
is described by the Langevin equation where for the nth bead
at time t,
m ¨r
n
(t) + γ ˙r
n
(t) + ∇
n
[V(r
n
(t)) + U] =
2γk
B
T ξ
n
(t), (1)
where m is the mass of a bead, γ the friction coefficient,
V (r
n
) the external potential, U the interaction potential
between adjacent beads, k
B
T the thermal energy, and
ξ
n
(t) a Gaussian random force satisfying ⟨ξ
n
(t)⟩ = 0 and
⟨ξ
n
(t)ξ
m
(t
′
)⟩ = δ(t − t
′
)δ
n, m
. The interaction between adjacent
beads is given by a harmonic potential function
U =
N −1
n=1
(K/2)(r
n
− r
n+1
)
2
. (2)
The resulting contribution to the force on bead n is
− ∇
n
U = −K(r
n−1
+ r
n+1
− 2r
n
). (3)
094901-3 Mökkönen, Ala-Nissila, and Jónsson J. Chem. Phys. 145, 094901 (2016)
FIG. 1. The external potential of Eq. (4). The maximum of height
∆V = ω
2
a
2
0
/4 ≈0.56 is located at x = 0 and the symmetric minima are lo-
cated at x = ±a
0
≈±1.22. The initial state, I, is confined to the left well x < 0
and the right well x > 0 corresponds to the final state, F.
The external potential, V (x), shown in Fig. 1, is a quartic
double well
V (x) = −
ω
2
2
x
2
+
ω
2
4a
2
0
x
4
, (4)
where ±a
0
gives the locations of the minima. The energy
has a maximum at x = 0 where the curvature of the potential
energy function is ω
2
. This is the same external potential that
was used in Refs. 36 and 25. The total potential energy of the
system is
N
n=1
V (r
n
) + U.
B. Direct simulations
From direct Langevin dynamics (LD) simulations starting
at the equilibrated initial state, the escape probability can be
evaluated by observing the time it takes for the system to
reach the final state in a number of statistically independent
trajectories. A trajectory is taken to have reached the final state
when the centre of mass is half-way between the maximum
and the final state minimum, X
0
> a
0
/2. If t
i
is the time of the
escape event occurring in trajectory i, then the thermally
averaged probability that a transition has occurred after
time t is
P
esc
(t) = (1/N
traj
)
N
traj
i=1
θ(t − t
i
), (5)
where N
traj
is the number of trajectories simulated, and θ(t − t
i
)
the Heaviside step function. The escape rate is then given by
R
LD
=
dP
esc
(t)
dt
, (6)
where the derivative is determined by fitting to the linear
region in the function P
esc
(t).
C. Forward flux sampling
Forward flux sampling is a class of methods based on
a series of hyperplanes, λ
0
, λ
1
, . . . , λ
n
placed between the
initial and final states.
40–42
The rate constant is calculated
by sampling the dynamics between the hyperplanes. We will
give a brief review of the method here. For a more detailed
description of the method, the reader is referred to the review
article by Allen et al.
42
The rate constant is obtained in FFS as
40
R
FFS
=
¯
Φ
I,0
¯
h
I
P(λ
n
, λ
0
), (7)
where
¯
Φ
I,0
/
¯
h
I
is the initial flux across the first plane λ
0
towards the final state and P(λ
n
|λ
0
) is the probability that the
system reaches plane λ
n
, given it was initially at λ
0
. The initial
flux is calculated by simulating a long trajectory in the initial
state for a time t
init
and counting the number of crossings q
of the first hyperplane, λ
0
, with the normal component of the
velocity pointing towards the final state. Therefore, the initial
flux is
¯
Φ
I,0
/
¯
h
I
= q/t
init
.
The configuration at each of the crossing events of the
λ
0
hyperplane serves as an initial configuration for a new
trajectory which is run until the next interface λ
1
is reached,
or the trajectory returns to the initial state by crossing λ
0
.
The probability P(λ
1
|λ
0
) is estimated as the fraction between
the number of successful trajectories and the number of all
trajectories initiated from λ
0
. The final configurations of the
successful trajectories are used as initial points to sample the
probability P(λ
2
|λ
1
), which is then a ratio of the trajectories
that reach λ
2
to those that return to the initial state by crossing
λ
0
. The procedure is repeated until all the hyperplanes have
been sampled and the probability
P(λ
n
|λ
0
) =
n−1
i=0
P(λ
i+1
|λ
i
) (8)
can be computed.
D. HTST and recrossing corrections
The evaluation of the HTST estimate of the escape
rate requires finding the first order saddle point on the
energy surface defining the transition state and evaluating
the vibrational frequencies from eigenvalues of the Hessian at
the saddle point and the initial state minimum. In order to find
the saddle point, the nudged elastic band (NEB) method
44–46
was used to determine the minimum energy path (MEP) for
the transition. The point of maximum energy along the MEP
is the relevant saddle point.
The Hessian matrix was evaluated at the minimum and
at the saddle point using finite differences of the forces on
the beads, and the eigenvalues were calculated. The HTST
estimate of the transition rate is
13,14
R
HTST
=
1
2π
√
µ
⊥
N
i=1
λ
0
i
N
i=2
λ
‡
i
e
−∆E/k
B
T
, (9)
where µ
⊥
is the reduced mass, and λ
0
i
and λ
‡
i
are the
eigenvalues of the Hessian matrices at the minimum and
at the saddle point, respectively. The activation energy, ∆E, is
the potential energy difference between the minimum and the
094901-4 Mökkönen, Ala-Nissila, and Jónsson J. Chem. Phys. 145, 094901 (2016)
saddle point. The negative eigenmode at the first order saddle
point is labeled as i = 1 and is omitted from the product
in the denominator. In HTST the transition state is chosen
to be the hyperplane containing the first order saddle point
and having a normal pointing in the direction of the unstable
mode, i.e., the eigenvector corresponding to the negative
eigenvalue.
The recrossing correction can be estimated by starting
trajectories at the transition state and observing recrossings
of the transition state until the trajectory ends up in either
the initial or final state. Voter and Doll
39
have described
a method for computing the correction factor, κ, from an
ensemble of such trajectories. The corrected transition rate
is R
HTST+VDDC
= κR
HTST
, where VDDC stands for dynamical
correction evaluated following the procedure of Voter and
Doll.
The escape rate of polymers has previously been studied
36
using HTST followed by recrossing corrections evaluated
using the method of Voter and Doll.
39
E. Hyperplane sequence for recrossing correction
We propose here an efficient method for calculating
the recrossing correction using a sequence of hyperplanes,
analogous to the formulation of the FFS method. Here,
however, the initial hyperplane is placed at the transition
state. The sequence of parallel hyperplanes then leads to a
final hyperplane near the minimum on the energy surface
corresponding to the final state. Instead of sampling the initial
flux, we equilibrate the system within the transition state
hyperplane to generate uncorrelated samples and afterwards
assign a random velocity from the Maxwellian distribution
P(v
i
) ∝ exp(−mv
2
i
/2k
B
T) to each degree of freedom. If the
net velocity
N
i=1
v
i
of the system is negative, the velocities
are reversed v
i
→ −v
i
to describe a trajectory heading towards
the final state. The dynamical correction factor is computed
according to Eq. (8) as
κ = P(λ
n
|λ
0
) =
n−1
i=0
P(λ
i+1
|λ
i
), (10)
where the initial hyperplane λ
0
is located at the transition state
and the final hyperplane λ
n
is near the final state minimum.
We will refer to this as FFDC method for calculating the
recrossing correction.
As illustrated below, the use of a hyperplane sequence and
short time trajectories between the hyperplanes can reduce the
computational effort involved in determining the recrossing
correction factor as compared with the use of long trajectories
that make it all the way from the transition state to the final
state.
F. Parameters and numerical methods
The Langevin trajectories of Eq. (1) were calculated using
the Brünger-Brooks-Karplus integration scheme
47
with a time
step of ∆t = 0.005. The number of beads N in the polymer
chains was in the range N ∈ {1, . . . , 64}. The parameters
were chosen to be γ = 1.0 and m = 1.0. The parameters for
the external potential of Eq. (4) were chosen to be ω
2
= 1.5
and a
2
0
= 1.5. The same values of the parameters were used in
Refs. 25 and 36.
If the units of length, mass, and energy are chosen to
be l
0
= 1.02 nm, m
0
= 1870 amu, corresponding to a double
stranded DNA, and the unit of energy is k
B
T at T = 300 K,
the unit of time becomes t
0
=
m
0
l
2
0
/k
B
T = 27.9 ps. In the
direct dynamical simulations to determine the escape rate, a
total of 1000 to 240 000 trajectories were used depending on
the chain length.
In the calculations of the MEP using the NEB method,
several images of the polymer were placed between the
initial and final states and connected with harmonic springs.
The energy was then minimised using the projected velocity
Verlet integration.
44
The precise location of the saddle point
was found by minimising the force acting on the highest
energy image using the Newton-Raphson method. The spring
constant used in the NEB calculations was k
NEB
= 8.2. The
number of images, P, was typically chosen to be between 9
and 19, but for larger values of N, P = N/2 was sometimes
used.
In the FFS method, the number of hyperplanes between
the initial and final states was n = 10. The initial hyperplane
was placed at λ
0
= −a
0
and the final hyperplane λ
n
= a
0
/2.
For the longer chains, additional hyperplanes were used,
up to 20. To obtain the escape rate, 100 000 initial points
were sampled and the number of trajectories started from
each one. To obtain the statistical error in the escape rate
for efficiency analysis, the FFS calculation was repeated
with 10 000 trajectories started from each hyperplane and the
standard deviation of the result computed. The statistical error
estimate (standard error of the mean) was then evaluated as the
standard deviation divided by
√
N
s
, where N
s
is the number
of samples. This was repeated until a similar level of accuracy
was reached as for other methods for estimating the escape
rate.
In the FFDC calculations of the recrossing factor using
a hyperplane sequence, 10 000 initial configurations were
generated by equilibrating the system within the transition
state. Subsequently, 10 000 trajectories were generated from
each plane to obtain the probability in Eq. (10). The same
error estimate (standard error of the mean) was used as with
the FFS calculation. The statistical error in κ was computed
by repeated runs until the desired level of accuracy was
reached.
III. RESULTS
The escape rate was calculated at two different
temperature values. At the higher temperature, T = 1.0, the
escape rate is high enough that direct dynamical simulations
are possible. This provides a good test for the accuracy
of the various methods. The lower temperature, T = 0.5,
is more representative of a practical situation where the
escape rate is so low that a direct dynamical simulation
is not practical. The FFS method also turns out to require
excessive computational effort in that case, much more than
the HTST+FFDC approach.
