New Stable and Fast Ring-Polymer Molecular Dynamics for Calculating Bimolecular Rate Coefficients with Example of OH + CH_4

The accurate and efficient calculation of the rate coefficients of chemical reactions is a key issue in the research of chemical dynamics. In this work, by applying the dimension-free ultra-stable Cayley propagator, the thermal rate coefficients of a prototypic high dimensional chemical reaction OH + CH4 ->H2O + CH3 in the temperature range of 200 K to 1500 K are investigated with ring polymer molecular dynamics (RPMD), on a highly accurate full-dimensional potential energy surface. Kinetic isotope effects (KIEs) for three isotopologues of the title reaction are also studied. The results demonstrate excellent agreement with experimental data, even in the deep tunneling region. Especially, the Cayley propagator shows high calculation efficiency with little loss of accuracy. The present results confirmed the applicability of the RPMD method, particularly the speed-up by Cayley propagator, in theoretical calculations of bimolecular reaction rates.


Introduction
The rate coefficients of radical-molecule reactions attract wide interest in both atmospheric and combustion chemistry, because of their importance in modeling the kinetics and revealing the mechanism. In particular, the kinetic isotope effects (KIEs) calculated for the same reaction with different isotopes shed light on the reaction mechanism and quantum effects. However, the precise measurement of rate coefficients and KIEs is often a difficult task and the data are not always available. [1] So, theoretical calculations are necessary, since they not only offer reliable predictions of rate coefficients and KIEs but also provide microscopic insight into the reaction mechanism.
Currently, there are several categories of methodologies in calculating rate coefficients. The first approach is accurate quantum dynamics (QD) with wave packets, but it suffers from high computational costs. The second method is quasi-classical trajectory (QCT), which is fast, but incapable of treating quantum effects such as tunneling. The third method is the transition state theory (TST), which requires information only near the barrier. However, TST treats the tunneling approximately and may miss important recrossing dynamics near the barrier. The last one is ringpolymer molecular dynamics (RPMD), which achieves a good balance between accuracy and efficiency. In the past decade, RPMD has evidenced a myriad of successful applications for calculating bimolecular reaction rate coefficients. [2] The reaction of OH + CH4 → H2O + CH3 is a prototypical polyatomic reaction with 15 degrees of freedom (DOF). It is a major channel for methane removal in the atmosphere at low temperatures and is also a key step in combustion chemistry at high temperatures. Extensive experiments have been performed on title reaction [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17], yielding thermal rate coefficients, KIEs, and even state-resolved rates over a large range of temperature. However, full dimensional QD calculations have not been reported yet due to the high dimensionality. On the other hand, due to the multidimensional tunneling of the transferring H atom, TST and QCT calculations could contain substantial inaccuracies at low temperatures.
Recently, a high accuracy full-dimensional neural-network potential energy surface (PIP-NN PES) [18] has been reported. With variational transition state theory (VTST) and quasi-classical trajectory (QCT), Li and Guo calculated the rate coefficients and KIEs and compare them with experimental data at temperatures ranging from 200 K to 1000 K [19]. Relatively large errors were discovered when the temperature drops below 300 K. This discrepancy presumably stems from the lack of multi-dimensional tunneling in both VTST and QCT. More specifically, the anharmonicity leads to a lower rate, and the lack of recrossing leads to a higher result, a combination of both caused error cancellation at high temperatures. This cancellation makes the calculated rate coefficients close to the experiments but the KIEs contain large errors. When temperature falls into the deep tunneling region, tunneling, recrossing, anharmonicity and zeropoint energy (ZPE) effects become substantial, resulting in larger deviations from experiments.
Considering the above problems, a proper method based on a correct description of rate coefficients should reproduce both quantum and recrossing effects. The RPMD based on the isomorphism between a set of ring-polymers in extended phase space and the original quantum system, using real-time dynamics to calculate time correlation functions for obtaining the rate coefficients [20], is a perfect candidate.
Previously, RPMD has been used for the investigation of title reaction on a less accurate analytical PES, [21,22] where the rate coefficients at low temperatures were severely overestimated by one order of magnitude. Given the reliability of RPMD in other bimolecular reactions, such deficiency is presumably attributable to the deficiency of the PES. In this work, we report RPMD calculations using the accurate PIP-NN PES. However, the PIP-NN PES is much more computationally expensive due to its complex functional form, hence numerical acceleration of the RPMD propagation becomes critical. [23,24] Recently, the Cayley propagator [23,24] was proposed for RPMD calculation to endow strong stability, ergodicity, and dimensionality freedom.
Such properties imply the ability to speed up the calculation by employing a larger time step without loss of accuracy. In this work, we deploy Cayley propagator in our RPMD calculations.
We organized this work as follows. The method is summarized in Sec. 2. All results are presented and discussed in Sec. 3, and the conclusion is drawn in Sec. 4.

