Locating Instantons in Calculations of Tunneling Splittings: The Test

Article pubs.acs.org/JCTC

Locating Instantons in Calculations of Tunneling Splittings: The Test Case of Malonaldehyde Marko T. Cvitaš*,† and Stuart C. Althorpe‡ †

Department of Physical Chemistry, Ruđer Bošković Institute, Bijenička Cesta 54, 10000 Zagreb, Croatia Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom

‡

ABSTRACT: The recently developed ring-polymer instanton (RPI) method [J. Chem. Phys. 2011, 134, 054109] is an efficient technique for calculating approximate tunneling splittings in high-dimensional molecular systems. The key step is locating the instanton tunnelingpath at zero temperature. Here, we show that techniques previously designed for locating instantons in finite-temperature rate calculations can be adapted to the RPI method, where they become extremely efficient, reducing the number of potential energy calls by 2 orders of magnitude. We investigate one technique that employs variable time steps to minimize the action integral, and two that employ equally spaced position steps to minimize the abbreviated (i.e., Jacobi) action integral, using respectively the nudged elastic band (NEB) and string methods. We recommend use of the latter because it is parameter-free, but all three methods give comparable efficiency savings. Having located the instanton pathway, we then interpolate the instanton path onto a fine grid of imaginary time points, allowing us to compute the fluctuation prefactor. The crucial modification needed to the original finite-temperature algorithms is to allow the end points of the zero-temperature instanton path to describe overall rotations, which is done using a standard quaternion algorithm. These approaches will allow the RPI method to be combined effectively with expensive potential energy surfaces or on-the-fly electronic structure methods.

1. INTRODUCTION Quantum tunneling can play an important role in molecular rate processes involving light atoms or at low temperatures.1,2 In a symmetric well system, tunneling manifests itself in the energy difference of the otherwise degenerate ground state level. Such “tunneling splittings” have been measured for a wide variety of systems, including molecular clusters, where the splittings give insight on rearrangement dynamics.3,4 Exact numerical approaches for calculating tunneling splittings include for small systems direct solution of the Schrödinger equation5−8 and for larger systems include the use of Monte Carlo sampling techniques (based on diffusion9,10 or path integral methods11,12). The Monte Carlo methods scale roughly linearly with system size and give the exact tunnelling splitting but are expensive since they require many evaluations of the molecular potential energy surface (PES). A number of approximate semiclassical methods have therefore been developed,13−18 which while not as accurate as the Monte Carlo method often give qualitative agreement with the experimental splittings. One such method is the instanton approach. It was introduced in the context of rate theory by Callan and Coleman19 and later applied to tunneling splittings by Vainshtein et al.20 In the formulation of statistical Feynman path integrals, the instanton represents the path with the dominant contribution to the quantum partition function. Thermal rate constants or tunneling splittings are obtained © 2016 American Chemical Society

from the quantum partition functions, whereby only the contribution from the dominant path is calculated and those from the harmonic fluctuations about it, which are evaluated analytically. In this way, the approach avoids the need to sample extensively the PES. Mil′nikov and Nakamura were the first to develop the instanton approach into a practical computational tool for calculating tunneling splittings,21 which they applied to a number of molecules,22 including malonaldehyde.21,23 They also demonstrated the close link between the instanton approach and WKB theory.21 Their method works in the internal molecular coordinates, which means that it is not trivial to adapt it to another system with a different geometry. The recently obtained ring-polymer instanton (RPI) method24,25 avoids this by using Cartesian coordinates. The instanton path is discretized into beads of a ring polymer (replicas or images of the system) that represent the molecular system at equally spaced imaginary time steps. The algebra to derive the instanton expression then simplifies and leads to a generally applicable algorithm, which consists of two time-consuming steps: a multidimensional minimization to find the instanton and a determinant evaluation for obtaining the contribution from harmonic fluctuations about the instanton. The RPI Received: November 10, 2015 Published: January 12, 2016 787

DOI: 10.1021/acs.jctc.5b01073 J. Chem. Theory Comput. 2016, 12, 787−803

Article

Journal of Chemical Theory and Computation method has been generalized to obtain the complete splitting patterns of systems with multiple symmetric wells.26 A number of new methods have been developed recently which make instanton calculations of reaction rates more efficient. Rommel et al.27 compared different optimization methods for finding stationary action paths and introduced unequal imaginary time steps28 into the instanton formalism to distribute beads along the path more uniformly and aid the convergence of the instanton search at low temperatures. Kryvohuz29,30 combined the instanton approach with other semiclassical rate expressions and developed an efficient numerical algorithm to calculate rates near the crossover temperature. Einarsdóttir et al.31 made the first use of chain-ofstates based methods to find instantons by applying the nudged elastic band method (NEB)32 developed initially for minimumenergy path (MEP) searches. To our knowledge, these methods have not yet been used to calculate tunnelling splittings (the only method we are aware of to increase the efficiency of the RPI method is the elegant approach of Kawatsu and Miura,33 which takes the zero-temperature limit analytically and is powerful if one wishes to converge the instanton results to better than 1%). Here, we adapt these finite-temperature instanton methods to the RPI method. We expect large increases in efficiency since these methods work by removing redundant clusters of beads near turning points, which become particularly dense in the zero-temperature limit. After summarizing the instanton theory in Section 2, we discuss in Section 3 how to apply to the RPI method the recent techniques for increasing the efficiency of instanton rate calculations; we focus on the method of Rommel et al.27 (which uses unequal time steps in the action integral) and on two methods (based on respectively the NEB and the string method) for minimizing the abbreviated (i.e., Jacobi) action integral, using equally spaced position steps.31 In Section 4, we test these methods on a two-dimensional model of malonaldehyde, which excludes overall rotations of the system. In Section 5, we describe how to include the effects of overall rotation (using quaternions to orient the end points of the instanton path) and carry out tests on a fully dimensional model PES for malonaldehyde. Section 6 concludes the article. We use atomic units throughout, with ℏ = 1.

N

S (x) =

i=1

(3)

In this paper, we characterize the system using the mass-scaled Cartesian coordinates xi,j (xi , j = mj × position of atom j), where the first subscript i denotes the ring-polymer bead number (or the imaginary time index), whereas the second subscript j denotes the degree of freedom, and f = 3 × (number of atoms) is the dimensionality of the system. In the instanton method,25 the partition function in the ringpolymer form eq 2 is evaluated by expanding the Euclidean action in the exponent to second order in x about its stationary point. For finite β, the stationary action is a first-order saddle point on the ring-polymer PES. In the tunneling splitting applications, the stationary action, with β → ∞, is a minimum. We carry out the Taylor expansion in the time- (and mass-) scaled Cartesian coordinates xi , j (di + 1 + di)/2 28 and integrate out the quadratic terms analytically by transforming to the normal modes of the ring polymer. The action minima of the full system (including tunneling) at zero temperature encompass paths that stay for the whole duration in one of the potential minima and also those that pass from one minimum to another and back any number of times. The minimum action path (MAP) that connects the two mimima is called an instanton, and it follows the classical orbit of the system with the period β in the upside-down potential. The contributions to the partition function from all the MAPs can be summed up analytically20 and expressed in terms of the single kink action Skink, which is the action of a single passage from one minimum to another, and the eigenvalues λl of its Hessian,28 Hi , j ; i ′ , j ′ =

lim

⎛ βΔ ⎞ Q (β ) ⎟ = cosh⎜ ⎝ 2 ⎠ Q 0(β)

(1)

Q (β) = Tr[e

] = lim

N →∞

∏ i

1 2πdi

∫ d xe

2 di + 1 + di di ′+ 1 + di ′

(4)

⎡ ∏ λ ⎤1/2 l Φ = ⎢ l 0⎥ ⎢⎣ ∏l λl ′ ⎥⎦ (5) ′ One zero eigenvalue in the full system always comes from the invariance under the imaginary time translation of the instanton. Other zero eigenvalues (if they exist) come from the translational and rotational invariance of the ring polymer (path). The final result of the instanton theory for the tunneling splitting20,25 is

where Q(β) is the partition function of the full system, and Q0(β) is the partition function of the system in the absence of tunneling. The expression in eq 1 is evaluated in the zerotemperature limit, with β = 1/kT, where k is the Boltzmann constant. The partition functions in eq 1 are converted into the ring-polymer form by splitting the exact partition function into N infinitesimal, and not necessarily equal,28 imaginary time steps di, where i = 1,...,N and ∑i N= 1 di = β, −βĤ

∂ 2S kink ∂xi , j∂xi ′ , j ′

in time- and mass-scaled coordinates. The matrix elements of the action Hessian on the r.h.s. of eq 4 are given in Appendix I in eq 25. If we denote the (time- and mass-scaled) Hessian eigenvalues of the action of the corresponding nontunneling system by λ0l , the final result can be expressed in terms of the ratio Φ of all the nonzero eigenvalues of the two Hessians,

2. INSTANTON THEORY In the RPI method,25 the tunneling splitting between two symmetric wells Δ is obtained from the ratio of partition functions, β→∞

2 ⎡ (x d + di + 1 ⎤ 1 i + 1, j − xi , j) ⎥ + V (xi ,1 , ..., xi , f ) i ⎢⎣ 2 ⎥⎦ di 2

∑⎢

Δ=

2 Φ

S kink −Skink e 2π

(6)

The above treatment neglects the overall rotation of the molecule and the anharmonicity perpendicular to the instanton tunneling path. The practical implementation of the instanton theory to estimate the tunneling splittings in molecules consists of two computationally time-consuming steps. In the first step, we

−S(x)

(2)

where S is recognized as the Euclidean action discretized on the imaginary time grid di as 788


Article

Journal of Chemical Theory and Computation

where the turning points of the orbit on either side of the barrier, as well as the orbit energy, are not known in advance. The instantons are then the first-order saddle points of the action in eq 7. Previous work on instantons exploited the algorithms developed for the purpose of finding transition states and applied them in the extended space of the ring polymer. In the minimum-mode following methods,37 the action is minimized along the Hessian eigenvector with the lowest eigenvalue, which lies along the instanton, using the dimer38 or Lanczos37 method and is maximized along the remaining directions. They were used in combination with the limited memory version of Broyden−Fletcher−Goldfarb− Shano (LBFGS)39 optimizer to locate instantons, which requires the evaluation of the action gradient at every iteration. All- mode following methods were also exploited in the instanton calculations previously.24,27 They include the partitioned rational function optimization40 or the trust radius-based methods,41 which use the approximate Hessians, updated according to Bofill’s scheme.41 Previous work also includes the modified Newton−Raphson optimizer,28 which does not necessarily converge to the saddle point but relies on the fact that the starting geometry is already close to it. The latter group of methods involve the calculation of the eigenvalues and eigenvectors of the (approximate) Hessian at every iteration and removal of the effect of its zero eigenvalues corresponding to the overall rotations and translations. Reference 27 compares some of the above methods in the context of the instanton search, whereas their efficiency in the transition-state search context has been tested earlier.37 In the above work on reaction rates, instanton was parametrized by molecular geometries (beads) at equally spaced time intervals di = β/N along the trajectory. It has been noted27 that such parametrization leads to the dense bead populations near the turning points, where the speed of the system is low. Clearly, it is wasteful to evaluate the potential and gradients at many similar configurations. In general, the lower the temperature is (the higher the period β), the closer the turning points of the trajectory are to the PES minimum (on either side of the barrier), where the potential is flat. The oversampling problem worsens at low temperatures. As a remedy, Rommel and Kästner28 developed the formalism of arbitrary time steps di that was introduced in the previous section. The prescription

