Should Thermostatted Ring Polymer Molecular Dynamics be used to calculate thermal reaction rates?
SShould Thermostatted Ring Polymer Molecular Dynamics be used tocalculate thermal reaction rates?
Timothy J. H. Hele and Yury V. Suleimanov
2, 3 Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW,UK. Computation-based Science and Technology Research Center, Cyprus Institute, 20 Kavafi Str., Nicosia 2121,Cyprus Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Ave., Cambridge,Massachusetts 02139, United States (Dated: 18 September 2018)
We apply Thermostatted Ring Polymer Molecular Dynamics (TRPMD), a recently-proposed approximatequantum dynamics method, to the computation of thermal reaction rates. Its short-time Transition-StateTheory (TST) limit is identical to rigorous Quantum Transition-State Theory, and we find that its long-time limit is independent of the location of the dividing surface. TRPMD rate theory is then applied toone-dimensional model systems, the atom-diatom bimolecular reactions H+H , D+MuH and F+H , and theprototypical polyatomic reaction H+CH . Above the crossover temperature, the TRPMD rate is virtuallyinvariant to the strength of the friction applied to the internal ring-polymer normal modes, and beneath thecrossover temperature the TRPMD rate generally decreases with increasing friction, in agreement with thepredictions of Kramers theory. We therefore find that TRPMD is approximately equal to, or less accuratethan, Ring Polymer Molecular Dynamics (RPMD) for symmetric reactions, and for certain asymmetric sys-tems and friction parameters closer to the quantum result, providing a basis for further assessment of theaccuracy of this method. Copyright (2015) American Institute of Physics. This article may be downloadedfor personal use only. Any other use requires prior permission of the author and the American Instituteof Physics. The following article appeared in J. Chem. Phys., (2015), 074107, and may be found athttp://dx.doi.org/10.1063/1.4928599.
I. INTRODUCTION
The accurate computation of thermal quantum ratesis a major challenge in theoretical chemistry, as a purelyclassical description of the kinetics fails to capture zero-point energy, tunnelling, and phase effects . Exactsolutions using correlation functions, developed by Ya-mamoto, Miller, and others are only tractable forsmall or model systems, as the difficulty of computationscales exponentially with the size of the system.Consequently, numerous approximate treatments havebeen developed, which can be broadly classed as thoseseeking an accurate description of the quantum statisticswithout direct calculation of the dynamics, and thosewhich also seek to use an approximate quantum dy-namics. Methods in the first category include instan-ton theory , “quantum instanton” , and varioustransition-state theory (TST) approaches . Of manyapproximate quantum dynamics methods, particularlysuccessful ones include the linearized semiclassical initial-value representation (LSC-IVR) , centroid moleculardynamics (CMD) , and ring polymer molecular dy-namics (RPMD) .RPMD has been very successful for the computationof thermal quantum rates in condensed-phase processes,due to the possibility of implementation in complex sys-tems such as (proton-coupled) electron transfer reac-tion dynamics or enzyme catalysis, and especially insmall gas-phase systems where comparison withexact quantum rates and experimental data has demon- strated that RPMD rate theory is a consistent and re-liable approach with a high level of accuracy. Thesenumerical results have shown that RPMD rate theoryis exact in the high-temperature limit (which can alsobe shown algebraically ), reliable at intermediate tem-peratures, and more accurate than other approximatemethods in the deep tunnelling regime (see Eq. (24)below), where it is within a factor of 2–3 of the exactquantum result. RPMD also captures zero-point energyeffects, and provides very accurate estimates for barri-erless reactions . It has been found to systematicallyoverestimate thermal rates for asymmetric reactions andunderestimate them for symmetric (and quasisymmet-ric) reactions in the deep tunnelling regime (Note thatzero-point energy effects along the reaction coordinatemust be taken into account when assigning the reactionsymmetry.) Recently a general code for RPMD cal-culations (RPMDrate) has been developed. Another appealing feature of RPMD rate theory isits rigorous independence to the location of the divid-ing surface between products and reactants , a prop-erty shared by classical rate theory and the exact quan-tum rate , but not by many transition-state theory ap-proaches. The t → + , TST limit of RPMD (RPMD-TST) is identical to true QTST: the instantaneous ther-mal quantum flux through a position-space dividing sur-face which is equal to the exact quantum rate in theabsence of recrossing . A corollary of this is thatRPMD will be exact for a parabolic barrier (where thereis no recrossing of the optimal dividing surface by RPMD a r X i v : . [ phy s i c s . c h e m - ph ] A ug dynamics or quantum dynamics, and QTST is thereforealso exact) .When the centroid is used as the dividing surface(see Eq. (25) below), RPMD-TST reduces to the ear-lier theory of centroid-TST , which is a goodapproximation for symmetric barriers but significantlyoverestimates the rate for asymmetric barriers at lowtemperatures . This effect is attributable to thecentroid being a poor dividing surface beneath the‘crossover’ temperature into deep tunnelling . In this‘deep tunnelling’ regime, RPMD-TST has a close rela-tionship to semiclassical “Im F” instanton theory ,which has been very successful for calculating rates be-neath the crossover temperature, though has no first-principles derivation and was recently shown to be lessaccurate than QTST when applied to realistic multidi-mensional reactions .Very recently, both CMD and RPMD have beenobtained from the exact quantum Kubo-transformed time-correlation function (with explicit error terms) via aBoltzmann-conserving “Matsubara dynamics” whichconsiders evolution of the low-frequency, smooth “Mat-subara” modes of the path integral . Matsubara dynam-ics suffers from the sign problem and is not presentlyamenable to computation on large systems. However,by taking a mean-field approximation to the centroiddynamics, such that fluctuations around the centroidare discarded, one obtains CMD. Alternatively, if themomentum contour is moved into the complex planein order to make the quantum Boltzmann distributionreal, a complex Liouvillian arises, the imaginary partof which only affects the higher, non-centroid, normalmodes. Discarding the imaginary Liouvillian leads tospurious springs in the dynamics and gives RPMD. Consequently, RPMD will be a reasonable approximationto Matsubara dynamics, provided that the timescale overwhich the resultant dynamics is required (the timescaleof ‘falling off’ the barrier in rate theory) is shorter thanthe timescale over which the springs ‘contaminate’ thedynamics of interest (in rate theory, this is usually cou-pling of the springs in the higher normal modes to themotion of the centroid dividing surface via anharmonicityin the potential).Both RPMD and CMD are inaccurate for the compu-tation of multidimensional spectra: the neglect of fluctu-ations in CMD leads to the “curvature problem” wherethe spectrum is red-shifted and broadened, whereas inRPMD the springs couple to the external potential lead-ing to “spurious resonances” . Recently, this problemhas been solved by attaching a Langevin thermostat to the internal modes of the ring polymer (which hadpreviously been used for the computation of statisticalproperties ), and the resulting Thermostatted RPMD(TRPMD) had neither the curvature nor resonance prob-lem.The success of RPMD for rate calculation, and theattachment of a thermostat for improving its compu-tation of spectra, naturally motivates studying whether TRPMD will be superior for the computation of ther-mal quantum rates to RPMD (and other approximatetheories) , which this article investigates. Given thatRPMD is one of the most accurate approximate meth-ods for systems where the quantum rates are availablefor comparison, further improvements would be of con-siderable benefit to the field.We firstly review TRPMD dynamics in section II A,followed by developing TRPMD rate theory in sec-tion II B. To predict the behaviour of the RPMD ratecompared to the TRPMD rate, we apply one-dimensionalKramers theory to the ring-polymer potential energysurface in section II C. Numerical results in section IIIapply TRPMD to the symmetric and asymmetric Eckartbarriers followed by representative bimolecular reac-tions: H+H (symmetric), D+MuH (quasisymmetrical),H+CH (prototypical polyatomic reaction) and F+H (asymmetric and highly anharmonic). Conclusions andavenues for further research are presented in section IV. II. THEORYA. Thermostatted Ring Polymer Molecular Dynamics
For simplicity we consider a one-dimensional system( F = 1) with position q and associated momentum p at inverse temperature β = 1 /k B T , where the N -beadring-polymer Hamiltonian is H N ( p , q ) = N − (cid:88) i =0 p i m + U N ( q ) (1)with the ring-polymer potential U N ( q ) = N − (cid:88) i =0 12 mω N ( q i − q i − ) + V ( q i ) (2)and the frequency of the ring-polymer springs ω N =1 /β N (cid:126) , where β N ≡ β/N . Generalization to further di-mensions follows immediately, and merely requires moreindices. The ring polymer is time-evolved by propagatingstochastic trajectories using TRPMD dynamics ,˙ p = − ∇ q U N ( q ) − Γp + (cid:115) m Γ β N ξ ( t ) (3)˙ q = 1 m p (4)where q ≡ ( q , . . . , q N − ) is the vector of bead positionsand p the vector of bead momenta, with ∇ q the gradoperator in position-space, ξ ( t ) a vector of N uniformGaussian deviates with zero mean and unit variance, and Γ the N × N positive semi-definite friction matrix .The Fokker-Planck operator corresponding to theTRPMD dynamics in Eqs. (3) and (4) is A N = − p m · ∇ q + U N ( q ) ←−∇ q · −→∇ p + ∇ p · Γ · p + mβ N ∇ p · Γ · ∇ p (5)(where the arrows correspond to the direction in whichthe derivative acts ) and for any Γ , TRPMD dynam-ics will conserve the quantum Boltzmann distribution( A N e − β N H N ( p , q ) = 0), a feature shared by RPMD andCMD but not some other approximate methods such asLSC-IVR . We then show in appendix A thatTRPMD obeys detailed balance, such that the TRPMDcorrelation function is invariant to swapping the opera-tors at zero time and finite time, and changing the signof the momenta.The time-evolution of an observable is given by theadjoint of Eq. (5), A † N = p m · ∇ q − U N ( q ) ←−∇ q · −→∇ p − p · Γ · ∇ p + mβ N ∇ p · Γ · ∇ p . (6)In the zero-friction limit, Γ = and A † N = L † N , where L N † is the adjoint of the Liouvillian corresponding todeterministic ring-polymer trajectories . B. TRPMD rate theory
We assume the standard depiction of rate dynamics,with a thermal distribution of reactants and a dividingsurface in position space. In what follows we assumescattering dynamics, with the potential tending to a con-stant value at large separation of products and reactants.The methodology is then immediately applicable to con-densed phase systems subject to the usual caveat thatthere is sufficient separation of timescales between reac-tion and equilibration.
The exact quantum rate can be formally given asthe long-time limit of the flux-side time-correlationfunction k QM ( β ) = lim t →∞ c QMfs ( t ) Q r ( β ) (7)where Q r ( β ) is the partition function in the reactant re-gion and c QMfs ( t ) = 1 β (cid:90) β dσ Tr (cid:104) e − ( β − σ ) ˆ H ˆ F e − σ ˆ H e i ˆ Ht/ (cid:126) ˆ he − i ˆ Ht/ (cid:126) (cid:105) (8)with ˆ F and ˆ h the quantum flux and side operators re-spectively, and ˆ H the Hamiltonian for the system. Thequantum rate can equivalently be given as minus the long-time limit of the time-derivative of the side-side cor-relation function, or the integral over the flux-flux corre-lation function .The TRPMD side-side correlation function is C TRPMDss ( t ) = 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q × e − β N H N ( p , q ) h [ f ( q )] h [ f ( q t )] (9)where (cid:82) d q ≡ (cid:82) ∞−∞ dq (cid:82) ∞−∞ dq . . . (cid:82) ∞−∞ dq N − and like-wise for (cid:82) d p , and q t ≡ q t ( p , q , t ) is obtained by evolu-tion of ( p , q ) for time t with TRPMD dynamics. The ringpolymer reaction co-ordinate f ( q ) is defined such thatthe dividing surface is at f ( q ) = 0, and that f ( q ) > f ( q ) < C TRPMDfs ( t ) = − ddt C TRPMDss ( t ) (10)= 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − β N H N ( p , q ) × δ [ f ( q )] S N ( p , q ) h [ f ( q t )] (11)where S N ( p , q ) is the flux perpendicular to f ( q ) at time t = 0, S N ( p , q ) = N − (cid:88) i =0 ∂f ( q ) ∂q i p i m . (12)We approximate the long-time limit of the quantum flux-side time-correlation function in Eq. (8) as the long-timelimit of the TRPMD flux-side time-correlation functionin Eq. (11), leading to the TRPMD approximation to thequantum rate as k TRPMD ( β ) = lim t →∞ C TRPMDfs ( t ) Q r ( β ) . (13)The flux-side time-correlation function Eq. (11) willdecay from an initial TST ( t → + ) value to a plateau,which (for a gas-phase scattering reaction with no frictionon motion out of the reactant or product channel) willextend to infinity. For condensed-phase reactions (andgas-phase reactions with friction in exit channels) a rateis defined provided that there is sufficient separation oftimescales between reaction and equilibration to define aplateau in C TRPMDfs ( t ), which at very long times (of theorder k TRPMD ( β ) − for a unimolecular reaction) tends tozero .Further differentiation of the flux-side time-correlationfunction (with the adjoint of the Fokker-Planck opera-tor in Eq. (6)) yields the TRPMD flux-flux correlationfunction C TRPMDff ( t ) = 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − β N H N ( p , q ) × δ [ f ( q )] S N ( p , q ) δ [ f ( q t )] S N ( p t , q t )(14)which, by construction, must be zero in the plateau re-gion, during which no trajectories recross the dividingsurface.Like RPMD rate theory, TRPMD has the appealingfeature that its short-time (TST) limit is identical to trueQuantum Transition-State Theory (QTST), as can beobserved by applying the short-time limit of the Fokker-Planck propagator e A † N t to f ( q ), yielding lim t → + C TRPMDfs ( t ) Q r ( β ) = k ‡ QM ( β ) (15)where k ‡ QM ( β ) is the QTST rate . In Appendix Bwe then show that the TRPMD rate in Eq. (13) is rig-orously independent of the location of the dividing sur-face. Consequently, the TRPMD rate will equal the exactquantum rate in the absence of recrossing of the opti-mal dividing surface (and those orthogonal to it in path-integral space) by either the exact quantum or TRPMDdynamics. We also note that Eq. (15) holds regardlessof the value of the friction matrix Γ and that recrossingof individual (stochastic) trajectories can only reduce theTRPMD rate from the QTST value, and hence QTST isan upper bound to the long-time TRPMD rate.In the following calculations we use a friction matrixwhich corresponds to damping of the free ring polymervibrational frequencies, and which has been used in previ-ous studies of TRPMD for spectra. For an orthogonaltransformation matrix T such that T T KT = m Ω (16)where K is the spring matrix in Eq. (2) and Ω ij = 2 δ ij sin( jπ/N ) /β N (cid:126) , the friction matrix is given by Γ = 2 λ TΩT T . (17)Here λ is an adjustable parameter, with λ = 1 giving crit-ical damping of the free ring polymer vibrations, λ = 0 . λ = 0 corresponding to zerofriction (i.e. RPMD). A crucial consequence of thischoice of friction matrix is that the centroid of the ringpolymer is unthermostatted, and the short-time errorof TRPMD from exact quantum dynamics is therefore O ( t ), the same as RPMD. C. Relation to Kramers Theory
To provide a qualitative description of the effect offriction on the TRPMD transmission coefficient, we ap-ply classical Kramers theory in the extended N F -dimensional ring polymer space, governed by dynamicson the (temperature-dependent) ring-polymer potentialenergy surface in Eq. (2). Since the short-time limit ofTRPMD rate theory is equal to QTST, and its long-timelimit invariant to the location of the dividing surface,TRPMD will give the QTST rate through the optimaldividing surface (defined as the surface which minimises k ‡ QM ( β )) , weighted by any recrossings of that dividingsurface by the respective dynamics. We express this us-ing the Bennett-Chandler factorization , k TRPMD ( β ) = k ‡∗ QM ( β ) lim t →∞ κ ∗ TRPMD ( t ) (18)where k ‡∗ QM ( β ) is the QTST rate, the asterisk denotesthat the optimal dividing surface f ∗ ( q ) is used and theTRPMD transmission coefficient is κ ∗ TRPMD ( t ) = (cid:82) d p (cid:82) d q e − β N H N ( p , q ) δ [ f ∗ ( q )] S ∗ N ( p , q ) h [ f ∗ ( q t )] (cid:82) d p (cid:82) d q e − β N H N ( p , q ) δ [ f ∗ ( q )] S ∗ N ( p , q ) h [ S ∗ N ( p , q )] (19)with analogous expressions to Eqs. (18) and (19) forRPMD. To examine the explicit effect of friction on theTRPMD rate we define the ratio χ λ ( β ) = k TRPMD ( β ) k RPMD ( β ) (20)and from Eq. (18) χ λ ( β ) = lim t →∞ κ ∗ TRPMD ( t ) κ ∗ RPMD ( t ) . (21)We then assume that the recrossing dynamics is dom-inated by one-dimensional motion through a parabolic saddle point on the ring-polymer potential energy sur-face, in which case the TRPMD transmission coefficientcan be approximated by the Kramers expression lim t →∞ κ ∗ TRPMD ( t ) (cid:39) (cid:113) α − α RP (22)where formally α RP = γ RP / ω RP , with γ RP the fric-tion along the reaction co-ordinate and ω RP the bar-rier frequency in ring-polymer space. For a general F -dimensional system finding f ∗ ( q ) and thereby comput-ing γ RP and ω RP is largely intractable. However, weexpect γ RP ∝ λ , and therefore define ˜ α RP = α RP /λ where the dimensionless parameter ˜ α RP is expected tobe independent of λ for a given system and tempera-ture, and represents the sensitivity of the TRPMD rateto friction. We further approximate that there is min-imal recrossing of the optimal dividing surface by the(unthermostatted) ring polymer trajectories such thatlim t →∞ κ ∗ RPMD ( t ) (cid:39) leading to χ λ ( β ) (cid:39) (cid:113) λ ˜ α − λ ˜ α RP . (23)Equation (23) relates the ratio of the TRPMD andRPMD rates as a function of λ with one parameter ˜ α RP ,and without requiring knowledge of the precise locationof the optimal dividing surface f ∗ ( q ). However, we canuse general observations concerning which ring-polymernormal modes contribute to f ∗ ( q ) to determine the likelysensitivity of the TRPMD rate to friction. Above thecrossover temperature into deep tunnelling, defined by β c = 2 π (cid:126) ω b (24)where ω b is the barrier frequency in the external potential V ( q ), the optimal dividing surface is well approximatedby the centroid f ∗ ( q ) = 1 N N − (cid:88) i =0 q i − q ‡ . (25)where q ‡ is the maximum in V ( q ). As the centroid is notthermostatted (since Ω = 0), in this regime γ RP = 0 =˜ α RP and we therefore predict from Eq. (23) that the ratewill be independent of λ , i.e. k TRPMD ( β ) (cid:39) k RPMD ( β ).Beneath the crossover temperature, the saddle point onthe ring-polymer potential energy surface bends into thespace of the first degenerate pair of normal modes. For symmetric systems, the optimal dividing surface isstill the centroid expression in Eq. (25) and (insofar as thereaction dynamics can be considered one-dimensional)˜ α RP (cid:39)
0, so k TRPMD ( β ) (cid:39) k RPMD ( β ).For asymmetric reactions, the optimal dividing surfaceis now a function of both the centroid and first degen-erate pair of normal modes (which are thermostatted) ,and we expect ˜ α RP >
0. From Eq. (23) the TRPMD ratewill decrease linearly with λ for small λ , for large frictionas λ − , and the ratio of the TRPMD to RPMD rates tobe a convex function of λ . This behaviour would alsobe expected for symmetric reactions beneath the secondcrossover temperature where the optimal dividing sur-face bends into the space of the second degenerate pairof normal modes. In all cases one would expect that in-creasing friction would either have no effect on the rate,or at sufficiently low temperatures cause it to decrease.It should be stressed that Eq. (23) is a considerablesimplification of the TRPMD dynamics and is not ex-pected to be reliable in systems where the ring poly-mer potential energy surface is highly anharmonic orskewed (such as F+H investigated below). In fact,even for a one-dimensional system, the minimum energy path on the N -dimensional ring polymer potential energysurface shows a significant skew beneath the crossovertemperature . The utility of Eq. (23) lies in its sim-plicity and qualitative description of friction-induced re-crossing. III. RESULTS
We initially study the benchmark one-dimensionalsymmetric and asymmetric Eckart barriers before pro-gressing to the multidimensional reactions H+H (sym-metric), D+MuH (quasisymmetrical), H+CH (asym-metric, polyatomic) and F+H (asymmetric, anhar-monic). A. One-dimensional results
The methodology for computation of TRPMD reac-tion rates is identical to that of RPMD , except for thethermostat attached to the internal normal modes of thering polymer, achieved using the algorithm in Ref. 79.The Bennett-Chandler factorization was employed, andthe same dynamics can be used for thermodymamic in-tegration along the reaction co-ordinate (to calculate theQTST rate) as to propagate trajectories (to calculate thetransmission coefficient). We firstly examine the symmetric Eckart barrier , V ( q ) = V sech ( q/a ) , (26)and to facilitate comparison with the literature ,use parameters to model the H+H reaction: V =0 . a = 0 . a , and m = 1061 m e , leading toa crossover temperature of k B β c = 2 . × − K − .The centroid reaction co-ordinate of Eq. (25) was usedthroughout. Results for a variety of temperatures andvalues of friction parameter λ are presented in Fig. 1,and values of ˜ α RP obtained by nonlinear least squares inTable I.Slightly beneath the crossover temperature ( k B β c =3 × − K − ), the TRPMD rate is indepedent of thevalue of friction (˜ α RP = 0), as predicted by Kramerstheory. Some sensitivity to λ is seen before twice thecrossover temperature, which is likely to be a breakdownof the one-dimensional assumption of Kramers theory;while the centroid is the optimal dividing surface, theminimum energy path bends into the space of the (ther-mostatted) lowest pair of normal modes . Beneath twicethe crossover temperature the friction parameter has asignificant effect on the rate, as to be expected from thesecond degenerate pair of normal modes becoming partof the optimal dividing surface . The functional formof χ λ ( β ) is also in accordance with the predictions ofKramers theory, monotonically decreasing as λ rises, andbeing a convex function of λ .Since RPMD underestimates the rate for this sym-metric reaction (and many others ), adding friction to TABLE I. Dimensionless friction sensitivity parameter ˜ α RP from Eq. (23), fitted by nonlinear least squares to simulation data.1D Eckart barriers Multidimensional reactionsSymmetric k B β/ − K − T /K 500 300 200 < β /a.u. 4 8 12 T /K 500 300 2000.00 0.06 0.17 H+CH –0.01 –0.01 0.00 k ( ) / m s - = 7 = 5 TRPMD( ) QM Kramers = 3 FIG. 1. Results for the symmetric Eckart barrier, showing theTRPMD result as a function of λ (red crosses), fitted Kramerscurve (green dashes) and quantum result (black line). β isquoted in units of k − − K − and the crossover temperatureis k B β c = 2 . × − K − . RPMD decreases its accuracy in approximating the quan-tum rate for this system.The asymmetric Eckart barrier is given by V ( q ) = A e − q/a + B cosh ( q/a ) (27)where A = − /π , B = 13 . /π and a = 8 / √ π in atomicunits ( (cid:126) = k B = m = 1), giving a crossover tempera-ture of β c = 2 π . To facilitate comparison with previousliterature the results are presented in Fig. 2 as = 12 = 8 c ( ) ) QM Kramers = 4 FIG. 2. Results for the asymmetric Eckart barrier quotedas c ( β ) [Eq. (28)], showing the TRPMD result as a functionof λ (red crosses), fitted Kramers curve (green dashes) andquantum result (black line). The crossover temperature is β c = 2 π a.u. the ratio c ( β ) = k ( β ) k clas ( β ) (28)and ˜ α RP values in Table I.Above the crossover temperature, TRPMD is invari-ant to the value of the friction parameter, and beneaththe crossover temperature, increasing λ results in a de-crease in the rate, such that TRPMD is closer to theexact quantum result than RPMD for all λ > λ is qual-itatively described by the crude Kramers approximation -3.0 -2.8 -2.6 -2.4 -2.2 -2.0020000400006000080000100000 c ( ) ‡ QTST RPMD TRPMD q FIG. 3. TRPMD (green dashes), RPMD (blue dots) andQTST (centroid dividing surface, red line) rates for the asym-metric Eckart barrier at β = 12, as a function of the dividingsurface q ‡ . (see Fig. 2), and it therefore seems that the improved ac-curacy of TRPMD could be a fortuitous cancellation be-tween the overestimation of the quantum rate by QTST,and the friction-induced recrossing of the optimal divid-ing surface by TRPMD trajectories. There is no particu-lar a priori reason to suppose that one value of λ shouldprovide superior results; from Fig. 2, at β = 8 a fric-tion parameter of λ = 1 .
25 causes TRPMD to equal thequantum result to within graphical accuracy, whereas at β = 12 this value of friction parameter causes overesti-mation of the rate, and further calculations (not shown)show that λ = 5 is needed for TRPMD and the quantumrates to agree.The numerical results also show a slightly higher cur-vature in k TRPMD ( β ) as a function of λ than Eq. (23)would predict, suggesting that the TRPMD rate reachesan asymptote at a finite value, rather than at zero asthe Kramers model would suggest. We suspect this isa breakdown of one-dimensional Kramers theory, sincein the λ → ∞ limit the system can still react via theunthermostatted centroid co-ordinate, but may have tosurmount a higher barrier on the ring polymer potentialenergy surface.We then investigate the effect of changing the loca-tion of the centroid dividing surface on the TRPMD rate.RPMD is already known to be invariant to the locationof the dividing surface , and we therefore choose a sys-tem for which TRPMD and RPMD are likely to differthe most, namely a low-temperature, asymmetric sys-tem where there is expected to be significant involvementof the thermostatted lowest degenerate pair of normalmodes in crossing the barrier. The asymmetric Eckartbarrier at β = 12 is therefore used as a particularlyharsh test, with the result plotted in Fig. 3. Although T = 200 KT = 300 K k ( T ) / c m s - TRPMD( ) QM Kramers T = 500 K
H+H FIG. 4. Results for the H+H reaction as a function of λ .Kramers is the fitted Kramers curve (see text). The crossovertemperature is 345K. the centroid-density QTST result varies by almost a fac-tor of six across the range of dividing surfaces considered( − ≤ q ‡ ≤ − . B. Multidimensional results
The results are calculated using adapted RPMDratecode , with details summarized in Table II. In the cal-culations reported below we used the potential energysurface developed by Boothroyd et al. (BKMP2 PES)for H+H and D+MuH, the Stark–Werner (SW) po-tential energy surface for F+H , and the PES-2008potential energy surface developed by Corchado et al.for H+CH . The computation of the free energy wasachieved using umbrella integration with TRPMDand checked against standard umbrella integration withan Andersen thermostat .H+H represents the simplest atom-diatom scatter-ing reaction and has been the subject of numerous TABLE II.
Input parameters for the TRPMD calculations on the H + H , D + MuH, and F + H Parameter Reaction ExplanationH + H D + MuH F + H H + CH Command line parameters
Temp
Nbeads
128 512 384 (200 K) 192 (250 K) Number of beads in the TRPMD calculations256 (300 K) 128 (300 K)64 (500 K) 64 (500 K)Dividing surface parameters R ∞
30 30 30 30 Dividing surface s parameter ( a ) N bonds N channel thermostat ’GLE/Andersen’ Thermostat for the QTST calculations λ
0; 0.25; 0.5; 0.75; 1.0; 1.5 Friction coefficient for the recrossing factor calculationsBiased sampling parameters N windows
111 111 111 111 Number of windows ξ -0.05 -0.05 -0.05 -0.05 Center of the first window dξ ξ N dt k i N trajectory
200 200 200 200 Number of trajectories t equilibration
20 20 20 20 Equilibration period (ps) t sampling
100 100 100 100 Sampling period in each trajectory (ps) N i × × × × Total number of sampling pointsPotential of mean force calculation ξ -0.02 -0.02 -0.02 -0.02 Start of umbrella integration ξ ‡ ∗ ∗ ∗ ∗ End of umbrella integration0.9904 (300 K) ∗ ∗ ∗ ∗ ∗ ∗ N bins dt t equilibration
20 20 20 20 Equilibration period (ps) in the constrained (parent)trajectory N totalchild t childsampling
20 20 20 20 Sampling increment along the parent trajectory (ps) N child
100 100 100 100 Number of child trajectories per oneinitially constrained configuration t child ∗ Detected automatically by RPMDrate. studies . The PES is symmetric and with a rel-atively large skew angle (60 ° ), and a crossover tempera-ture of 345K. The results in Fig. 4 show that the rate isessentially invariant to the value of λ above the crossovertemperature. At 300K there is a slight decrease in therate with increasing friction from 0 to 1.5 ( ∼
25 %), andthis is far more pronounced at 200K where the λ = 1 . λ = 0 (RPMD) result.D+MuH is “quasisymmetrical” since DMu and MuHhave very similar zero-point energies, and one wouldtherefore expect the RPMD rate to underestimate the ex-act quantum rate . Since it is Mu-transfer the crossovertemperature is very high (860 K) and therefore this reac- tion can be considered as a stress test for the deep tun-neling regime. The results in Fig. 5 show that friction inthe TRPMD dynamics causes further underestimation ofthe rate, especially at low temperatures; for λ = 1 . ∼
40% over the range of λ explored here.As an example of a typical asymmetric reaction, re-sults for H+CH are plotted in Fig. 6, which has acrossover temperature of 341K. RPMD is well-known tooverestimate the quantum rate for this system at lowtemperatures. Fig. 6 shows that above the crossovertemperature (500K) the friction parameter has a negligi-
T = 200 KT = 300 K k ( T ) / c m s - TRPMD( ) QM Kramers T = 500 K
D+MuH
FIG. 5. As for Fig. 4, but for the D+MuH reaction with ahigh crossover temperature of 860K. ble effect on the rate. As the temperature is decreasedbelow the crossover temperature (300K and 250K), thefriction induces more recrossings of the dividing surfaceand, as a result, the TRPMD rate approaches the exactquantum rate with increasing the friction parameter.Thus far, Kramers theory has been surprisingly suc-cessful at qualitatively explaining the behaviour of theTRPMD rate with increasing friction. Present resultswould suggest that TRPMD would therefore improveupon RPMD for all asymmetric reactions, where RPMDgenerally overestimates the rate beneath crossover .We then examine another prototypical asymmetric reac-tion, F+H , with a low crossover temperature of 264K.Fig. 7 shows that at 500K and 300K, the TRPMD rateis in good agreement with the quantum result, but in-creases very slightly with λ causing a spurious small neg-ative value of ˜ α RP in Table I. Beneath crossover, at 200Kthe rate is virtually independent of lambda, apart froma very slight increase around λ = 0 .
5. Consequently,TRPMD fares no better than RPMD for this system,contrary to the H+CH results and the predictions ofKramers theory. This is likely attributable to a highlyanharmonic and exothermic energy profile, and a veryflat saddle point in ring-polymer space .As can be seen from the graphs, the simple Kramers T = 250 KT = 300 K k ( T ) / c m s - TRPMD( ) QM Kramers T = 500 K
H+CH FIG. 6. Results for the H+CH reaction. The crossover tem-perature is 341K. prediction is in surprisingly good qualitative agreementwith the numerical results (apart from F+H beneathcrossover), even for the multidimensional cases, which isprobably attributable to those reactions being dominatedby a significant thermal barrier which appears parabolicon the ring-polymer potential energy surface, meaningthat the one-dimensional Kramers model is adequate forcapturing the friction-induced recrossing. In Table I the˜ α RP values, fitted to the numerical data, show that for agiven reaction ˜ α RP (cid:39) α RP increases asthe temperature is decreased. This can be qualitativelyexplained as the optimal dividing surface becoming moredependent on the thermostatted higher normal modesas the temperature is lowered . Not surprisingly, thehighest value of ˜ α RP is observed for the highly quantummechanical D+MuH reaction at 200K with ˜ α RP = 0 . T = 200 KT = 300 K k ( T ) / c m s - TRPMD( ) QM Kramers T = 500 K
F+H FIG. 7. Results for the anharmonic and asymmetric F+H reaction with a crossover temperature of 264K. IV. CONCLUSIONS
In this paper we have, for the first time, applied Ther-mostatted Ring Polymer Molecular Dynamics (TRPMD)to reaction rate theory. Regardless of the applied fric-tion, the long-time limit of the TRPMD flux-side time-correlation function (and therefore the TRPMD rate) isindependent of the location of the dividing surface, andits short-time limit is equal to rigorous QTST .In section II C we use Kramers theory to predictthat, above the crossover temperature, the RPMD andTRPMD rates will be similar, and beneath crossover theTRPMD rate for asymmetric systems will decrease with λ , and the same effect should be observed for symmetricsystems beneath half the crossover temperature.TRPMD rate theory has then been applied to thestandard one-dimensional model systems of the symmet-ric and asymmetric Eckart barriers, followed by the bi-molecular reactions H+H , D+MuH, H+CH and F+H .For all reactions considered, above the crossover tem-perature the TRPMD rate is virtually invariant to thevalue of λ and therefore almost equal to RPMD, as pre-dicted by Kramers theory. Beneath the crossover temper-ature, most asymmetric reactions show a decrease in the TRPMD rate as λ is increased, and in qualitative agree-ment with the Kramers prediction in Eq. (23). A simi-lar trend is observed for symmetric reactions, which alsoshow some diminution in the rate with increasing frictionabove half the crossover temperature ( β c < β < β c ),probably due to the skewed ring-polymer PES caus-ing a breakdown in the one-dimensional assumption ofKramers theory. For the asymmetric and anharmoniccase of F+H , beneath the crossover temperature thereis no significant decrease in the rate with increased fric-tion, illustrating the limitations of Kramers theory.These results mean that beneath the crossover tem-perature TRPMD will be a worse approximation to thequantum result than RPMD for symmetric and qua-sisymmetrical systems (where RPMD underestimates therate ), and TRPMD will be closer to the quantumrate for asymmetric potentials (where RPMD overesti-mates the rate). However, the apparent increase in accu-racy for asymmetric systems appears to be a cancellationof errors from the overestimation of the quantum rate byRPMD which is then decreased by the friction in thenon-centroid normal modes of TRPMD, and there is noa priori reason to suppose that one effect should equalthe other for any given value of λ .Although the above results do not advocate the useof TRPMD rate theory as generally being more accu-rate than RPMD, TRPMD rate calculation above thecrossover temperature may be computationally advanta-geous in complex systems due to more efficient samplingof the ring-polymer phase space by TRPMD trajecto-ries than RPMD trajectories . TRPMD may there-fore provide the same accuracy as RPMD rate calcula-tion at a lower computational cost, and testing this inhigh-dimensional systems where RPMD has been suc-cessful, such as complex-forming reactions , sur-face dynamics , and enzyme catalysis would be a use-ful avenue of future research.Future work could also include non-adiabaticsystems , applying a thermostat to thecentroid to model a bath system , and generalizationsto non-Markovian friction using Grote-Hynes theory .In closing, present results suggest that TRPMD canbe used above the crossover temperature for thermallyactivated reactions, and beneath crossover further testingis required to assess its utility for asymmetric systems. V. ACKNOWLEDGEMENTS
TJHH acknowledges a Research Fellowship from Je-sus College, Cambridge, and helpful comments onthe manuscript from Stuart Althorpe. YVS acknowl-edges support via the Newton International AlumniScheme from the Royal Society. YVS also thanksthe European Regional Development Fund and theRepublic of Cyprus for support through the Re-search Promotion Foundation (Project Cy-Tera NEAΓΠO∆OMH/ΣTPATH/0308/31).1
Appendix A: Detailed Balance
For a homogeneous Markov process such as TRPMDfor which negative time is not defined , detailed balanceis defined as P ( p (cid:48) , q (cid:48) , t | p , q , ρ s ( p , q )= P ( − p , q , t | − p (cid:48) , q (cid:48) , ρ s ( p (cid:48) , q (cid:48) ) (A1)where ρ s ( p , q ) = e − β N H N ( p , q ) is the stationary distribu-tion and P ( p (cid:48) , q (cid:48) , t | p , q ,
0) is the conditional probabilitythat a ring polymer will be found at point ( p (cid:48) , q (cid:48) ) at time t , given that is was at ( p , q ) at time t = 0.To demonstrate that Eq. (A1) is statisfied, we rewritethe Fokker-Planck operator Eq. (5) as A N = − N − (cid:88) j =0 (cid:18) ∂∂q j a ( p , q ) j + ∂∂p j b ( p , q ) j (cid:19) + 12 N − (cid:88) j =0 N − (cid:88) j (cid:48) =0 ∂∂p j ∂∂p j (cid:48) C ( p , q ) jj (cid:48) (A2)where the vectors a ( p , q ) = p /m , b ( p , q ) = − U N ( q ) ←−∇ q − Γ · p and the matrix C ( p , q ) = 2 m Γ /β N .Note that the derivatives in Eq. (A2) act on a ( p , q ), b ( p , q ) or C ( p , q ) and whatever follows them which isacted upon by A N .The necessary and sufficient conditions for detailedbalance [Eq. (A1)] to hold, in addition to ρ s ( p , q ) be-ing a stationary distribution, are then given by a ( − p , q ) ρ s ( p , q ) = − a ( p , q ) ρ s ( p , q ) (A3) − b ( − p , q ) T ρ s ( p , q ) = − b ( p , q ) T ρ s ( p , q )+ ∇ p · C ( p , q ) ρ s ( p , q ) (A4) C ( − p , q ) = C ( p , q ) (A5)Condition Eq. (A3) is trivially satisfied. Provided thatthe friction matrix is even w.r.t. momenta (satisfied hereas Γ is not a function of p ) Eq. (A5) will be satisfied.Eq. (A4) becomes( Γ · p ) T ρ s ( p , q ) = − mβ N ∇ p · Γ ρ s ( p , q ) (A6)which is satisfied with ρ s ( p , q ) = e − β N H N ( p , q ) and thefriction matrix used here.Given that Eq. (A1) is satisfied, for an arbitrary cor-relation function one can then show C TRPMD AB ( t ) = 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q (cid:90) d p (cid:48) (cid:90) d q (cid:48) e − β N H N ( p , q ) A ( p , q ) P ( p (cid:48) , q (cid:48) , t | p , q , B ( p (cid:48) , q (cid:48) ) (A7)= 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q (cid:90) d p (cid:48) (cid:90) d q (cid:48) e − β N H N ( p , q ) A ( − p (cid:48) , q (cid:48) ) P ( p (cid:48) , q (cid:48) , t | p , q , B ( − p , q ) (A8)and for the Langevin trajectories considered here, whichare continuous but not differentiable, this means C TRPMD AB ( t ) = 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − β N H N ( p , q ) × A ( p , q ) ¯ B ( p t , q t ) (A9)= 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − β N H N ( p , q ) × ¯ A ( − p t , q t ) B ( − p , q ) (A10)where q t ≡ q t ( p , q , t ) is the vector of positions stochas-tically time-evolved according to Eqs. (3) and (4), and¯ B ( p t , q t ) = (cid:90) d p (cid:48) (cid:90) d q (cid:48) P ( p (cid:48) , q (cid:48) , t | p , q , B ( p (cid:48) , q (cid:48) )(A11)with ¯ A ( p t , q t ) similarly defined. Appendix B: Independence of k TRPMD ( β ) to the dividingsurface location We use a similar methodology to that which Craig andManolopoulos employed for RPMD , and give the mainsteps here. We firstly differentiate the side-side correla-tion function in Eq. (9) w.r.t. the location of the dividingsurface q ‡ (or any other parameter specifying the natureof the dividing surface), giving ddq ‡ C ss ( t ) = 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − β N H N ( p , q ) × ∂f ( q ) ∂q ‡ { δ [ f ( q )] h [ f ( q t )] + h [ f ( q )] δ [ f ( q t )] } . (B1)Since TRPMD dynamics obeys detailed balance (asshown in appendix A), and the dividing surface is onlya function of position, the second term on the RHS of2Eq. (B1) is identical to the first, ddq ‡ C ss ( t ) = 2(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − β N H N ( p , q ) × ∂f ( q ) ∂q ‡ δ [ f ( q )] h [ f ( q t )] . (B2)Differentiation of Eq. (B2) w.r.t. time using Eq. (6),and relating the side-side and flux-side functions usingEq. (10), yields ddq ‡ C fs ( t ) = − π (cid:126) ) N (cid:90) d p (cid:90) d q e − β N H N ( p , q ) × ∂f ( q ) ∂q ‡ δ [ f ( q )] δ [ f ( q t )] S N ( p t , q t ) . (B3)Equation (B3) corresponds to a trajectory commencingat the dividing surface at time zero and returning to itat time t with non-zero velocity S N ( p t , q t ). At finitetimes while there is recrossing of the barrier, there willbe trajectories satisfying these conditions, but after theplateau time when no trajectories recross the barrier [cf.Eq. (14)], these conditions are clearly not satisfied, andthe rate will be independent of the location of the divid-ing surface. This proof is valid for any friction matrix which satis-fies the detailed balance conditions of appendix A, anddoes not require the presence of ring-polymer springs inthe potential, so is valid for any classical-like reactionrate calculation using Langevin dynamics. E. Pollak and P. Talkner,
Chaos (2005), 026116. P. H¨anggi, P. Talkner and M. Borkovec,
Rev. Mod. Phys. (1990), 251. T. Yamamoto,
J. Chem. Phys. (1960), 281. W. H. Miller,
J. Chem. Phys. (1974), 1823. W. H. Miller, S. D. Schwartz and J. W. Tromp,
J. Chem. Phys. (1983), 4889. W. H. Miller,
Acc. Chem. Res. (1993), 174. W. H. Miller,
J. Chem. Phys. (1975), 1899. S. Chapman, B. C. Garrett and W. H. Miller,
J. Chem. Phys. (1975), 2710. I. Affleck,
Phys. Rev. Lett. (1981), 388. G. Mills, G. Schenter, D. Makarov and H. J´onsson,
Chem. Phys.Lett. (1997), 91 . G. Mil’nikov and H. Nakamura,
Phys. Chem. Chem. Phys. (2008), 1374. S. Andersson, G. Nyman, A. Arnaldsson, U. Manthe andH. J´onsson,
J. Phys. Chem. A (2009), 4468. J. O. Richardson and S. C. Althorpe,
J. Chem. Phys. (2009), 214106. S. C. Althorpe,
J. Chem. Phys. (2011), 114104. T. Kawatsu and S. Miura,
J. Chem. Phys. (2014), 024101. J. B. Rommel, T. P. M. Goumans and J. K¨astner,
J. Chem.Theor. Comput. (2011), 690. Y. Zhang, J. B. Rommel, M. T. Cvitas and S. C. Althorpe,
Phys. Chem. Chem. Phys. (2014), 24292. J. Van´ıˇcek, W. H. Miller, J. F. Castillo and F. J. Aoiz,
J. Chem.Phys. (2005), 054108. W. H. Miller, Y. Zhao, M. Ceotto and S. Yang,
J. Chem. Phys. (2003), 1329. H. Eyring,
J. Chem. Phys. (1935), 107. H. Eyring,
Chem. Rev. (1935), 65. E. Wigner,
Trans. Faraday Soc. (1938), 29. D. G. Truhlar, B. C. Garrett and S. J. Klippenstein,
J. Phys.Chem. (1996), 12771. D. G. Truhlar and B. C. Garrett,
Acc. Chem. Res. (1980),440. G. A. Voth, D. Chandler and W. H. Miller,
J. Chem. Phys. (1989), 7749. J. Liu and W. H. Miller,
J. Chem. Phys. (2009), 074113. J. Liu,
Int. J. Quantum Chem. (2015), published online, doi:10.1002/qua.24872. Q. Shi and E. Geva,
J. Chem. Phys. (2003), 8173. J. Cao and G. A. Voth,
J. Chem. Phys. (1993), 10070.J. Cao and G. A. Voth, J. Chem. Phys. (1994), 5093.J. Cao and G. A. Voth,
J. Chem. Phys. (1994), 5106.J. Cao and G. A. Voth,
J. Chem. Phys. (1994), 6157.J. Cao and G. A. Voth,
J. Chem. Phys. (1994), 6168.J. Cao and G. A. Voth,
J. Chem. Phys. (1994), 6184. G. A. Voth,
J. Phys. Chem. (1993), 8365. G. A. Voth,
Path-Integral Centroid Methods in Quantum Statis-tical Mechanics and Dynamics , Adv. Chem. Phys., John Wiley& Sons, Inc. (1996). Q. Shi and E. Geva,
J. Chem. Phys. (2002), 3223. T. D. Hone, P. J. Rossky and G. A. Voth,
J. Chem. Phys. (2006), 154103. S. Jang and G. A. Voth,
J. Chem. Phys. (1999), 2371. S. Jang and G. A. Voth,
J. Chem. Phys. (1999), 2357. I. R. Craig and D. E. Manolopoulos,
J. Chem. Phys. (2004),3368. I. R. Craig and D. E. Manolopoulos,
J. Chem. Phys. (2005),084106. I. R. Craig and D. E. Manolopoulos,
J. Chem. Phys. (2005),034102. S. Habershon, D. E. Manolopoulos, T. E. Markland and T. F.Miller,
Annu. Rev. Phys. Chem. (2013), 387. N. Boekelheide, R. Salom´on-Ferrer and T. F. Miller,
PNAS (2011), 16159. J. S. Kretchmer and T. F. Miller III,
J. Chem. Phys. (2013),134109. A. R. Menzeleev, N. Ananth and T. F. Miller III,
J. Chem.Phys. (2011), 074106. A. R. Menzeleev, F. Bell and T. F. Miller,
J. Chem. Phys. (2014), 064103. R. Collepardo-Guevara, Y. V. Suleimanov and D. E.Manolopoulos,
J. Chem. Phys. (2009), 174713. R. Collepardo-Guevara, Y. V. Suleimanov and D. E.Manolopoulos,
J. Chem. Phys. (2010), 049902. Y. V. Suleimanov, R. Collepardo-Guevara and D. E.Manolopoulos,
J. Chem. Phys. (2011), 044131. Y. V. Suleimanov,
J. Phys. Chem. C (2012), 11141. Y. V. Suleimanov, W. J. Kong, H. Guo and W. H. Green,
J.Chem. Phys. (2014), 244103. R. P´erez de Tudela, Y. V. Suleimanov, J. O. Richardson,V. Saez-Rabanos, W. H. Green and F. J. Aoiz,
J. Phys. Chem.Lett. (2014), 4219. Y. Li, Y. V. Suleimanov, M. Yang, W. H. Green and H. Guo,
J. Phys. Chem. Lett. (2013), 48. Y. Li, Y. V. Suleimanov, J. Li, W. H. Green and H. Guo,
J.Chem. Phys. (2013), 094307. Y. V. Suleimanov, R. P. de Tudela, P. G. Jambrina, J. F.Castillo, V. Saez-Rabanos, D. E. Manolopoulos and F. J. Aoiz,
Phys. Chem. Chem. Phys. (2013), 3655. J. W. Allen, W. H. Green, Y. Li, H. Guo and Y. V. Suleimanov,
J. Chem. Phys. (2013), 221103. R. P´erez de Tudela, F. J. Aoiz, Y. V. Suleimanov and D. E.Manolopoulos,
J. Phys. Chem. Lett. (2012), 493. R. P. de Tudela, Y. V. Suleimanov, M. Menendez, J. F. Castilloand F. J. Aoiz,
Phys. Chem. Chem. Phys. (2014), 2920. Y. Li, Y. V. Suleimanov, W. H. Green and H. Guo,
J. Phys.Chem. A (2014), 1989. J. Espinosa-Garcia, A. Fernandez-Ramos, Y. V. Suleimanov and J. C. Corchado,
J. Phys. Chem. A (2014), 554. E. Gonzalez-Lavado, J. C. Corchado, Y. V. Suleimanov, W. H.Green and J. Espinosa-Garcia,
J. Phys. Chem. A (2014),3243. Y. Li, Y. V. Suleimanov and H. Guo,
J. Phys. Chem. Lett. (2014), 700. Y. V. Suleimanov and J. Espinosa-Garcia,
J. Phys. Chem. B (2015), published online, doi: 10.1021/acs.jpcb.5b02103. K. M. Hickson, J.-C. Loisin, H. Guo and Y. V. Suleimanov,
Phys. Rev. Lett. (2015), (submitted). Y. Suleimanov, J. Allen and W. Green,
Comp. Phys. Comm. (2013), 833 . T. J. H. Hele and S. C. Althorpe,
J. Chem. Phys. (2013),084108. S. C. Althorpe and T. J. H. Hele,
J. Chem. Phys. (2013),084115. T. J. H. Hele and S. C. Althorpe,
J. Chem. Phys. (2013),084116. M. J. Gillan,
J. Phys. C (1987), 3621. M. J. Gillan,
Phys. Rev. Lett. (1987), 563. E. Geva, Q. Shi and G. A. Voth,
J. Chem. Phys. (2001),9209. J. Richardson,
Ring-Polymer Approaches to Instanton Theory ,Ph.D. thesis, Cambridge University (2012). Y. Zhang, T. Stecher, M. T. Cvitaˇs and S. C. Althorpe,
J. Phys.Chem. Lett. (2014), 3976. R. Kubo,
J. Phys. Soc. Jpn. (1957), 570. T. J. H. Hele, M. J. Willatt, A. Muolo and S. C. Althorpe,
J.Chem. Phys. (2015), 134103. T. J. H. Hele, M. J. Willatt, A. Muolo and S. C. Althorpe,
J.Chem. Phys. (2015), 191101. T. Matsubara,
Progress of Theoretical Physics (1955), 351. A. Witt, S. D. Ivanov, M. Shiga, H. Forbert and D. Marx,
J.Chem. Phys. (2009), 194510. S. D. Ivanov, A. Witt, M. Shiga and D. Marx,
J. Chem. Phys. (2010), 031101. G. Bussi and M. Parrinello,
Phys. Rev. E (2007), 056707. M. Rossi, M. Ceriotti and D. E. Manolopoulos,
J. Chem. Phys. (2014), 234116. M. Ceriotti, M. Parrinello, T. E. Markland and D. E.Manolopoulos,
J. Chem. Phys. (2010), 124104. H. Kramers,
Physica (1940), 284 . R. P. Feynman and A. R. Hibbs,
Quantum Mechanics and PathIntegrals , McGraw-Hill (1965). R. Zwanzig,
Nonequilibrium statistical mechanics , Oxford Uni-versity Press, New York (2001). D. Chandler,
J. Chem. Phys. (1978), 2959. We note that the flux-side function in the gas phase was orig-inally derived as an asymmetric-split trace and later as asymmetric-split trace , both of which are equivalent to the Kubotransformed trace given here in the t → ∞ limit. See e.g. Fig. 5 of Ref. 47. T. J. H. Hele,
Quantum Transition-State Theory , Ph.D. thesis,University of Cambridge (2014). M. Rossi, H. Liu, F. Paesani, J. Bowman and M. Ceriotti,
J.Chem. Phys. (2014), 181101. B. J. Braams and D. E. Manolopoulos,
J. Chem. Phys. (2006), 124105. D. Frenkel and B. Smit,
Understanding Molecular Simulation ,Academic Press (2002). A. Nitzan,
Chemical Dynamics in Condensed Phases , OxfordUniversity Press, New York (2006). E. Pollak,
J. Chem. Phys. (1986), 865. In practice, both the location of the dividing surface and thedegree of recrossing through it are very difficult to determinenumerically, though Richardson and Althorpe found that forthe symmetric and asymmetric Eckart barrier at almost twicethe inverse crossover temperature ( β (cid:126) ω b = 12) fewer than 20%of trajectories recrossed. C. Eckart,
Phys. Rev. (1930), 1303. E. Pollak and J.-L. Liao,
J. Chem. Phys. (1998), 2733. S. Jang and G. A. Voth,
J. Chem. Phys. (2000), 8747. This proof was originally for a gas-phase system and a cen-troid dividing surface, but this can easily be extended to thecondensed phase provided there is the necessary separation oftimescales . A. I. Boothroyd, W. J. Keogh, P. G. Martin and M. R. Peterson,
J. Chem. Phys. (1996), 7139. K. Stark and H.-J. Werner,
J. Chem. Phys. (1996), 6515. J. C. Corchado, J. L. Bravo and J. Espinosa-Garcia,
J. Chem.Phys. (2009), 184314.
J. K¨astner and W. Thiel,
J. Chem. Phys. (2005), 144104.
J. K¨astner and W. Thiel,
J. Chem. Phys. (2006), 234106.
H. C. Andersen,
J. Comput. Phys. (1983), 24 . J. Jankunas, M. Sneha, R. N. Zare, F. Bouakline, S. C. Al-thorpe, D. Herr´aez-Aguilar and F. J. Aoiz,
PNAS (2014),15.
T. Stecher,
Benchmark Studies of Ring Polymer Molecular Dy-namics Rate Theory , Ph.D. thesis, Downing College, Universityof Cambridge (2010).
M. Rossi and D. E. Manolopoulos, private communication,(2015).
T. J. H. Hele,
An electronically non-adiabatic generalization ofring polymer molecular dynamics , Master’s thesis, University ofOxford (2011).
N. Ananth,
J. Chem. Phys. (2013), 124102.
N. Ananth and T. F. Miller,
J. Chem. Phys. (2010), 234103.
J. O. Richardson and M. Thoss,
J. Chem. Phys. (2013),031102.
J. O. Richardson and M. Thoss,
J. Chem. Phys. (2014),074106.
R. F. Grote and J. T. Hynes,
J. Chem. Phys. (1980), 2715. C. Gardiner,