RPMD
The theoretical method is summarized thoroughly elsewhere [2,25], so we only give a brief description. For the title reaction, the quantum Hamiltonian is: (1) where ˆi p and ˆi q are the momentum and position operators of the i th atom with its mass i m . The RPMD Hamiltonian can be written down based on the isomorphism between the quantum system and the corresponding classical ring-polymer system: here each atom is presented as a ring formed by a set of P beads. The cyclic boundary is used: 0 P = qq , and the force constant of harmonic potential connecting adjacent beads is , with the reciprocal temperature of the system b P = b / P and 6 Making use of the Bennett-Chandler factorization [26,27], the RPMD rate coefficient is given The first term ‡ QTST ( ; ) kT  is the quantum transition state theory (QTST) rate coefficients, which is evaluated at the peak of the free energy barrier in the reaction coordinate ‡  .
The reaction coordinate (RC) of the title reaction is constructed from a pair of dividing surfaces defined in terms of the centroid of the ring polymer. The first dividing surface is defined by a distance R  , which is the distance between the centers of mass (COM) of both reactants OH and CH4, and is large enough making the interaction between them negligible: [29,30] s 0 q where r AB is the centroid position vector of group AB's COM And the centroid position vector is obtained from averaging the position vectors of beads The second dividing surface is located at the transition state region, [29,30] (6) In practice, ‡