need to locate the optimal tunneling path by minimizing the action in eq 3 of the kink. In the second step, the Hessians of the potential are evaluated along that path in order to construct the action Hessian in eq 4, which is then diagonalized to obtain the fluctuation prefactor Φ in eq 5. The evaluation of Φ in the second step can be achieved in three ways. The now obvious way25 is an explicit diagonalization of the action Hessian in eq 4 to obtain its eigenvalues. Reference 33 develops the approach further to improve its numerical accuracy and converegence properties near solution. The second approach, by Mil′nikov and Nakamura,21 develops a semiclassical formula in the β → ∞ limit based on the Jacobi fields and the Hamilton−Jacobi equation. The third approach is based on the evaluation of stability parameters of the instanton orbit or the eigenvalues of the monodromy matrix.34,35 This paper however concentrates on the methods for finding the MAPs in the first step in the following sections.

3. OPTIMIZATION METHODS An instanton is an orbit with stationary action on the upsidedown molecular potential with a set imaginary time period β. The action can be taken to be the full action S or the abbreviated (i.e., Jacobi) action SA (see below). This gives rise to two types of instanton methods, involving discretization in either imaginary time or position. We consider here both types of methods, summarizing how they are used to calculate thermal rate constants and then explaining how we must modify them to calculate tunneling splittings. The full action is given by S=

∫0

β

L(x(τ ), x(̇ τ ))dτ

(7)

where L is the imaginary time Lagrangian with the upside-down potential. If the end points of the instanton orbit are presumed fixed (as well as the end times), the stationary value of S corresponding to an instanton orbit is a minimum. Equation 7 is the time-continuous version of the discretized form in eq 3. The instanton orbit is unstable, and it traverses the potential back and forth along the one-dimensional curve between the two turning points in the full configuration space. Thus, only half of the orbit is sufficient to characterize the path and be used as active variables in the optimization. The potential at the turning point equals the total energy E of the classical orbit. The abbreviated (i.e., Jacobi) action is given by36 f

SA =

∑∫ j=1

0

β

pj(τ )xj̇ (τ )dτ =

∫x

xN

1

di = ri

pdx

2(V (x i) − E)

2(V[i − 1, i] − E)

(10)

that derives from eq 9 was used. The V[i−1,i] is the potential between the beads i−1 and i, which can be evaluated by interpolation. The Euclidean distances ri = |xi − xi−1| are taken to be equal to r for all i, and r is then used as a normalization constant, so that the di values add up to β. The energy of the trajectory is not known in advance; however, it is suggested that E = 0, the potential at the PES minimum, is used. This choice leads to somewhat denser bead population near the turning points, which aids the accuracy of the Φ calculation, if the same grid is used for representing the Hessian along the instanton. A scheme that can continuously switch between equal time steps and (approximately) equal distances along the path was also proposed,28

(8)

where the time dependence has been factored out in the last equality (and x ̇ = p in the mass-scaled coordinates). If the end points of the instanton orbit are presumed fixed and the total energy E conserved, then the stationary value of S A corresponding to an instanton orbit is a minimum. The momentum along the trajectory in the f-dimensional space can be obtained from the potential as pi =

1

(9)

Substitution of the f-dimensional momentum from eq 9 in eq 8 leads to the Jacobi’s form of least action principle. 3.1. Full Action Instanton Methods. 3.1.1. Thermal Rates. In rate calculations, we are faced with the problem of finding the classical over-the-barrier orbit of a finite period β,

di = αr 789

1 2(V[i − 1, i] − E)

+ (1 − α)

β N

(11)


Article

Journal of Chemical Theory and Computation where α is a parameter in the range [0,1]. Other representations of the instanton trajectories use expansion in Fourier series29,30 cos(2π i τ/β) of period β. The Newton−Raphson optimizer was then employed to find the optimal expansion coefficients. The method used exact Hessians calculated along the path at every iteration. The advantage of this approach was that the shape of the path was described in terms of analytic functions, and for the purpose of evaluating the action integral, any suitable quadrature can be used. In this way, the problem with accumulation of beads at turning points is avoided. References 29 and 30 use equally spaced beads in terms of Euclidean distance in the fdimensional mass-scaled configuration space and interpolate the potential at the quadrature points (at equal imaginary time steps). 3.1.2. Tunneling Splittings. In tunneling-splitting calculations, the end points of the trajectory are known and lie in the two symmetric minima set at V = 0. Thus, the energy of the orbit is E = 0. The period β is evaluated in the β →∞ limit; for practical purposes, β needs to be finite but large enough to reach convergence in order that the classical motion has enough (imaginary) time to connect the end points. All previous tunneling splitting calculations4,25,26,33 in Cartesian coordinates were based on S-minimization in eq 3 to find instantons and employed the LBFGS optimizer. The parameters of each minimization run are the duration of the kink β and the number of beads N used to discretize the orbit in equal imaginary time intervals di = β/N (along with the starting guess for the path). The minimization runs are repeated until a suitable set of parameters is found that converges the tunnelling splittings. The turning points of the orbit are kept as minimization variables. Although their geometry is known, their relative orientation with respect to the intermediate beads is not. The above-described problem of bead accumulation near the PES minima is at its most severe here. In order to improve the sampling efficiency, we are free to choose the length of imaginary time steps di in eq 11 to obtain a more uniform distribution of beads along the kink. We explore this approach numerically in Section 4. 3.2. Abbreviated-Action Instanton Methods. Although the SA-based methods are formally equivalent to the S-based methods just discussed, the methodology that must be used to optimize SA is rather different from that used to optimize S. This is because the period β and the finite imaginary time steps di are not fixed in SA-based methods but must be determined during the calculation in such a way that the beads end up distributed along the instanton trajectory (instead of collapsing to a single point describing no time evolution). To ensure that this happens, one needs to employ the sorts of methods used to calculate minimum-energy paths between known end points. We first summarize such methods (Section 3.2.1) before describing how they are applied in SA-based instanton methods (Sections 3.2.2 and 3.2.3). 3.2.1. Methods for Calculating Minimum-Energy Paths. Minimum-energy paths have widely been used to describe the transitions in physical and chemical processes from one state to another in order to understand the mechanism or as a starting point for studying kinetics. Therefore, a number of computational methods have been developed for locating MEPs. We concentrate here on the algorithms that assume the knowledge of final states of the path, the two minima of the PES, and approximate the path by a chain of states. Notable examples

that belong to this class of algorithms are the nudged elastic band (NEB),32 doubly nudged elastic band (DNEB),42 string method,43 and simplified string method.44 In the nudged elastic band method, a string of beads (images or replicas) is used to describe the path. Each bead i represents the f-dimensional molecular system xi. Images are connected by springs that ensure equal spacing of beads along the path. The calculation starts from an initial path which connects the known final states x1 and xN. Images are then evolved in an iterative process toward the MEP under the following force scheme. Spring forces are designed to act along the path, FiS = k[|x i + 1 − x i| − |x i − x i − 1|]τî

(12)

where k is a chosen force constant, and τ̂i is the unit vector tangent to the path at the image i. The force due to potential is projected to act perpendicular to the path, F⊥i = −∇Vi + (∇Vi ·τi)̂ τî

(13)

The total force on an image is Fi = F⊥i + FiS

(14)

A particular definition of tangents, the so-called upwinding tangents,45 in which the tangent on a bead is parallel to the vector that points to the neighboring bead with the higher potential, was found to improve the stability of the optimization. We use them in this work, and the definitions are given in Appendix II for completeness. Originally, the band was evolved over the PES using a damped molecular dynamics (MD) scheme, which involved zeroing the component of the velocity along the band and keeping only the component along the force if they point in the same direction at the current iteration. Other optimizers can also work efficiently with NEB, despite the fact that the projection scheme used in NEB results in nonconservative forces and a nonhermitian Hessian matrix.46,47 A comparison of different optimizers is given in ref 46. The most efficient optimizer was found to be the global LBFGS (without the line minimization)42 in which all images along the band are minimized with a single instance of the optimizer. Another choice that worked extremely well was FIRE (fast inertial relaxation engine),48 an MD algorithm with a set of rules that allows for a variable time step in the MD simulation and keeps some components of the velocity that are not parallel to the force. A modification of the original NEB algorithm, called DNEB, was later introduced42 to improve the stability for some more involved paths. In this method, a component of the spring force perpendicular to the path and the potential force is kept and multiplied by a switching factor, which zeros it near the MEP. The string method is similar to the NEB method, the only difference being the way in which images are kept equally spaced. In the string method, the images are repositioned after each step. A comparative study46 showed that a similar amount of work is needed for both the NEB and string methods to converge the MEP. More recently, a modification to the string method44 was introduced in which images are evolved in the full potential force, without the need of calculating tangent to the path, and redistributed after each step using cubic splines. Convergence is based on how much the images move in an iteration. This method was found to perform well in high accuracy 790


Article

Journal of Chemical Theory and Computation calculations,46 when the path is represented with a large number of beads, while it is unstable with small bead numbers. 3.2.2. Thermal Rates. When calculating thermal rates, the minimum-energy path methods just described cannot be applied directly to minimize SA since the end points are unknown. To get round this problem, Einarsdóttir et al. have developed the NEB method,31 such that the end points are not kept fixed but are constrained by artificial springs. Each end bead is connected by a spring to the adjacent bead, and its component parallel to the potential gradient is set to zero. In this way, the spring acts to minimize the Euclidean distance between the last two beads at either end by moving the end beads along the potential isocontour. Parallel to the potential gradient, a force is introduced that pulls the end beads toward the right isocontour at the preset energy E(∝V(x1/N) − E). The abbreviated action in eq 8 is discretized31 as N−1

SA =

⎡ pi + pi + 1 ⎤ ⎥ri + 1 ⎦ 2

∑ ⎢⎣ i=1

