Ring Polymer Molecular Dynamics of the C(D)+H2 Reaction on the Most Recent Potential Energy Surfaces†

Ring polymer molecular dynamics (RPMD) calculations for the C(D)+H2 reaction are performed on the Zhang-Ma-Bian ab initio potential energy surfaces (PESs) recently constructed by our group, which are unique in very good descriptions of the regions around conical intersections and of van der Waals (vdW) interactions. The calculated reaction thermal rate coefficients are in very good agreement with the latest experimental results. The rate coefficients obtained from the ground ã1A′ ZMB-a PES are much larger than those from the previous RKHS PES, which can be attributed to that the vdW saddles on our PESs have very different dynamical effects from the vdW wells on the previous PESs, indicating that the RPMD approach is able to include dynamical effects of the topological structures caused by vdW interactions. The importance of the excited b̃1A′′ ZMB-b PES and quantum effects in the title reaction is also underscored.


I. INTRODUCTION
Complex-forming reactions, which are characterized by deep potential energy wells supporting intermediate complexes, are very important in many fields such as the combustion, atmospheric and interstellar chemistry. While the kinetic and dynamical pictures of direct reactions have been well established [1][2][3][4][5][6][7], the understanding of complex-forming reactions is still far from complete nowadays. Quantum dynamics (QD) studies on complex-forming reactions are very difficult due to † Part of Special Issue "John Z.H. Zhang Festschrift for celebrating his 60th birthday".
One of the prototypical complex-forming reactions is the C( 1 D)+H 2 reaction, and extensive experimental and theoretical efforts have been devoted to this reaction and its isotopic variants. Experimentally, its thermal rate coefficients were measured with different methods; the early measurements were carried out at room temperature, and recent studies focused on the temperature dependence of rate coefficients. In particular, the measured room-temperature rate coefficients reported by different groups [12][13][14] varied from (2.0±0.6)×10 −10 cm 3 ·s −1 to (3.7±0.2)×10 −10 cm 3 ·s −1 , and the latest value of (3.20±0.33)×10 −10 cm 3 ·s −1 was reported by Hickson et al. [14] using the supersonic flow reactor technique. In addition, by using a crossed-beam tech-nique, Liu [15] measured the excitation function of the C( 1 D)+D 2 reaction, from the shape of which it was concluded that little temperature dependence of rate coefficients was anticipated. The crossed-beam experiments were also used to investigate product angular and state distributions [15,16]. On the theoretical side, before our studies, three ab initio potential energy surfaces (PESs) for the groundã 1 A ′ state were constructed by other groups, denoted as the BHL [17], RKHS [18] and DMBE [19] PESs, respectively, using which various kinetic and dynamical calculations were performed [20][21][22][23][24]; as for the excitedb 1 A ′′ state, only one ab initio PES (referred to as BJHL) was reported [25]. Although the subsequent calculations on the BJHL PES displayed that the addition of the contribution of this PES worsened the agreement with experimental results [26], indicating the deficiency of this PES, a series of calculations were performed by using this PES [14,[27][28][29][30], including ring polymer molecular dynamics (RPMD) calculations [14,30].
The studies by our group [31][32][33] on the C( 1 D)+H 2 reaction and its isotopic variants can be traced back to the investigation of the PES intersection, in which Bian and coworkers [31] systematically revealed various PES intersection seams among the five lowest singlet electronic states in the C( 1 D)+H 2 reactive system, providing us a basis for further PES construction and dynamical studies. Later, a highly accurate global ab initio PES for theã 1 A ′ state, referred to as the Zhang-Ma-Bian-a (ZMB-a) surface [34], was constructed, which is unique in the very good description of the regions around conical intersections (CIs) and of van der Waals (vdW) interactions [34]. The accurate QD [35,36], quasiclassical trajectory (QCT) [37], statistical quantum mechanical (SQM) [38] and mean potential phase space theory (MPPST) [38] calculations on this surface for the C( 1 D)+H 2 reaction have yielded lots of interesting results. For example, remarkable nonstatistical behaviors in the differential cross sections were revealed. Recently, a highly accurate global ab initio PES (called ZMBb) [8,9] for the excitedb 1 A ′′ state was constructed. Accurate QD calculations on our ZMB-a and ZMB-b surfaces yielded various dynamical quantities for the C( 1 D)+D 2 reaction in very good agreement with results from crossed-beam experiments [9]. More importantly, the concept of vdW saddle in a chemical reaction was proposed, which is expected to become a useful concept in the study of complex-forming reactions; the dynami-cal importance of the vdW saddle and the excitedb 1 A ′′ PES was also underscored [9]. Most recently, a combined theoretical and experimental study [8] was performed for the C( 1 D)+HD reaction, and excellent agreement was achieved between theory and experiment for rate coefficients and product CD/CH branching ratios; meanwhile, a kind of conical intersection regulated (CI-R) intermediate was uncovered, which has the effect of bond-selective activation and leads to two unusual reaction mechanisms. It was also shown that the CI-R intermediate exists in other systems, and its dynamical effects may have general significance [8]. Moreover, various dynamical and kinetic quantities obtained from QCT and accurate QD calculations for the C( 1 D)+HD reaction were reported, which are in very good agreement with the available experimental data [39,40].
The RPMD method [41][42][43] is able to correctly and efficiently describe quantum effects such as tunneling and zero-point energy in the reaction process. A series of recent studies [44][45][46][47] have proven that the RPMD calculations can yield reliable estimates of thermal rate coefficients for a variety of reactions over wide temperature ranges. The RPMD method has been successfully applied to several complex-forming reactions such as the O( 1 D)+H 2 [45,48], S( 1 D)+H 2 [49], C( 1 D)+H 2 [14,49,50], and Si( 1 D)+H 2 reactions [44]. For the C( 1 D)+H 2 reaction, the most recent ZMB-a and ZMBb PESs [8,9,34] are much more accurate than the previous PESs, however, the RPMD calculations on these PESs are very limited, and the available RPMD calculations on ZMB-a only involve a rather narrow temperature range (200−400 K) [50]. To gain more insights into the kinetics and dynamics of the C( 1 D)+H 2 reaction, it is desirable to perform a complete RPMD study on the ZMB-a and ZMB-b PESs.
In the present work, we perform detailed RPMD calculations for the C( 1 D)+H 2 reaction on the most recent ZMB-a and ZMB-b surfaces constructed by our group. The thermal rate coefficients are calculated at temperatures ranging from 50 K to 900 K, which are compared with the latest experimental measurements. The roles of the excitedb 1 A ′′ PES and quantum effects in the title reaction are also investigated.
This paper is organized as follows. The theory and details of the calculations are given in Sec. II. The results are presented and discussed in Sec. III. Finally, a brief summary is given in Sec. IV.

A. Theory
Unlike accurate QD approaches [51][52][53][54][55] which are based on solving the nuclear Schrödinger equation, the RPMD method is based on the imaginary-time path integral formalism. The detailed explanation of the RPMD rate theory and computational procedure can be found in many publications [42,56]. Here we just outline the theory with some important formulas men-tioned in this work. The ring polymer rate coefficient is is the n-bead path integral approximation to the quantum mechanical partition function of the reactants per unit volume and c here β n =β/n with β=1/(k B T ) is the appropriate reciprocal temperature, f =3N is the total number of physical degrees of freedom, and H n (p, q) is the Hamiltonian of a system of N three-dimensional n-bead harmonic ring polymers with an external potential of V (q 1 , . . ., q N ) acting on each bead, with ω n =1/(β n ) and q is the velocity along this reaction coordinate, and h[s(q t )] is a heaviside step function which counts ring polymer trajectories (p t , q t ) that are on the product side of the dividing surface at time t.
The method begins through introducing two dividing surfaces defined in terms of the ring polymer centroid variables [57]. The first dividing surface is located in the asymptotic reactant valley where R is the centroid of the vector that connects the centers of mass of the reactant molecules and R ∞ is an adjustable parameter that is chosen to be sufficiently large as to make the interaction between the reactants negligible. The second dividing surface is located in the transition state region and is defined in terms of the bond-breaking and bond-forming distances, where N bonds is the number relevant combinations of breaking and forming bonds, and where N channel is the number of equivalent product channels and is the dividing surface for each individual combination.
Here, one bond breaks between atoms i and j, and one bond forms between atoms j and k.q ij denotes the vector that connects the centroids of atoms i and j, and q ‡ ij is the corresponding interatomic distance at the transition state saddle point. The reaction coordinate ξ is taken to be an interpolating function that connects these dividing surfaces, such that ξ→0 as s 0 →0 and ξ→1 as s 1 →0.
With the Bennett-Chandler factorization [58,59], the rate calculation has been split into two separate stages: the evaluation of the centroid-density quantum transition-state theory (QTST) rate coefficient k QTST DOI:10.1063/1674-0068/cjcp2110197 c ⃝2021 Chinese Physical Society and the evaluation of the transmission coefficient κ for a dividing surface located in the transition state region.
k QTST can be expressed in terms of the centroid potential of mean force (PMF) along the reaction coordinate, W (ξ ‡ ), as Refs. [42,57] is the dynamical correction accounting for recrossing of the dividing surface at t→t p where t p is a "plateau" time [59], which ensures that the resulting RPMD rate coefficient k RPMD (T ) will be independent of the choice of dividing surfaces [60]. It is given by [42] where and the angular brackets in Eq.(11) denote a canonical average.

B. Calculation details
The present RPMD calculations are performed using the RPMD rate code [56], and the ZMB PESs [8,9,34] constructed by our group are used. The simulation parameters are given in Table I. The dividing surfaces are carefully defined in the entrance channel. For a given temperature, test calculations are performed to find the reaction coordinate position of the free-energy barrier, then the dividing surface is defined with the geometry corresponding this reaction coordinate position to ensure that the free-energy barrier is located around reaction coordinate ξ=1. Convergence has been also carefully checked by performing RPMD calculations with different numbers of beads. It should be noted that the choice of dividing surfaces in this work is different from that in the previous RPMD calculations [14,49], in which the second dividing surface was placed at the potential wells. As shown in FIG. 1, on the RKHS PES, the RPMD rate coefficients with our choice of dividing surfaces coincide well with the previous RPMD results [14], proving that the choice of dividing surfaces has negligible effects on RPMD rate coefficients. In addition, the electronic states of C( 1 D)+H 2 are fivefold degenerate and thus the electronic degeneracy factor of 0.2 has been taken into account in the rate coefficients calculations.

III. RESULTS AND DISCUSSION
The convergences of k QTST , the plateau value of ring polymer transmission coefficients (κ plateau ), and k RPMD with respect to the number of beads (N bead ) on the ZMB-a surface for the title reaction are shown in FIG. 2. As can be seen, the rate coefficients with N bead =16 are somewhat close to those with N bead ≥32 but much larger than those with N bead =1, and the difference between RPMD results with N bead =32 and those with N bead =64 is very small, meaning that RPMD calculations at N bead =32 can yield converged results. In addition, the role of quantum effects in chemical reactions and molecular systems is an issue of considerable interest [61][62][63][64], the RPMD rate coefficient with N bead =1 is close to the classical limit, which provides us a good chance to evaluate the contribution of quantum effects to the reactivity by comparing the results with N bead =1 and those with N bead =64. We find that quantum effects are very important in the title reaction. For instance, quantum effects account for about 34% of the reactivity at 50 K, and the proportion decreases to about 17% at 500 K, meaning that quantum effects are more impor-   reaction in which the PMF curve displays two barriers [44]. The free-energy barrier height is very small, for example, at 300 K the free-energy barrier is 0.014 eV on the ZMB-a surface, and the corresponding value on the ZMB-b surface is 0.019 eV. On both the ZMB-a and ZMB-b PESs, the free-energy barrier height increases with the increasing temperature, and the shapes of the PMF curve do not change much with the temperatures, which is very similar to the behaviors displayed in the O( 1 D)+H 2 and S( 1 D)+H 2 reactions [45,49]. In addition, we can see that each PMF curve in FIG. 3 shows a sharp drop when the maximum of the PMF value is reached, thus the location of the free-energy barrier is close to the PES regions around the entrance-channel vdW saddle, indicating that the vdW saddles on the ZMB-a and ZMB-b surfaces play an important role in the dynamics of the title reaction.
The transmission coefficients for the C( 1 D)+H 2 re- action are plotted in FIG. 4. As can be seen, on both the ZMB-a and ZMB-b surfaces, the curves of transmission coefficients reach the plateau value after a long propagation time, which is consistent with the complexforming mechanism. In addition, the propagation time decreases with the increasing temperature. For example, on the ZMB-a surface, the transmission coefficient takes about 1300 fs to converge at the temperature of 50 K, whereas the corresponding value at 900 K is around 500 fs. This is understandable in light of the fact that the intermediate complex is on average much more long-lived at low temperatures. It seems that the difference of propagation time between the ZMB-a and ZMB-b surfaces is negligible at temperatures larger than 500 K, whereas the propagation time on ZMB-b is a bit larger than that on ZMB-a at low temperatures. This may be attributed to the dynamical effects of the CI-R intermediate on the ZMB-b surface. As demonstrated in Ref. [8], the CI between theb 1 A ′′ andd 1 A ′′ states induces a small barrier and a shallow precursor well (which supports a CI-R intermediate) before the deep well on the ZMB-b PES, and obviously such complex topological structure needs more time to converge the transmission coefficient. In addition, the values of κ plateau at different temperatures are listed in Table II. It can be seen that κ plateau roughly drops with the increase of temperature. This negative temperature dependence of κ plateau is consistent with those observed previously in other complex-forming reactions [44,65].
The RPMD rate coefficients obtained from the ZMBa PES in the temperature range of 50−300 K are shown    5, compared with the previous RPMD results from the RKHS surface. As can be seen, although both rate coefficients obtained from ZMB-a and RKHS show slight positive temperature dependence, our rate coefficients from ZMB-a are much larger than those from RKHS in the whole temperature range. In particular, the rate coefficient of the present work at 300 K is 1.95×10 −10 cm 3 ·s −1 , whereas the corresponding value from RKHS is 1.6×10 −10 cm 3 ·s −1 . This phenomenon results from the dynamical effects of the vdW saddle on the ZMB-a PES, as shown in our previous work [9]. FIG. 5(B) displays the values of ∆k, which is defined as the difference between the rate coefficient from ZMB-a and that from RKHS (i.e., k ZMB-a −k RKHS ). As can be seen, the ∆k values obtained from the RPMD calculations are in good accord with those obtained from the accurate QD calculations, and the reproduction of such  [14], the quantum dynamics (QD) results on the ZMB-a surface are from Ref. [35], and the corresponding QD results on the RKHS surface are from Ref. [66]. behavior by means of the RPMD calculation indicates that the RPMD method is able to include dynamical effects of the topological structures caused by vdW interactions.
FIG . 6 shows the temperature dependence of thermal rate coefficients (summing the results from ZMB-a and ZMB-b) for the C( 1 D)+H 2 reaction in the temperature range of 50−900 K, compared with the latest experimental results by Hickson's group [14]. As can be seen, the rate coefficients are nearly temperatureindependent in the studied temperature range, and the agreement between the present RPMD results and experimental values is very good. In particular, the thermal rate coefficient of the present work at room temperature is 3.44×10 −10 cm 3 ·s −1 , which is within the error bar of Hickson's experimental value of (3.20±0.33)×10 −10 cm 3 ·s −1 , while the previous RPMD calculations [14] reported a lower value of 2.8×10 −10 cm 3 ·s −1 at the same temperature. Moreover, the present RPMD rate coefficient shows a slight positive temperature dependence below 300 K, which well reproduces the behavior of experimental values. More interestingly, in the present RPMD calculations, the excited b 1 A ′′ state makes a contribution of around 43% to the overall reactivity. Considering that theb 1 A ′′ state contribution to the overall reactivity is about 38% in our previous QCT and QD calculations on the ZMB-a and ZMB-b surfaces, it seems that the present RPMD calculations slightly overestimate the contribution of thẽ b 1 A ′′ ZMB-b PES, which may be due to the complex topological structures in the regions around the CI-R intermediate on the ZMB-b PES [8]. Nevertheless, the present RPMD calculations underscored the important role of the excitedb 1 A ′′ ZMB-b PES in the title reaction.

IV. CONCLUSION
In this work, the kinetics and dynamics of the typical complex-forming C( 1 D)+H 2 reaction are investigated, which is achieved by detailed RPMD calculations on the most recent ZMB-a and ZMB-b PESs constructed by our group. The interesting features in PMF curves and transmission coefficients are displayed and discussed, and the contribution of quantum effects to the reactivity is also evaluated. More importantly, the obtained RPMD thermal rate coefficients are in very good agreement with the latest experimental results, and the importance of the excitedb 1 A ′′ ZMB-b PES is underscored. Last but not least, the present RPMD calculations also indicate that the RPMD method is able to include dynamical effects of the topological structures caused by vdW interactions. We hope that the present work will stimulate further research interests in the title reaction.