QTST ( ; )
kT  is determined from the centroid potential of the mean force (PMF), The term ‡ in the exponential is the free-energy difference between both dividing surfaces, usually obtained by using umbrella integration (UI) [31,32] along  . And  is the reduced mass of reactants. The transmission coefficient, ‡ ) (; t  → , provides dynamical correction evaluated by the ratio between long-time limit and zero-time limit of the flux-side correlation function: which accounts for recrossing at the bottleneck ( x ‡ ) with ZPE.
The number of beads determines the accuracy of the simulation. More beads would lead more accuracy, but less computational efficiency. And there is also a formula to predict the minimum number of beads needed for a given system: (9) where max  is the largest vibrational frequency of the system. In this work, it is the frequency of OH bond stretching, 1 3757 cm − . [18] In each system at each temperature, P min is calculated first and is rounded up to be power of 2. It worth to point out, at low temperatures, the number of beads needed would become large, making calculation time-consuming.
Another key parameter, the cross-over temperature, also needs considering with: where b i is the imaginary frequency of the reaction system at the bottleneck. The system below c T is considered at the so-called deep-tunneling region, where RPMD would give results with larger error, and much more number of beads than min P needed to obtain results accurate enough.
The crossover temperature for the title reaction is T c = 326 K with . [18] The final rate coefficients should be corrected by a ratio of the electronic partition function TS elec. reactants elec.

Cayley Propagator
The traditional propagation of MD is Verlet and its various. These propagation schemes were proven to be not stable and ergodic in RPMD calculations.
In the original path integral molecular dynamics (PIMD), the integration scheme propagating (13) where and are the coordinates and momentum for P beads respectively, and i L is the classical time evolution operator of the i th atom. Using the Trotter expansion, the evolution induced by the free ring-polymer part, and the external potential part can be separated with an 2 ) ( t  residual error, where A , B and O are the time evolution operator of pure free ring-polymer, the time evolution operator of the system on potential energy surface, and thermostat. This expansion scheme is denoted as "BAOAB", and is recommended in a recent study. [34] A recent work [23] revealed that such a "BAOAB" integrator is not strongly stable with respect to the external potential, which causes poor ergodicity when using a large timestep and a large number of beads.
The "BCOCB" integrator is then introduced in the following study [24] to solve this issue by replacing the operator "A" from the exponential propagator to the Cayley propagator: Such "BCOCB" integrator is shown to be not only strongly stable with respect to the external potential but also provides the exact equilibrium position marginal distribution in the phase space for harmonic external potentials. It is later shown that the "BCOCB" integrator is the unique one in the "BAOAB" family that can achieve such good properties. [35] As mentioned above in equation (3), (7) and (8)

Computational Details
The globally full-dimensional PIP-NN PES with high accuracy developed by Li et al. [18] (PIP-NN PES) was used in our RPMD calculations. It's constructed based on about 135,000 points at the level of UCCSD(T)-F12a/AVTZ, whose root-mean-square fitting error is less than 0.09 kcal/mol.
All simulations about rate coefficients reported here used RPMDrate [29] and with our customizations.
The number of beads used in simulations is based on Eq.(9), and rounded up to the power of 2.
During the sampling, the force constants needed for windows near the transition state . So, we reduced the size of windows in that region ( 0.005

 =
). In each window, the system was equilibrated for 0.2 ps, following a production run of 10 ps per trajectory. Here 30 trajectories were employed in each window and equivalent to 300 ps for the system totally. Such choice is a balance between accuracy and efficiency, and the convergence validation is shown in SI. The Andersen thermostat was used in all NVT ensemble simulations. The time step used for traditional ring-polymer propagator at 500K T  is 0.1 fs, 200 K to 300 K is 0.2 fs, and for Cayley propagator is 0.5 fs.
In the calculation of transmission coefficients, a long parent trajectory constrained at the peak of PMF was simulated 10 ns using the SHAKE algorithm, with the purpose of generating initial configurations that were sampled every 0.2 ps. For each of these constrained configurations, 100 child trajectories were then spawned at the random initial momentum from the Boltzmann distribution. Each child trajectory was then propagated without constraint for 0.5 ps to reach the "plateau" time.
At each temperature, the RPMD calculations were firstly performed with one bead, which provides the classical limit, and then with a sufficient number of beads ensuring convergence.
along the MEP were also obtained using POLYRATE. It is the vibrational zero-point energy of the generalized normal-mode analysis.

Results and Discussion
The minimum energy path (MEP) and the ground state adiabatic energy G a V for all four isotopic reactions are shown in Figure 1. For MEP, all the 4 reactions are exothermic and with highly similar shape, since in this plot, the mass-scaled reaction coordinate [39]

OH + CH4
The rate coefficients for the OH + CH4 reaction and its isotopes, calculated by RPMD and other theoretical methods, and measured from experiments, are summarized in Table 1. The rate coefficients for the reaction OH+CH4 are also depicted in Figure 2, with results from other theoretical methods and experimental results. In detail, all collected experimental results are plotted in Figure 2(a), and the theoretical ones are shown in Figure 2(b). First of all, at low temperature, the reactants are trapped in the vdW well at first, forming a complex with orientation favoring to the reaction, leaving the rate coefficients larger than those predicted by Arrhenius' law.
Another feature of the title reaction is that the recrossing is significant in all the four isotopic reactions, which can be observed in the table. This phenomenon would be attributed to the small skew angle of the title reaction. Here the word recrossing is used loosely since it contains not only the classical recrossing over the barrier, but also the quantum tunneling effect which is reflected by the smaller transmission coefficient value at lower temperature. This is also the characteristic feature of the heavy-light-heavy reactions, as discovered in O/Cl+CH4 reactions [41,42].
Furthermore, the small skewing angle would also lead to a quantum reactive scattering resonance [42], which would be another source of deviation between RPMD results and experimental one. Unfortunately this system is too large for time-dependent quantum dynamics calculations, it cannot elucidate the effect of quantum resonance in the title reactions till now.
Another interesting feature of the RPMD results is the position of the peak of PMF ( x ‡ ), which is increasing with temperature, in contrast to the reactions of Cl/O/H+CH4 [41][42][43]. This can be understood as a significant van der Waals (vdW) well on PES before the TS in the reactant region.
This trend is also evidenced in the reactions of HCl + OH [44] and S + H2 [45].
From both the Table 1 and Figure 1, the RPMD rate results in this work agree perfectly with experimental results among the entire range of temperatures investigated, 200 K to 1500 K.
Considering when the temperature is below c T (326 K for title reaction), the RPMD would overestimate the rate coefficients due to the asymmetric barrier [46], our results coincide with this prediction. Even the error becomes larger at 200 K, the relative error to the experimental values are still less than 20%, proving again RPMD is a highly reliable method for obtaining rate coefficients.
The RPMD results based on other PESs have great errors when the temperature goes below c T . To better understand the above results, it's helpful to recall the details of other PESs used in previous works. The PES-2014 is fitted thoroughly from ab initio points, but the PES-2000 is constructed with both ab initio points and experimental data. [47] The height of the energy barrier (6.8 kcal/mol, from PES-2014) from the left van der Waals (vdW) well is lower than the value (7.53 kcal/mol) calculated by UCCSD(T)-F12a/AVTZ [18], which would make the rate coefficients higher than the experimental value at low temperatures. On the other hand, the height of TS (6.4 kcal/mol, also from PES-2014) is higher than the ab initio value (6.3 kcal/mol), which leads to lower rate coefficients at high temperatures, since higher temperature corresponds to higher kinetic energy, which will smear the vdW well and the height of barrier plays the determinate role.   Table   1, mentioned above.
The transmission coefficients as a function of time are shown in the right panel of Figure 3. From the plots, the "plateau" values of k(t) are relatively small as mentioned above. Especially at low temperatures, the plateau values are 0.383 for 200 K and 0.482 for 300 K, similar as observed in a previous study for H + CH4. [43] All the transmission coefficients converged to their plateau values after a short period of fluctuation, which is from the choice of the dividing surfaces, and with less physical meaning [48].
reducing with increasing of temperature, means the results are closer to the classical limit. Figure   4 also shows comparison of the  from RPMD and from classical MD in the right panel.  from RPMD is smaller than its one-bead counterparts. The difference between  from single-bead and

OD + CH4
The rate coefficients of RPMD and other theoretical and experimental values are drawn in Figure   5. The RPMD rate coefficient at 300 K matches well with the experimental value, with an error of about 17% [49]. But the rates on the whole range of temperature calculated on PES-2014 using RPMD have large deviation, which is higher than the experiment at lower temperatures but lower  well, which led to great error at the whole temperature range. The results from CVT/μOMT on PIP-NN PES display oscillation at the temperature range from 500 K to 1000 K, and show deficiency at temperatures lower than 500 K. These results reflect that VTST cannot describe the title reaction correctly. RPMD results in this work show abnormal KIE at temperature less than 700 K, and normal KIE at higher temperatures. Since there is no experimental result at temperatures larger than 700 K, we hope future experimental works confirm our finding.

OH + CD4
The isotope reaction OH + CD4 was also investigated to provide a thorough understanding of the mechanism. The rate coefficients and KIEs were collected in the lower panel of Figure 5 and However, the RPMD rates based on the PES-2014 overestimate much the results at temperature lower than 400 K. The rate coefficients and KIEs from CVT/μOMT on PIP-NN PES also near the experimental results in general, but systematically lower than experimental values, might stem from lack of recrossing in VTST methods or other error accumulation.

OH + 13 CH4
The rate coefficients of isotope of carbon, OH + 13  respectively. The average value of 1.0123 of them is close but larger than the experimental value.
The 12 C/ 13 C KIEs for OH + CH4 is a stringent test, whose results depend strongly on the treatment of the torsional anharmonicity of the lowest vibrational mode. [21] The Cayley propagator RPMD includes all vibrational degree of freedom, which correctly accounts for all the anhamonicity effects. Such great agreement shows again that the accurate rate coefficients must be obtained using both RPMD and high-fidelity PES.

Conclusion
We investigated the rate coefficients of the title reaction, OH + CH4, and three of its isotope The kinetic isotope effects of deuterated reactants OD and CD4 are also investigated. It is valuable to explore the quantum effects using KIE, but the traditional TST and its extensions will 20 lead to a huge error. Our KIE using RPMD are also in good agreement with the experimental values.
The utilizing of ultra-stable and symplectic propagator, Cayley propagator, shows considerable speed up of the calculations without accuracy loss, with deviation less than 10% comparing with the traditional propagator RPMD. With such efficiency, we obtained successfully the KIE of 12 C/ 13 C. We hope that the Cayley propagator RPMD can be further used in more complicated reactions with even larger degrees of freedom.

Author Contributions
Xiongfei Gui and Wenbin Fan contributed equally to this work.

Supporting Information
The MRP was searched using the gradient descent method, with a wider range of windows from 0.  Fig. S1 shows the MRP. Figure S1. The minimum reaction path in RPMD reaction coordinates. All energies are shifted to make the energy at 0  = zero. The black scatter-line is searched from transition state ( 1  = ) to both sides. Note that the reaction coordinates tend to slide in the negative direction of the gradient.
The red scatter-line shows the reaction coordinates interpolated from the components of configuration in each window. The subplot zoomed the transition state region with drop lines.    amu Bohr ). There are three light dashed lines for both MEP and G a V , representing the heights of reactant limit, production limit, and peak. There are three double-headed arrows with values for both MEP and G a V , representing the heights of reactants to peak, production to peak, and production to reactants. There is an additional double-headed arrow from the peak of MEP to the G a V vertically.
39 Figure S6. The PMFs and transmission coefficients of standard RPMD propagator and Cayley propagator for the OH + CH4 at 200 K, 300 K, and 500 K.