beads can additionally rotate, which can further deteriorate the stability of the calculation if the starting path ends are inconveniently rotated, one with respect to the other. To our knowledge, the other minimum-energy path methods discussed in Section 3.2.1 have not been used to minimize SA in thermal rate calculations. However, they would require similar modifications and would suffer similar shortcomings to the application of the NEB methods just discussed. 3.2.3. Tunneling Splittings. All the above shortcomings are avoided when minimum-energy path methods are applied to calculating tunneling splittings. The zero temperature limit in eq 1 corresponds to the trajectory of zero energy. Momenta pi are defined everywhere, and the ends are fixed at the two symmetric PES minima, apart from the rotation and translation, which will be addressed in Section 5. We also note that at zero energy, the abbreviated action in eq 8 and the action in eq 3 are equal, S = SA, at the instanton. This suggests that the minimumenergy path methods are likely to be highly effificient when used to calculate tunnelling splittings by minimizing SA. In Section 4, we investigate this numerically for the NEB and the string methods. We note that the use of such methods to calculate tunnelling splittings has not been reported previously in the literature, although Mil′nikov and Nakamura21,22 have developed an SA-based method that uses a different approach (based on expressing the path in terms of internal coordinates and expanding it in terms of a Legendre series).

(15)

where ri = |x i − x i − 1|

(16)

and pi is given by eq 9. The intermediate beads were optimized under the action of forces in eq 14, composed of the spring force parallel to the path and the gradient of the abbreviated action

4. TUNNELING SPLITTINGS IN MALONALDEHYDE: TWO-DIMENSIONAL MODEL POTENTIAL We carried out numerical tests to determine how the various methods just described improve the efficiency of the RPI method. Some of the methods could be ruled out immediately since they require evalulation of the Hessian at every imaginary time step, which is expensive. This leaves us with three types of methods: the use of variable imaginary time steps when minimizing S (Section 3.1) and two types of minimum-energy path methods for minimizing SA (Section 3.2) based on either the NEB method or the string method. We tested the methods on a two-dimensional model of malonaldehyde. This system was chosen because it excludes the complications that arise from overall molecular rotation and translation (which we deal with in Section 5 using tests on fully dimensional malonaldehyde). We measure computational efficiency by counting the total number of potential and gradient evaluations required to converge the instanton orbit. To achieve maximal efficiency, one needs to minimize both the number of “beads” (i.e., discrete steps at which the potential and gradients are evaluated) and the number of iterations needed to reach convergence. All other contributions to the computational overhead are considered to be small. The convergence criterion ϵ that we employ is the maximum absolute value of the gradient (in mass-scaled Cartesian coordinates) at a bead used in LBFGS optimization, ϵ = maxi(|∇i S |) (or S⊥ instead of S in the string method, etc.). We use the model potential by Bosch et al.49 The x coordinate in the model represents the linear reaction path coordinate that is the shortest path between the two symmetric minima of malonaldehyde. The displacement of the transferring H atom is the dominant motion in the x coordinate, and the barrier along the x coordinate (for y = 0) is 7289 cm−1. The y coordinate mainly represents the skeletal rearrangement, where a significant component comes from the stretching between hydrogen-donor and hydrogen-acceptor oxygen atoms. The y

1 1 ∇i SA = (ri + ri + 1)∇i V (x i) − (pi + pi + 1 )ri ̂+ 1 2pi 2 +

1 (p + pi − 1 )ri ̂ 2 i

(17)

where

ri ̂ =

xi − xi−1 |x i − x i − 1|

(18)

is projected perpendicular to the path, in place of the potential gradient in eq 13. A damped MD algorithm was employed to move the NEB toward MAP. Once the instanton was found, the duration of the orbit was calculated by quadrature and the beads placed at the equally spaced time intervals using interpolation in order to evaluate the rate expression. The calculated period β determined the temperature to which the rate pertains. In this way, they turned the action saddle point search for the thermal rate calculations into a minimum search. There are several minor drawbacks to the above method. First, the rate is usually desired at a certain temperature, not at a particular energy. In the above method, the temperature (period of the MAP) is not known before the MAP has been determined. Perhaps the method is best suited for the case that the rate is needed for a range of temperatures. Second, a potential problem is that there are regions in configuration space for which momentum pi is undefined. This is not expected to be a major issue since those regions are far away from the instanton orbit, apart from the end points. Finally, the NEB method works best with the fixed end points x1 and xN, whereas the above algorithm simultaneously searches for the optimal end points as well. When the end points are far from their final positions in a search, the algorithm may become unstable. The above method, MAP NEB,31 was applied to a surface scattering problem. In molecular applications, the end 791


Article

Journal of Chemical Theory and Computation value is zero at both minima, but the coordinate provides the distortions that lead to a path with a lower energy barrier of 2732 cm−1 for hydrogen crossing. The 2D potential is depicted in Figure 1 together with the linear path connecting the two

Figure 2. Convergence of tunneling splittings Δ with the number of beads N using standard RPI method (S minimization with LBFGS). Various shades of blue/red lines are results for different values of parameter α in eq 11 with the orbit duration β = 5000 and 8000 a.u., respectively. Dotted line refers to Δ calculated using an alternative formula for action (β = 8000 a.u., α = 0), which makes use of V = T at instanton, with more details given in the text. Gray line indicates the converged Δ.

Figure 1. Contour plot of the two-dimensional model potential of malonaldehyde by Bosch et al.49 Black line represents MEP; red line represents converged instanton orbit. Black circles represent the initial guess for the instanton path composed of eight beads; blue circles represent the converged eight-bead instanton path.

480 beads for β = 5000 and 8000 a.u., respectively, to converge Δ within 1% of the exact instanton result of 1.85 cm−1. Potential and gradients were evaluated at all 320/480 geometries at every iteration for the respective minimization runs. Using the convergence criterion ϵ = 10−4 a.u. with about 300 beads results in a statistical error in Δ of roughly 1%. For ϵ = 10−5 a.u. and 10−6 a.u., the error in Δ roughly drops to 0.1% and 0.01%, respectively. Empirically, it was found that the convergence of tunneling splittings Δ can be monitored by comparing Δ calculated using the action S in eq 7 to the action calculated as S = ∫ β02Tdτ, where T is the kinetic energy, shown with a dotted line in Figure 2 for β = 8000 a.u. (and α = 0). The test checks that the trajectory is Newtonian and works for a sufficiently large β to connect the two minima within the accuracy determined by the convergence criterion ϵ. 4.2. S-Minimization with Variable Time Steps. In Figure 2, we compare the above results to those obtained by using variable time steps di. At the start of the calculation, we calculate the di values at the starting guess using eq 10, with V[i−1,i] = (Vi + Vi−1)/2 and the normalization constant r adjusted to reproduce the desired β. The di values then remain constant during the LBFGS optimization. We also set α in eq 11 to 0.5 and 1 (equal spacing between beads at the starting geometry). It is shown in Figure 2 that the tunneling splittings with α = 0.5 converge faster than with α = 0 (equal time steps) initially. The 1% convergence is reached using N = 300 and 320 beads for β = 5000 and 8000 a.u., respectively. With α = 1, the convergence of Δ with N is worse than with the equal time steps di (α = 0). The reason for that is that the calculation of the prefactor Φ in eq 5 requires many discretization points (beads) near the instanton turning points. Using α = 0.5 is a compromise between the equal spacings that represent the path more faithfully and the requirements on the quadrature accuracy for calculating Φ near the turning points. Nevertheless, if a high accuracy result for the splitting is required, there are no significant gains in using α ≠ 0.

potential minima, used as a starting guess in most calculations below, the instanton path (MAP) and the MEP. The cornercutting effect and the skeletal rearrangement are both shown to have a significant effect on the tunneling dynamics. Reference 49 provides the tunneling splitting value of 1.54 cm−1 for the model potential, obtained by means of an exact quantum approach. 4.1. RPI Calculations. We first calculate the tunnelling splitting for the 2D malonaldehyde model using the standard RPI method25 (in which one minimizes S in eq 3 with equal imaginary time steps di = β/N). The starting guess for the instanton (kink) is taken to be the linear interpolation between the two potential minima, shown in Figure 1, and the path is represented using equally spaced beads along the starting guess. All N beads are kept as minimization variables and LBFGS is employed to carry out the minimization. The kink duration β and the number of beads N are converegence parameters of each minimization run. In order to find the optimal β and N, a table is constructed in ref 25, where the results are shown for various β (columns) and N (rows), with the imaginary time step β/N kept constant along the diagonals. In this way, β is found, along the diagonal, that gives sufficient time for the trajectory to connect the two minima, and subsequently, N is increased to reach the converegence with respect to the time step. The procedure obviously involves many minimization runs. We chose two representative orbit durations β = 5000 and 8000 a.u., both sufficient to reach the minima within our convergence criterion ϵ and varied the number of beads N starting from N = 4 in steps of 4. The converged path with N beads was taken as the starting path in the next minimization with N + 4 beads, along which the beads were again distributed equally spaced at the start. The convergence of tunneling splitting Δ with the number of beads N is shown in Figure 2 (with α = 0). It is shown that Δ converges slower for β = 8000 a.u., as more beads are needed to reduce the time step β/N to converge the action S and prefactor Φ. We needed N = 320 and 792


Article

Journal of Chemical Theory and Computation It can be noticed in Figure 2 that the dependence of Δ on N for α = 1 is less sensitive to the value of β. This does not come as a surprise because time steps di adjust to represent the orbit irrespective of its duration. In Figure 3, the number of iterations to converge Δ in above calculations is shown as a function of the number of beads N.

in the middle of the orbit. The beads then redistribute so that those near the turning points absorb the additional (imaginary) time. All the beads then initially move out of the position, although they lie on top of the converged instanton. Similarly, an initially uniform distribution of beads for α = 0 does not correspond to equal time steps di, and the beads move out of the position at the start of the optimization. We tried to correct this behavior by using di in eq 10 with the calculated ri = r at the starting guess, rather than determining r by scaling to reproduce β. We also placed the first and the last bead in the minima and did not include them in optimization. In this way, we stopped the growth of the number of iterations with N at large N, as is shown by black lines in Figure 3. To verify that the slow convergence of Δ with N in Figure 2 is indeed due to slow convergence of the prefactor Φ, we plot the action convergence in Figure 4. The action S converges

Figure 3. Number of iterations to converge tunneling splittings Δ using tolerance ϵ in the standard RPI calculations for different values of α in eq 11 and orbit durations β of 5000 a.u. (blue) and 8000 a.u. (red), as a function of number of beads N. Parameters are defined in the text. The initial guess in the optimizations is the converged path from the previous calculation with N − 4 beads, apart from the results shown with red circles, which use linear path connecting PES minima (empty circles) and MEP (full circles) as the starting guesses. Black lines use estimates of di based on the starting path without scaling to β, as described in the text. Figure 4. Minimum action S calculated using the converged N-bead path and eq 3. Solid red/blue lines are based on S minimization with orbit duration β of 5000 and 8000 a.u., respectively. Results for α = 0 and 1 in eq 11 are shown. Dashed lines show fully converged action calculated on the interpolated N-bead path, using cubic splines, and a dense grid on the interpolated potential, as described in the text. Inset shows the convergence of the path measured by the r.m.s. error in eq 19, with the number of beads N.

Solid lines of different shades of blue show how the number of iterations vary for different convergence criteria ϵ. It is shown that by increasing α to 0.5 and 1, the number of iterations increases for the same N and ϵ. We attribute this to the imbalance in the magnitude of “spring” constants connecting the beads, ∝ 1/di, in eq 3, which produces the disparity of different components i in the action gradient ∇iS. The increase in the number of iterations to reach convergence for α = 0.5 neutralizes, in terms of the number of potential evaluations, the faster convergence of Δ with the number of beads N (depending on the exact value of β). If highly accurate results are required, using the prescription of ref 28 to calculate di gives almost no advantage in tunneling splitting calculations. Again for α = 1, the dependence of the number of iterations to converge Δ on the orbit duration β becomes less pronounced. Perhaps surprisingly, the number of iterations grows with N for all the dependencies shown in Figure 3 (in color), although the starting guess for a calculation with a large N is already very close to the exact instanton. In fact, the dark red circles in Figure 3 show that roughly the same number of iterations are required if one uses alternative starting guesses, either a line connecting the PES minima (linear path) or the MEP, i.e., using a better starting guess does not provide an advantage in the α = 0 (or α = 1 (not shown)) calculations. We determined the values of di based on the starting path (for α = 1) and then scaled it in order to obtain the desired β. If the period β at the starting guess, calculated using ∑i r / 2V[i − 1, i] , is smaller than the desired β in our calculation, the di values are scaled to larger (than appropriate) values, which results in spreading of beads

faster for α = 1 than for α = 0 (unlike Δ in Figure 2). The inset in Figure 4 shows the convergence of the path with N, where we defined the root-mean-square deviation as (

∫0

1

|x(t ) − x inst(t )|2 dt )1/2

(19)

where t ∈ [0,1] parametrizes the trajectory between x1 and xN uniformly in terms of Euclidean distance. The integral in eq 19 thus measures the distance between the final path x with N beads and the converged instanton xinst (with N large). We evaluated the intermediate positions along instanton using the natural cubic spline interpolation in parameter t for each Cartesian coordinate separately. It is shown that the convergence of the r.m.s. deviation is even faster than that of the action S, which suggests that we should converge the action integral using the interpolated coordinates and the interpolated potential on a dense grid along the instanton. The resulting interpolated action integrals are depicted in Figure 4 using dashed lines. We needed 12 and 16 beads in order to fully converge S within 0.1% for β = 5000 and 8000 a.u., respectively. Indeed, there is no reason to use the same discretization points for the path search and the calculation of the prefactor Φ 793


Article


the starting guess is used to determine the length of time steps di for the minimization. In this way, the final path is obtained with the time steps di ideally suited to the starting guess. Moreover, using the above strategy, we need to repeat the minimization process in order to improve on the values of di, which are based on the starting guess, to be able to represent the path with smaller number of beads. It would be preferable to adapt the times steps in course of the minimization process. For this purpose, we devised an algorithm in which we recalculate the di values at every iteration in the following way. The turning points x1 and xN are set at the PES minima and are not used as optimization variables. The time steps d1 and dN+1 are set to zero, whereas d 2 = 2(V1 + V2)/2 = V2 , with dN defined likewise. The remaining di values are chosen in such a way as to zero the component of action gradient ∇iS along the path, τ̂i∇iS = 0. Using eq 3, we obtain

either. Once we have converged the path with a given number of beads N, we can use the one-dimensional cut through the PES along the final path to calculate the action S and the prefactor Φ with the larger number of beads. In order to redistribute new beads along the one-dimensional (1D) cut, it is not sufficient to use the interpolation in imaginary time as it was done in ref 31 for rate calculations. The clustering of beads near turning points requires very high accuracy for bead positions. Instead, we perform another instanton search with equal imaginary time steps (α = 0) along the 1D cut in the interpolated potential. This requires no additional potential evaluations. In order to represent the potential along the 1D cut, we construct a parabola at each turning point that has a minimum at the PES minimum and reproduces the potential at its nearest bead. The potential between the two beads, nearest to the PES minima, are then joined using cubic splines, where end point derivatives are set to match those of the two parabolas. In this way, the potential is extrapolated by the parabolas to bring back any bead that escapes the region between the minima during the minimization in the 1D potential. The prefactor Φ is then calculated by diagonalization using interpolated Hessian matrix elements from the Hessians evaluated at N original bead positions. Convergences of tunneling splittings evaluated in this way with the number of beads N are shown in Figure 5. When

⎡ xi , j − xi − 1, j

∑ τi ,j⎢ ⎣

j

di

+

xi , j − xi + 1, j di + 1

+ ∇j Vi

di + di + 1 ⎤ ⎥=0 2 ⎦ (20)

where τi,j is the j component of the unit vector tangent to the path at bead i. For each i > 2, the condition in eq 20 forms a quadratic equation used to calculate di+1 in terms of di. If eq 20 has no solution for a di, we either use the di values from the previous step or calculate the di values by minimizing the overall action gradient along the path, th

∑ |τi∇̂ i S|2 i

= min (21)

The time steps di calculated in one of the above ways are then used to calculate action gradients using eq 3, and the step is taken using the LBFGS optimizer. After we have taken the step, the beads are redistributed along the path to restore the desired spacing of beads along the path. We then evaluate the potential and gradient at beads and restart the process of determining the new di values. The reason why we have not used the obvious choice for time steps di = Vi − 1 + Vi above is that it is not entirely consistent with equal beads spacing, which is forced through interpolation, so it fails to converge. The above algorithm introduces nonconservative forces in the LBFGS minimization, similar to the chain-of-states based methods46 that we discuss in the next subsection. The tests have shown that time steps need indeed be changed at every iteration, so that the object function changes as smoothly as possible. In fact, since we use ordered set of beads that define the path, as well as path tangents to define the time step sizes di, the above algorithm can be classed as a chain-of-states based method. We applied the above algorithm to the 2D model potential of malonaldehyde to calculate tunneling splittings in Figure 6. The starting path for each S minimization with N beads was a straight line connecting two PES minima. The action S and the prefactor Φ were calculated by performing a full calculation on a 1D cut through the PES along the instanton using the interpolated potential and Hessian matrix elements in the manner described above for α = 0. This requires no additional potential evaluations. We compared the efficiency of distributing beads (i) equally spaced at every iteration (black line), (ii) equally spaced only when a bead spacing differs by more than 10% from average (red), and (iii) equally spaced apart from the first and last

Figure 5. Convergence of tunneling splittings Δ with the number of beads N. Converged N-bead instanton path is obtained by S minimization based on the starting guess from the previous optimization with N − 4 beads (linear path for N = 4). Results using various values of the orbit duration β and parameter α in eq 11 are shown. Tunneling splittings are obtained by a subsequent standard RPI calculation on the interpolated path (and potential), based on the converged N-bead path, using 2000 beads and the orbit duration of β = 8000 a.u. Circles show the tunneling splittings calculated based on the N-bead path obtained by S minimization with time steps di calculated at the starting guess without scaling to β, as described in the text. Starting path for optimization with N = 4 beads is MEP (empty circles) or linear path (full circles).

α = 1 is used, the tunneling splitting Δ converges within 1% for N = 16 for β = 5000 and 8000 a.u. If we use the calculated time steps di without the scaling to produce the desired β, the convergence is fast albeit a little slower (shown with circles in Figure 5). The convergence criterion ϵ = 10−4 a.u. with 20−30 beads results in a statistical error in Δ of roughly 0.3%. The above examples use the converged path from the previous calculation (with N − 4 beads) as a starting guess, and 794


Article


Figure 6. Convergence of tunneling splittings Δ with the number of beads N. Tunneling splittings are calculated by S minimization using variable time steps di. Beads are redistributed: (i) equally spaced at every iteration (black line), (ii) equally spaced only when a bead spacing differs by more than 10% from average, or (iii) equally spaced, apart from the first and last spacing that is set to half the value of the remaining spacings. If eq 20 for calculating di fails, either old di values are kept (solid lines) or new di values are calculated by minimizing eq 20 (circles). Inset shows the number of iterations to reach convergence in the S minimization using LBFGS. Solid, dashed, and dotted black lines show the exact instanton result and 1% and 3% deviation limits, respectively.

Figure 7. Number of iterations to convergence Δ with tolerance ϵ (defined in the text). Results obtained by S minimization with time steps di determined at the starting guess, taken from the previous calculation with N − 4 beads, are shown with lines in shades of blue and red for β = 5000 and 8000 a.u., respectively (α = 1 in eq 11). Results obtained by S minimization with variable time step sizes di are shown with lines in shades of gray (case (i) in the text).

In summary, the proposed algorithm makes use of the knowledge of the molecular geometry at PES minima. It does not require the knowledge of imaginary time duration of the kink β, although it relies on S minimization, and the final path, for a given accuracy ϵ and the number of beads N, does not depend on the starting path (through the choice of di values). This algorithm perhaps comes at a cost of degraded stability of LBFGS optimization, introduced through nonconservative forces, but it does not show in simple applications such as the calculation of tunneling path of malonaldehyde. Tunneling paths are usually simple orbits, and this is not expected to present a problem in future applications. 4.3. SA Minimization. We next determine the instanton orbit based on SA minimization in eq 8. The energy is set to zero in eq 9, and the end points of the orbits x1 and xN remain fixed in the PES minima. The remaining N − 2 beads are placed equally spaced in terms of the mass-scaled Euclidean distance along the line connecting the minima. The time is factored out from eq 8 and β is not a parameter in the optimization process. In order to evaluate the integral in eq 8, we use the trapezium rule in eq 15, with the potential evaluated at bead positions. We aim to keep beads approximately equally spaced throughout the minimization to maintain a good description of the path with small N. Perhaps the simplest way to establish this is to add an artificial harmonic potential, which possesses a minimum at equal bead spacings, to SA and minimize

spacing that is set to half the distance of the remaining spacings (blue) in order to improve accuracy near turning points. Solid lines in Figure 6 are calculated by keeping old di values if the calculation of new di values, using eq 20, fails. Circles are calculated by determining di values using the condition in eq 21, if using eq 20 fails. It is shown that convergence of Δ within 3% is achieved using (i) 13, (ii) 12, and (iii) 11 beads, irrespective of the way we calculate di values. The 1% convergence is reached using (i) 16, (ii) 16, and (iii) 15 beads. The inset shows the number of iterations to reach convergence for a given number of beads N. The different variants of the method all work with similar efficiency. The calculation with N = 16 beads converges in 22−27 iterations for ϵ = 10−5 a.u. The number of potential evaluations for locating instanton in that case is between 22 × (16 − 2) = 308 and 378. We assumed that the geometry of the PES minimum is known. The orbit duration β is not a convergence parameter, and the value obtained from optimization varies between 3530 a.u. for 13 beads and 4780 a.u. for 40 beads. Hessians were evaluated at N bead positions, and interpolated matrix element values were used in the calculation of prefactor Φ. In terms of the number of potential evaluations, the above algorithm improves by 2 orders of magnitude on the S minimization with equal di values and Φ calculated on the same grid, as shown in Figures 2 and 3. It is difficult to make a fair comparison with the S minimization using α = 1 and Φ calculated separately on a dense grid, shown in Figure 5, as those calculations use a better starting path (from a previous optimization). Nevertheless, in Figure 7, it is shown that even though the above algorithm (case (i) is shown) uses an inferior starting guess, for the same accuracy ϵ and number of beads N, the number of iterations for the above algorithm is consistently lower, irrespective of β, whereas the number of beads N needed for convergence of Δ is approximately equal.

N−1

S NEB =

⎡ pi + pi + 1 ⎤ ⎥ri + 1 + ⎦ 2

∑ ⎢⎣ i=1

N−1

∑ i=2

1 k(ri + 1 − ri)2 2

(22)

using LBFGS. The spring constant should be big enough to ensure equal bead spacings and small enough not to interfere with the minimization of the first term in eq 22. The alternative is to use the chain-of-states based methods, described in the previous section, to minimize the abbreviated action SA in eq 15 directly, while redistributing beads using interpolation, in the string method, or introducing the artificial spring force, eq 12, in the NEB method. In the string method, we use linear interpolation, rather than cubic splines, to redistribute beads after each iteration in order to add to the 795


Article

Journal of Chemical Theory and Computation stability of the method. In both the string and NEB methods, we use the upwinding tangents45 in the definition of forces, eqs 12 and 13, for improved stability. We determined the instanton orbit in our 2D model potential by minimizing SNEB in eq 22 and by minimizing SA in eq 15 using the string and NEB methods with the LBFGS optimizer. Subsequently, we calculated tunneling splittings using the standard ring-polymer instanton method on the interpolated 1D potential and Hessian matrix elements, as described in the previous section. The convergence of Δ with the number of beads N is shown in Figure 8 and compared to

Figure 9. Number of iterations to converge tunneling splittings Δ using tolerance ϵ with (defined in the text) with LBFGS vs number of beads N. Results obtained using SNEB-minimization in eq 22 are shown for spring constant k = 0.003 a.u. (green line) and the “optimal” spring constant (dark green), as defined in the text. Results obtained by SA minimization using NEB method (k = 0.03 a.u.) are shown using solid red lines. Circles (lines) show results obtained by applying the convergence criterion to the perpendicular (full) NEB gradients. NEB results using the “optimal” spring constants are shown using solid dark red lines. The “optimal” spring constant in the NEB method is shown above the frame for each N. Results obtained by SA minimization using string method are shown using solid lines in shades of blue. Blue crosses show results obtained by S minimization with variable time steps, as in Figure 6 (case (i)).

Figure 8. Convergence of tunneling splittings Δ with the number of beads N. Results using the minimization of SNEB in eq 22 (green line), SA minimization using NEB (red line) and string (blue circles) methods, and S minimization with variable di values (purple) are shown. Solid, dashed, and dotted black lines show the exact instanton result and 1% and 3% deviation limits, respectively.

10−4−10−3 a.u. in steps of 10−4 a.u., k = 10−3−10−2 a.u. in steps of 10−3 a.u., and so on. The results which converge to unequal bead spacings (deviation of more than 10% from average) were disregarded as the spring constant was deemed to be too low to serve its purpose. Results shown with the red line are obtained using k = 0.03 a.u., which is near optimal over the whole range of N studied. The convergence criterion of ϵ = 10−5 a.u. was used as a limit for the full NEB gradient at beads. Red circles show the number of iterations with the criterion ϵ applied to only the perpendicular component of the SA gradients for a direct comparison with the string method. The number of iterations in the SNEB minimization is considerably higher than in the NEB method. Results using the optimal spring constant k (dark green) and k = 0.003 a.u. (green) for all N are both shown. We also implemented the improved string method,44 which avoids the use of tangents by propagating the string using the full gradients with LBFGS. The beads are kept at equal spacings by interpolation. We used cubic splines with arc length as a parameter. As noted in ref 44, the method is unstable for low numbers of beads or performs worse than the standard string method. At high bead densities, in our case at ≈40 beads and above, it achieves similar efficiency to the standard string method. Since we often obtain satisfactory results for N < 40, we prefer to use the standard string method. The string method is parameter free and came out as a preferable choice. Other optimizers, apart from LBFGS, have been tested in the past in connection with the NEB and string methods in the MEP finding applications. Initial applications of NEB employed a damped molecular dynamics scheme.32 More recently, ref 46 found that FIRE48 is another good choice for

the S minimization results with variable di values, obtained in Figure 6 (case (i) is shown). We needed 15 and 17 beads to converge Δ within 3% and 1%, respectively, for all three minimization schemes proposed above. Figure 9 compares the number of iterations needed to converge Δ using the above methods based on S A minimization. The S minimization (with variable di values) results (crosses) are also shown for comparison. The string method and S minimization method converge comparably fast. For 1% accuracy, the string method converges with (17 − 2) × 25 = 375 potential evaluations. The NEB results perform significantly worse, with 1% convergence limit reached with 17 beads in 38 iterations (570 potential evaluations). We believe that the reason for the S minimization and string method to work better than NEB is that the redistribution of beads at every iteration makes the minimum search for each bead onedimensional in the direction perpendicular to the path. In the NEB, there is an additional computational effort to search for minimum along the path in order to place beads at equal distances. In the multidimensional case, the redistribution of beads does no longer specify uniquely the direction in which to search for the minimum, and the advantage is lost, as is shown below. The number of iterations obtained using NEB with the optimal spring constant k are shown in dark red in Figure 9. The optimal spring constant k is taken to be the one that results in the fewest number of iterations to converge Δ. It is plotted in the bar graph in Figure 9 as a function of N. We tried using k = 796


Article

Journal of Chemical Theory and Computation the optimizer in the MEP search, albeit somewhat less efficient than LBFGS. We tested the above molecular dynamics-based optimizers and used them to minimize the action SA in the string method. In QuickMin (QM),46 a damped MD scheme, the system is evolved in time, in steps of Δt, using the action gradients as forces. After each step, the component of the velocity that is parallel to the string is zeroed at every bead. Whenever the beads become unevenly distributed, for example, when the distances deviate by more than 10% from the average spacing, the beads (and velocities) are repositioned to equal spacings by linear interpolation. Additionally, if the velocity is pointing in the direction opposite to the force, it is zeroed. Here, we used the collective Nf-dimensional velocity v and force to check for negative sign in the dot product ∇S⊥v. Applying the velocity criterion to every bead separately proved to be extremely inefficient. Similarly, propagating the string with FIRE algorithm that uses variable time steps proved to be inefficient or unstable. If we apply the FIRE criteria for varying the time step to each bead separately, the scheme becomes extremely inefficient, whereas if we implement it globally for Nfdimensional vectors, the time step increases to values that lead to failure of the method. The Jacobi action is a product of pi that favors regions of lower potential and pushes the string toward MEP and ri that favors shorter paths and pushes the string toward the linear path between the two minima. The total force on a particular bead depends on two adjacent beads, and the beads need to move in accord with each other to keep the two opposing components of the force balanced throughout the propagation. If the time step becomes too big, in our case at about 5 a.u., the string buckles and leads to failure of the method. If the velocity criterion is applied to every bead, instead of collectively, the velocities are zeroed and time step reduced too frequently. In the MEP searches, forces at all beads point toward the lower potential and are independent of the adjacent beads, apart from defining the tangents. Both QM and FIRE work well in MEP searches. In Figure 10, the efficiency of MD-based optimizers is compared to LBFGS in both MEP and instanton searches. The starting guess is a linear path (of length 70 m1/2 e a0) connecting two minima. The r.m.s. distances, defined in eq 19, of the starting paths to the converged instanton and MEP are 26.0 and 46.1 m1/2 e a0, respectively. The lengths of the instanton and MEP paths are 120 and 174 me1/2a0, respectively. The convergence criterion based on the norm of the gradient at the beads is inappropriate for the comparison of MEP and instanton searches, as gradients have different units for the two cases. Instead, we define ϵ as the average Euclidean distance to the fully converged path for the particular number of beads N. We experimented with the QuickMin optimizer using different sizes of time steps. The results in Figure 10 are shown for Δt = 3 a.u. in the instanton search. Using a greater time step leads to instabilities (for some N), whereas smaller time steps lead to inefficient searches (in terms of the number of iterations). The MEP results are shown for the FIRE method (with FIRE criteria applied to each bead separately). The FIRE method is not very sensitive to the initial time step, in comparison with the QM optimizer. Results are shown for the initial time step Δt = 90 a.u., and a similar time step leads to the best results with the QM optimizer, which is then of similar efficiency to FIRE for this system. The main advantage of using FIRE over QM for this system is that the choice of the initial

Figure 10. Number of iterations needed to converge tunneling splittings and MEPs using string method and tolerance ϵ, based on the Euclidean mass-scaled distance to the exact N-bead path, as defined in the text, as a function of number of beads N. Solid lines refer to instantons; dashed lines refer to MEP. Results shown in shades of red are obtained using LBFGS optimizer; results shown in shades of blue are obtained using QM (for instanton path; using Δt = 3 a.u.) or FIRE (for MEPs).

time step is not crucial then, as it is difficult to predict the optimal size in advance. Figure 10 shows that it is more efficient to use LBFGS than an MD-based optimizer for this simple system for both instanton and MEP searches. For an MD-based optimizer, finding the instanton is costlier than finding the MEP, although the distance of the starting path to the instanton is smaller. Therefore, the slower convergence in the instanton search must be due to the more involved force scheme. On the other hand, the LBFGS optimizer finds instantons more efficiently, perhaps due to the smaller initial distance to the solution. An important conclusion is that, for simple systems, it is as easy to find instantons as it is to find MEPs.

5. TUNNELING SPLITTINGS IN MALONALDEHYDE IN FULL DIMENSIONALITY In this section, we test the methods introduced in the previous section on a three-dimensional application. The additional multidimensional effects include the overall system rotations and translations, as well as the treatment of “frozen” ends in the chain-of-states based methods. The efficiency of the methods is again measured in terms of the number of potential evaluations. The convergence criterion ϵ is defined in terms of the maximal action gradient at a bead. We compare the methods on a full-dimensional PES of malonaldehyde by Sewell et al.50 with the planar part of the potential modified as in ref 21. The potential is given in an analytical form, so that the gradients are obtained analytically, which makes it particularly convenient for testing purposes. The energy barrier along the MEP is 3497 cm−1. 5.1. Standard RPI Calculations. In the standard ringpolymer instanton calculations,25 the action S in eq 3 is minimized using LBFGS with respect to all 27 × N Cartesian coordinates. Overall rotations and translations do not interfere with the minimization.25 The linear path connecting two symmetric minima, with hydrogen atom at either oxygen, was used as the starting guess, and it was represented by N beads at equal Euclidean distances in mass-scaled Cartesian coordinates. 797


Article

Journal of Chemical Theory and Computation The coordinates of the end beads are connected by a reflection through a plane containing two principal axes of the malonaldehyde molecule. Following the action minimization, the prefactor Φ is obtained from Hessians evaluated at N converged bead positions. Reference 25 obtained the tunneling splitting Δ = 51 cm−1 using the standard RPI method. We repeated the ring-polymer instanton calculations of tunneling splittings in malonaldehyde for the two representative kink durations of β = 6000 and 10,000 a.u. and varied the number of beads N in steps of 10. Two sets of time steps di were used in eq 3: equal imaginary time steps di = β/N (α = 0) and the prescription from ref 28 with α = 0.5 in eq 11. The convergence of tunneling splittings Δ with the number of beads N is shown in Figure 11. Similar trends are at play as in 2D

Figure 12. Number of iterations to converge tunneling splittings Δ in the standard RPI calculations, shown in Figure 11, using tolerance ϵ = 10−5 a.u., and α = 0 and 0.5, in eq 11, for orbit durations β = 5000 a.u. (blue) and 8000 a.u. (red), as a function of number of beads N. Parameters are defined in the text. Initial guess in the optimizations is the converged path from the previous calculation with N− 10 beads.

we present below), for example, for β = 6000 a.u., the convergence of tunneling splittings to 51 cm−1 using α = 0 and 0.5 (and ϵ = 10−5) requires 310 × 558 = 172,980 and 250 × 472 = 118,000 potential evaluations, respectively. 5.2. Chain-of-States Based RPI Calculations. Next, we would like to apply the three methods that gave best results in terms of efficiency in the 2D calculations: S minimization with variable di values, SA minimization using NEB and string methods. In those methods, the end beads were fixed at PES minima in 2D calculations. In 3D molecular applications, it is only possible to fix the molecular geometry at PES minima. The relative orientation of molecules at beads with respect to its neighbors along the path remains to be adjusted in order to minimize the action in eq 15. The action gradient at the first and the last bead, see eq 17, points along the upwinding tangent and would be entirely projected out in the projection scheme that we use for the intermediate beads (eq 13). In thermal rate calculations, Einarsdóttir et al.31 project the action gradient at end beads along the potential isocontour by subtracting the component that lies along the potential gradient. At the PES minimum, the direction of the potential gradient is undefined. A potential resolution of the problem is to rotate the end beads according to the “torque” exerted by the adjacent bead. Instead, we simply reorient the end beads in the following way. The geometry of end bead 1 is fixed at the minimum, so p1 in eq 15 remains zero throughout the minimization. Dependence of SA on x1 is only through r2 = |x2 − x1|. After we have moved all the beads xi, with i = 2,...,N − 1, apart from end beads, in the current iteration, we readjust the end bead orientations so as to minimize the Euclidean distance in mass-scaled Cartesian coordinates to their nearest neighbor,

Figure 11. Convergence of tunneling splittings Δ with the number of beads N in a standard RPI calculation based on S minimization with LBFGS. Results are shown for the orbit durations β = 6000 and 10,000 a.u. and for the values of parameter α in eq 11 of 0 and 0.5. Gray line indicates the converged Δ.

results from the previous section, shown in Figure 2. For low N and for rough estimates of Δ, it is beneficial to use α = 0.5. For high accuracy results, the advantage of using α > 0 is lost. The precision of the results is difficult to judge as the diagonalization of the 27 × N-dimensional square matrix quickly becomes very costly. For β = 6000 a.u., the tunneling splittings round off to 51 cm−1 using 250 and 310 beads for α = 0 and 0.5, respectively. For β = 10,000 a.u., those bead numbers increase to 280 and 410, respectively. Using the larger orbit duration, β is justified only in high accuracy calculations (with large N). The number of iterations to converge Δ in LBFGS optimization as a function of the number of beads N for ϵ = 10−5 a.u. is shown in Figure 12. In the α = 0.5 calculation, the converged path from the previous calculation, with N−10 beads, is used as the starting guess. The trends in Figure 12 are similar to what we observed in 2D calculations in the previous section. The number of iterations increases with the number of beads used to describe the path. The calculations that use equal time steps (α = 0) require lower number of iterations to converge the instanton. The number of iterations for a particular value of β and N is related to the numerically calculated duration of the orbit using the quadrature with N points. At small N, β = 6000 a.u. leads to lower number of iterations, and at large N, the situation reverses. The calculations are computationally costly (compared to those

|x i + 1 − x i |min , for i = 1 and N − 1

(23)

This is a well-studied least-squares fit problem which has an analytic and elegant solution, with notable examples of Kabsch51 and Karney.52 In this work, we use the procedure by Karney,52 outlined in Appendix III, to reorient the ends beads toward their nearest neighbors, which solves eq 23. The reorientation of end beads is performed at every iteration in course of minimization. The rotation of “frozen” molecular end 798


Article


of Δ from our “best result”, below, we are studying the convergence of Δ on the shape of the path, that is, the convergence of the “tangent” error resulting from a finite number of beads N used to represent the path in chain-of-states based methods of our choice. Convergence of tunneling splittings within 3% of our best result is achieved using 18 beads and within 1% using 25/27 (for S/SA minimization) beads. The S minimization using variable di values results in orbit durations β between 5320 a.u. for 18 beads and 6730 a.u. for 40 beads. We also tried using cubic spline tangents instead of upwinding tangents and obtained the convergence within 1% using 18 beads. Using the string method, followed up by the hermite interpolation of path in terms of arc length for later use in the standard RPI optimization on a 1D path, instead of cubic splines, results in the tunneling splittings that are of similar quality but differ by ≈3.5% with N = 18 beads and by less than 0.1% for N > 24. The tangent error, originating from the finite-N representation of the path, is translated into the interpolation scheme. In Figure 14, we compare the number of iterations to converge Δ within ϵ using LBFGS for the three methods. Two

beads at PES minima is represented in terms of quaternions. The optimal quaternion is obtained by constructing and diagonalizing a 4-by-4 matrix, and it defines the rotation matrix that multiplies the coordinates of each atom of the rigid molecule at the minimum geometry to orient it toward the molecular geometry of the adjacent bead in the current guess of the instanton path. In this procedure, an approximation is made by ignoring the dependence of x1/N on x2/N−1 in the calculation of ∇2/N−1SA, which is difficult to obtain analytically. The convergence rate did not noticeably deteriorate as a result. Related problems with bead orientation in MEP optimizations using NEB have been discussed by Gonzalez-Garcia et al.53 and Bohner et al.47 Those problems persist in MAP optimizations based on the chain-of-states representation. The initial orientation of end beads is important in order to have a physically meaningful starting guess, when it is constructed by interpolation between end beads.53 Numerical contamination by the bead moves along directions that are projected out in our projection scheme,47 such as the relative bead orientations and translations, can potentially lead to errors in action, but the problem is less severe than in MEP optimizations due to the presence of straightening force in eq 17. In the present calculations that involve small number of beads, no related difficulties were observed. Dependence of tunneling splittings on the number of beads N in range [4−40] obtained by S minimization with variable di values and SA minimization using string method is shown in Figure 13. The NEB results are indistiguishable from the string

Figure 14. Number of iterations to converge tunneling splittings Δ using tolerance ϵ (defined in the text) vs number of beads N. Results obtained by SA minimization using NEB method (k = 0.03 a.u.) are shown using solid red lines. Red circles (lines) show results obtained by applying the convergence criterion to the perpendicular (full) NEB gradients. NEB results using the “optimal” spring constants are shown using solid dark red lines. The “optimal” spring constant in the NEB method is shown above the frame for each N. Results obtained by SA minimization using string method are shown using solid lines in shades of blue. Blue crosses are results obtained by S minimization with variable time steps.

Figure 13. Convergence of tunneling splittings Δ with the number of beads N. Results using SA minimization with NEB and string methods are shown using blue line (indistiguishable). Results obtained by S minimization with variable time steps di are shown with red crosses. Solid, dashed, and dotted black lines show the exact instanton result and 1% and 3% deviation limits, respectively.

sets of NEB results, shown as lines, are obtained using the spring constant k = 0.03 a.u. and using the optimal spring constant determined in the same way as in 2D calculations above (shown in Figure 9). The value of the optimal spring constant is shown in the bar graph in Figure 14. The convergence criterion ϵ pertains to the full NEB gradient (including the spring force) in results shown by lines. The results shown by circles are obtained using the perpendicular NEB gradient in ϵ. All three methods work with similar efficiency. For 18 beads, we need 117, 105, and 97 iterations using S minimization, NEB and string method, respectively, which leads to 1872, 1680, and

results (and are not shown). A standard RPI calculation on the 1D interpolated potential was performed using α = 0.5, β = 10,000 a.u., and N = 600 for each converged N-bead path. This is not sufficient for full convergence of Δ. Our estimate is that the results obtained using this choice of parameters are within ≈3% of the exact instanton result. Using larger β and N would have been too time consuming, as the calculation of Φ in eq 5 involves a diagonalization of a 27 × N-dimensional banded SHessian matrix. Our “best result” here is obtained by the standard RPI with the above choice of parameters on the converged 1D instanton path. By measuring the discrepancies 799


Article


times steps lead to failure of the method (for some N), whereas smaller time steps lead to a less efficient search. Instanton search using FIRE algorithm48 with a variable time step size proved unstable, presumably due to same reasons as in the 2D model above. The MEP search using FIRE algorithm, where the criteria for changing the time step size were applied to each bead separately, proved to be stable, if we limited the size of the time step to Δtmax = 100 a.u. The number of iterations to reach convergence did not vary significantly with the size of the initial time step in results obtained using FIRE. On the contrary, the efficiency of MEP search using QM varies drastically with the initial time step size. Best results were achieved using Δt = 80− 90 a.u., and the efficiency savings are then similar to FIRE results shown in Figure 15. Since it is difficult to predict the optimal step size in advance, FIRE is preferred over QM in MEP searches. Again, as in Figure 10 (2D model), we conclude that LBFGS outperforms the MD-based optimizers for locating simple paths (MEP and instanton) in malonaldehyde. For damped MDbased optimizers, the instanton search is more demanding than the MEP search, probably due to the more complicated force scheme. When using LBFGS, the situation reverses and the instanton search is slightly more efficient, probably because the path has a shorter distance to travel to the solution. 5.3. Tunneling Splittings on a Realistic Potential. Finally, we apply the instanton method to a recent realistic potential by Wang et al.10 The PES is a weighted least-squares fit of 11,147 near basis-set-limit frozen-core CCSD(T) electronic energies to the permutationally invariant polynomials in the Morse-type functions of internuclear distances. The barrier at transition state is 1438 cm−1. The tunneling splittings on this potential have been obtained using MCTDH, which is formally exact. The values from two different sources54,55 are 23.8 and 23.4 cm−1. Viel et al.56 obtained the tunneling splittings of 25.7 and 3.21 cm−1 for H and D tunneling atoms in malonaldehyde by means of an exact method, respectively, but using a different potential energy surface.57 The experimental values for H and D tunneling splittings in malonaldehyde are 21.6 cm−158,59 and 2.92 cm−1,60 respectively. The fitted PES by Wang et al. is inherently bumpy, which can give rise to problems during optimization based on gradients. We calculated gradients and Hessians numerically using a step of 10−3 a0. With low quality potentials, the methods that involve fewer beads for representing the path, such as the string and NEB methods we use here, potentially have an advantage as they are less likely to encounter rough regions of PES that contain artifacts of the fit/interpolant. We used the string method with the linear starting path (of length 72 m1/2 e a0 or 1.47 a0 without mass scaling) connecting two symmetric minima. We needed only nine beads to converge the tunneling splittings to within 2%; 14 beads were needed for convergence to within 1%. The instanton path was ≈1.60 a0 long for both H and D tunneling (104 and 129 m1/2 e a0 using mass-scaled distances, respectively), which means that the path was almost straight. The r.m.s. distance, defined in eq 19, 1/2 was 24 m1/2 e a0 (28 me a0) for H (D) tunneling path. The number of iterations to converge the H (D) tunneling path within ϵ = 10−5 a.u. was 60 (64) and 86 (78), similar to the 27D model considered above. Total number of potential evaluations to converge the H tunneling splittings in malonaldehyde within 1% was (14 − 2) × 85 = 1020. We used β = 10,000 a.u. and 1200 beads in the subsequent RPI

1552 potential evaluations. The potential evaluations at 16 beads are independent and can be computed in parallel. For 25 beads, the number of iterations is, respectively, 143, 138 and 147, leading to 3289, 3174, and 3381 potential evalulations. We note here that the NEB and string methods work equally well here, although the NEB method was outperformed by the string method in the 2D model of malonaldehyde, as shown in Figure 9. We still prefer to use the string method because it is parameter-free. We also tried including a perpendicular spring force with a switching function that turns it off near the solution in the way it was done in ref 32 (eqs 9 and 10) and ref 46 (eqs 14 and 15). The purpose of introducing the additional spring force proved to be necessary to stabilize the NEB and string methods in some more involved applications.42 It is not needed for the simple path such as the tunneling path in malonaldehyde, but we learned from the test that the introduction of the additional straightening force does not change the rate of convergence in terms of the number of iterations. The full-dimensional calculations presented here are almost 5-fold computationally more expensive than in the 2D model potential above but outperform (in terms of the number of potential evaluations) the standard RPI method with fixed β and di values by almost 2 orders of magnitude. Next, we use the string method to compare the damped MDbased optimizers to LBFGS optimizer in the instanton and MEP searches in the full-dimensional malonaldehyde. The r.m.s. distances of the linear starting paths to the instanton and MEP paths are 19 and 69 m1/2 e a0, respectively, whereas the corresponding path lengths are 100 and 251 m e1/2 a 0 , respectively. The convergence criterion ϵ is modified for the purpose of comparing instanton and MEP results in the same manner as in the 2D model above (shown in Figure 10). ϵ is defined as the average Euclidean distance of a bead to the corresponding bead on the fully converged path for a given N. Results are shown in Figure 15. Instanton search using the QuickMin algorithm46 with a fixed time step of Δt = 3 a.u. proved to be stable and efficient. Larger

Figure 15. Number of iterations needed to converge tunneling splittings and MEPs using string method and tolerance ϵ, based on the Euclidean mass-scaled distance to the exact N-beads path, as defined in the text, as a function of number of beads N. Solid lines refer to instantons; dashed lines refer to MEPs. Results shown in shades of red are obtained using LBFGS optimizer; results shown in shades of blue are obtained using QM (for instanton path, using Δt = 3 a.u.) or FIRE (for MEPs). 800


Article


with expensive potential energy surfaces and on-the-fly electronic structure calculations.

calculation on the 1D interpolated potential and obtained the ground state tunneling splittings of 25 and 3.4 cm−1 for tunneling splittings of hydrogen and deuterium in malonaldehyde, respectively. Our instanton result is ≈5% off from results obtained by MCTDH.54,55 The errors in tunneling splittings obtained by the instanton method are similar to the current discrepancies between results obtained using different ab initio potentials and also between theory and experiment, but are, with ≈1,000 potential evaluations, much less costly to obtain.

■

APPENDIX I

Action Gradient and Hessian

The methods for locating instantons in this paper are based on the minimization of action in eq 3. The minimization algorithms that we use here require gradients with respect to intermediate bead positions xi, with i = 2,...,N − 1. The gradient components are given by xi , j − xi − 1, j xi , j − xi + 1, j ∂V (x i) di + di + 1 ∂S kink = + + ∂xi , j di di + 1 ∂xi , j 2

6. CONCLUSIONS We have proposed three efficient algorithms for locating instanton paths in calculations of tunneling splittings. We concentrated on methods that represent the tunneling path by a string of beads in Cartesian coordinates, where each bead represents the molecular geometry along the tunneling path. The three methods represent the tunneling path by equally spaced beads in terms of Euclidean distance in mass-scaled coordinates. They also make use of the knowledge of molecular geometry at PES minimum and use the LBFGS optimizer to iteratively improve on the initial guess until convergence is reached. Tunneling splittings are obtained from the converged tunneling path by a standard RPI calculation on a 1D interpolated potential along the path using a dense imaginary time grid. Potential and Hessians on a dense imaginary time grid are obtained by interpolation. The method based on Euclidean action minimization with variable imaginary time steps is a modification of the method by Rommel and Kästner28 for rate calculations. The ends beads are placed at PES minima, and the action is minimized using LBFGS. At every iteration, the beads are redistributed equidistantly by interpolation, and the time steps, as well as the duration of the orbit, are determined by the requirement that the action gradient points perpendicular to the path. In this way, duration of the orbit is not a convergence parameter, and beads do not accumulate near PES minima. The other two algorithms are based on the string method and NEB method, introduced originally for finding MEPs. The end beads are again placed at two symmetric minima equidistantly, and abbreviated action in eq 15 is minimized. The perpendicular action gradient is used to evolve string/NEB on the PES, instead of the perpendicular potential gradient. Additional nudging forces that improve stability of the calculations for more involved paths can easily be added to both methods. The NEB method requires the spring constant as a parameter, whereas the string method is parameter free. The crucial step in applying all three methods on a realistic 3D molecular system was to reorient the end beads toward its nearest neighbor at every iteration using quaternions to solve eq 23. The methods can easily be parallelized by distributing gradient evaluations at each bead along the path to a different processor at every iteration. The methods were successfully applied to malonaldehyde. All three methods have similar computational cost and require almost 2 orders of magnitude fewer potential evalulations than a standard RPI method. Overall, we recommend using the string method to minimize the reduced action since this approach is parameter-free and does not require computation of time steps. The cost of finding the tunneling path in terms of the number of potential evaluations was similar to the cost of finding a minimum-energy path. These approaches should thus enable the RPI method to be used effectively in combination

(24)

where j labels Cartesian coordinates of molecular system. Once the MAP has been located, the tunneling splittings are calculated using eq 6, which contains the prefactor Φ. One way to calculate Φ in eq 5 uses eigenvalues of the S-hessian matrix. The S-hessian matrix elements are given by ⎤ ⎡⎛ 1 ∂ 2S kink 1 ⎞ 1 1 = ⎢⎜ + δi + 1, i ′⎥δj , j ′ ⎟δi , i ′ − δi , i ′− 1 − ⎥⎦ di + 1 ⎠ di di + 1 ∂xi , j∂xi ′ , j ′ ⎣⎢⎝ di +

∂ 2V (x i) di + di + 1 δii ′ 2 ∂xi , j∂xi ′ , j ′

(25)

for 1 < i < N and 1 < j < f.

■

APPENDIX II

Upwinding Tangents for NEB and String Methods

The NEB and string methods require specification of tangents τi at every bead xi, i = 1,...,N, along the path, which are used in the force projection scheme, as described in Section 3.2.1. We use, unless otherwise specified, the so-called upwinding tangents, defined in ref 45 as ⎧ τi+ if Vi + 1 > Vi > Vi − 1 ⎪ ⎪ τi− if Vi + 1 < Vi < Vi − 1 ⎪ − min ⎪ + max τi = ⎨ τi ΔVi + τi ΔVi if Vi + 1 > Vi − 1 > Vi ⎪ or Vi > Vi + 1 > Vi − 1 ⎪ + − min max ⎪ τi ΔVi + τi ΔVi if Vi + 1 < Vi − 1 < Vi ⎪ or Vi < Vi + 1 < Vi − 1 ⎩

where Vi = V(xi), and

τ+i

= xi+1 −

xi, τ−i

= xi − xi−1, τ1 =

(26)

τ+1 ,

τN = τ−N,

ΔVimax = max(|Vi + 1 − Vi |, |Vi − 1 − Vi |) ΔVimin = min(|Vi + 1 − Vi |, |Vi − 1 − Vi |)

(27)

The unit tangent vectors τ̂i in eqs 12 and 13 are obtained by normalizing τi in eq 26.

■

APPENDIX III

Least-Squares Fit of Rotations

Let ai be the position vector of atom i of the end bead x1/N with Cartesian components ai , j = x1,3(i − 1) + j/ mi , where j = 1−3 and i = 1−n, where n is the number of atoms in the molecule (of f = 3n degrees of freedom). The position vector bi is defined likewise as the position vector (without mass scaling) of the adjacent bead x2/N−1. We assume that both molecules (beads 801


Article


(4) Richardson, J. O.; Wales, D. J.; Althorpe, S. C.; McLaughlin, R. P.; Viant, M. R.; Shih, O.; Saykally, R. J. J. Phys. Chem. A 2013, 117, 6960−6966. (5) Althorpe, S. C.; Clary, D. C. J. Chem. Phys. 1994, 101, 3603− 3609. (6) Hammer, T.; Coutinho-Neto, M. D.; Viel, A.; Manthe, U. J. Chem. Phys. 2009, 131, 224109. (7) Neff, M.; Rauhut, G. Spectrochim. Acta, Part A 2014, 119, 100− 106. (8) Marquardt, R.; Sagui, K.; Zheng, J.; Thiel, W.; Luckhaus, D.; Yurchenko, S.; Mariotti, F.; Quack, M. J. Phys. Chem. A 2013, 117, 7502−7522. (9) Coker, D. F.; Watts, R. O. J. Phys. Chem. 1987, 91, 2513−2518. (10) Wang, Y.; Braams, B. J.; Bowman, J. M.; Carter, S.; Tew, D. P. J. Chem. Phys. 2008, 128, 224314. (11) Ceperley, D. M.; Jacucci, G. Phys. Rev. Lett. 1987, 58, 1648− 1651. (12) Marchi, M.; Chandler, D. J. Chem. Phys. 1991, 95, 889−894. (13) Miller, W. H. J. Phys. Chem. 1979, 83, 960−963. (14) Makri, M.; Miller, W. H. J. Chem. Phys. 1989, 91, 4026−4036. (15) Ben-Nun, M.; Martinez, T. J. J. Phys. Chem. A 1999, 103, 6055− 6059. (16) Wales, D. J. J. Am. Chem. Soc. 1993, 115, 11191−11201. (17) Takahashi, M.; Watanabe, Y.; Taketsugu, T.; Wales, D. J. J. Chem. Phys. 2005, 123, 044302. (18) Smedarchina, Z.; Siebrand, W.; Fernandez-Ramos, A. J. Chem. Phys. 2012, 137, 224105. (19) Callan, C. G.; Coleman, S. Phys. Rev. D: Part. Fields 1977, 16, 1762−1768. (20) Vainshtein, A. I.; Zakharov, V. I.; Novikov, V. A.; Shifman, M. A. Sov. Phys. Usp. 1982, 25, 195−215. (21) Mil’nikov, G. V.; Nakamura, H. J. Chem. Phys. 2001, 115, 6881− 6897. (22) Nakamura, H.; Mil′nikov, G. V. Quantum Mechanical Tunneling in Chemical Physics; CRC Press, Inc, 2013. (23) Mil’nikov, G. V.; Yagi, K.; Taketsugu, T.; Nakamura, H.; Hirao, K. J. Chem. Phys. 2004, 120, 5036−5045. (24) Richardson, J. O.; Althorpe, S. C. J. Chem. Phys. 2009, 131, 214106. (25) Richardson, J. O.; Althorpe, S. C. J. Chem. Phys. 2011, 134, 054109. (26) Richardson, J. O.; Althorpe, S. C.; Wales, D. J. J. Chem. Phys. 2011, 135, 124109. (27) Rommel, J. B.; Goumans, T. P. M.; Kästner, J. J. Chem. Theory Comput. 2011, 7, 690−698. (28) Rommel, J. B.; Kästner, J. J. Chem. Phys. 2011, 134, 184107. (29) Kryvohuz, M. J. Chem. Phys. 2011, 134, 114103. (30) Kryvohuz, M. J. Chem. Phys. 2012, 137, 234304. (31) Einarsdóttir, D. M.; Arnaldsson, A.; Ó skarsson, F.; Jónsson, H. Path optimization with application to tunneling. Lecture Notes in Computer Science 2012, 7134, 45−55. (32) Jónsson, H.; Mills, G.; Jacobsen, K. W. Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions. In Classical and Quantum Dynamics in Condensed Phase Simulations; World Scientific, 1998; pp 385−404. (33) Kawatsu, T.; Miura, S. J. Chem. Phys. 2014, 141, 024101. (34) Miller, W. H. J. Chem. Phys. 1975, 62, 1899−1906. (35) Gutzwiller, M. C. J. Math. Phys. 1971, 12, 343−358. (36) Goldstein, H. Classical Mechanics; Addison-Wesley, 1950; pp 365−371. (37) Olsen, R. A.; Kroes, G. J.; Henkelman, G.; Arnaldsson, A.; Jónsson, H. J. Chem. Phys. 2004, 121, 9776−9792. (38) Henkelman, G.; Jónsson, H. J. Chem. Phys. 1999, 111, 7010− 7022. (39) Zhu, C.; Byrd, R. H.; Lu, P.; Nocedal, J. ACM Trans. Math. Softw. 1997, 23, 550−560. (40) Banerjee, A.; Adams, N.; Simons, J.; Shepard, R. J. Phys. Chem. 1985, 89, 52−57. (41) Bofill, J. M. J. Comput. Chem. 1994, 15, 1−11.

x1/N) have center of mass at the origin. If this is not the case, the translations can be added/subtracted trivially. The condition given by eq 23 transcribes to 1 M

n

∑ mi |bi − R(a i)| = min i=1

(28)

where M = ∑imi is the total mass of the molecule. We wish to determine the rotation R that minimizes eq 28 and, in that, we follow ref 52. Alternatives that avoid description of rotations in terms of quaternions are available51 and lead to equivalent results. We define 4 × 4 skew matrices Ai = A(bi + ai,bi − ai), with ⎛ 0 −f −f −f ⎞ 1 2 3 ⎜ ⎟ ⎜ f 0 − e3 e 2 ⎟ 1 ⎟ A(e, f) = ⎜ ⎜ f e3 ⎟ − 0 e 1 ⎜ 2 ⎟ ⎜ f −e e 0 ⎟⎠ 2 1 ⎝ 3

(29)

where e and f are spatial vectors with Cartesian components ej and f j, respectively. Then, the symmetric matrix B is formed as B=

1 M

n

∑ mi A iTA i i=1

(30)

and diagonalized. All eigenvalues are real and positive, and the smallest eigenvalue gives the minimum value of the expression in eq 28. The optimal rotation is obtained by setting the quaternion q = (q1,q2,q3,q4) to the corresponding eigenvector of B. The rotation matrix defined by the quaternion q is ⎛q 2 + q 2 − q 2 − q 2 2(q2q3 − q1q4) 2(q1q3 + q2q4)⎞ 2 3 4 ⎟ ⎜ 1 ⎟ ⎜ 2 2 2 2 q1 − q2 + q3 − q4 2(q3q4 − q1q2)⎟ 2(q1q4 + q2q3) R q = ⎜⎜ ⎟ 2(q1q2 + q3q4) q12 − q22 − q32 ⎟ ⎜ 2(q2q4 − q1q3) ⎟ ⎜ + q42 ⎠ ⎝

(31)

The new optimized Cartesian coordinates of the positions of end bead atoms aopt are obtained by a rotation using Rq, i namely, aopt = R(ai) = Rq ai. i

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work has been supported in part by Croatian Science Foundation under project IP-2014-09-7540. The authors acknowledge the assistance of the COST Action CM1405 “Molecules in Motion” (MOLIM).

■

REFERENCES

(1) DeVault, D. Quantum Mechanical Tunneling in Biological Systems. Q. Rev. Biophys. 1980, 13, 387−564. (2) Benderskii, V. A.; Makarov, D. E.; Wight, C. A. Chemical Dynamics at Low Temperatures; Advances in Chemical Physics Series; John Wiley & Sons, Inc.: Hoboken, NJ, 1994. (3) Keutsch, F. N.; Saykally, R. J. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10533−10540. 802


Article

Journal of Chemical Theory and Computation (42) Trygubenko, S. A.; Wales, D. J. J. Chem. Phys. 2004, 120, 2082− 2094. (43) E, W.; Ren, W.; Vanden-Eijnden, E. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 66, 052301. (44) E, W.; Ren, W.; Vanden-Eijnden, E. J. Chem. Phys. 2007, 126, 164103. (45) Henkelman, G.; Jónsson, H. J. Chem. Phys. 2000, 113, 9978− 9985. (46) Sheppard, D.; Terrell, R.; Henkelman, G. J. Chem. Phys. 2008, 128, 134106. (47) Bohner, M. U.; Meisner, J.; Kästner, J. J. Chem. Theory Comput. 2013, 9, 3498−3504. (48) Bitzek, E.; Koskinen, P.; Gähler, F.; Moseler, M.; Gumbsch, P. Phys. Rev. Lett. 2006, 97, 170201. (49) Bosch, E.; Moreno, M.; Lluch, J. M.; Bertrán, J. J. Chem. Phys. 1990, 93, 5685−5692. (50) Sewell, T. D.; Guo, Y.; Thompson, D. L. J. Chem. Phys. 1995, 103, 8557−8565. (51) Kabsch, W. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1976, 32, 922−923. Kabsch, W. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1978, 34, 827−828. (52) Karney, C. F. F. J. Mol. Graphics Modell. 2007, 25, 595−604. (53) Gonzalez-Garcia, N.; Pu, J.; Gonzalez-Lafont, A.; Lluch, J. M.; Truhlar, D. G. J. Chem. Theory Comput. 2006, 2, 895−904. (54) Schröder, M.; Gatti, F.; Meyer, H.-D. J. Chem. Phys. 2011, 134, 234307. (55) Hammer, T.; Manthe, U. J. Chem. Phys. 2011, 134, 224305. (56) Viel, A.; Coutinho-Neto, M. D.; Manthe, U. J. Chem. Phys. 2007, 126, 024308. (57) Yagi, K.; Taketsugu, T.; Hirao, K. J. Chem. Phys. 2001, 115, 10647−10655. (58) Firth, D. W.; Beyer, K.; Dvorak, M. A.; Reeve, S. W.; Grushow, A.; Leopold, K. R. J. Chem. Phys. 1991, 94, 1812−1819. (59) Baba, T.; Tanaka, T.; Morino, I.; Yamada, K. M. T.; Tanaka, K. J. Chem. Phys. 1999, 110, 4131−4133. (60) Baughcum, S. L.; Duerst, R. W.; Rowe, W. F.; Smith, Z.; Wilson, E. B. J. Am. Chem. Soc. 1981, 103, 6296−6303.

803


Locating Instantons in Calculations of Tunneling Splittings: The Test

Recommend Documents