Reaction dynamics for the Cl($^2$P) + XCl $\to$ XCl + Cl($^2$P) (X = H, D, Mu) reaction on a high-fidelity ground state potential energy surface

Globally accurate full-dimensional ground state potential energy surface (PES) for the Cl($^2$P) + XCl $\to$ HCl + Cl($^2$P) reaction, a prototypical heavy-light-heavy abstract reaction, is developed using permutation invariant polynomial neural network (PIP-NN) method and embedded atom neural network (EANN) method, with the corresponding total root mean square error (RMSE) being only 0.043 and 0.056 kcal/mol, respectively. The saddle point of this reaction system is found to be nonlinear. A full-dimensional approximate quantum mechanical method, ring-polymer molecular dynamics (RPMD) with Cayley propagator, is employed to calculate the thermal rate coefficients and kinetic isotopic effects of title reactions Cl($^2$P) + XCl $\to$ XCl + Cl($^2$P) (X = H, D, Mu) on both new PESs. The results reproduce the experimental results at high temperatures perfectly, but with moderate accuracy at lower temperatures. The similar kinetic behavior is supported by quantum dynamics using wave packet calculations as well.

as 8.3 kcal/mol. It's found that BCMR PES gave highly rotational excitation products, by studies of quasi-classical trajectory (QCT) and quantum scattering dynamics. 7 Later there are a pair of PESs based on ab initio calculations are released. The first one is from D. Truhlar's group. 1 was fitted by rotated-Morse-oscillator-spline (RMOS) based on ~5500 ab initio points using the polarization configuration-interaction (POL-CI) method. 9 Although the TS from POL-CI calculation was not collinear, it was set to be collinear on the PES, with a potential barrier as 10 shell coupled cluster singles doubles with perturbation triples (RCCSD-T) 10 and 3 multireference configuration interactions (MRCI), 11,12 then fitted by the rotated-Morse cubic-spline function. The DCBKS PES contains three electronic states, and its TS is nonlinear. This nonlinear feature was attributed as the repulsion of both p-polarized orbitals on Cl atoms. But the potential barrier of TS was also scaled by a factor of 0.815 to match the calculated rate coefficients from QCT to the experimental values. 13 Recently, our group 14 calculated the thermal rate coefficients and KIE using an approximate full-dimensional quantum mechanical method named ring polymer molecular dynamics (RPMD) method on LEPS PES 15,16 . The results of RPMD are also consistent with those of other theoretical approaches such as ICVT and quantum dynamics (QD). For Cl( 2 P) + DCl reactions, the RPMD rate coefficients at higher temperatures are very accurate compared with the experimental results. However at 312.5 K the RPMD results are slightly lower than the experimental values, although it still close to all values from other theoretical methods. This may stem from the inaccuracy of the LEPS PES, since during our previous experience, the quality of PES used is essential for RPMD calculations. And since both the geometry and energetics of TS are not well defined from above discussion, it's needed to prepare accurate PES from high-level quantum chemical method.
To unveil the dynamics of the title reaction, it is usually essential to build an accurate global PES, which can be achieved by fitting a large number of high-level ab initio energy points. According to the previous work of G. Schatz 8 , the MRCI is necessary to calculate the energy of sampled configurations, due to the multiple electronic states features of the reaction system. And for fitting the PES, recently there is a novel method named embedded atom neural network (EANN) proposed by Bin Jiang's group 17 . It is extendable for high dimensional bimolecular reactions when with active learning technique. So, we choose the EANN in this work. To test the performance of the EANN method on constructing potential energy surfaces of bimolecular reactions, we have also adopted the standard permutation invariant polynomial-neural network (PIP-NN) method [18][19][20] , which has already been demonstrated to be suitable in fitting polyatomic bimolecular reactions and widely used, such as OH + H2O, 21  for isotope H is also validated by quantum dynamics using wave packets. This work is organized as follows. The PESs employed in the current work and the related theories and calculation details are introduced in section II. The results are presented and discussed in Section III. The final summary is contained in Section IV.

II.A PIP-NN PES
All electronic calculations in this work were performed using MOLPRO2015. 27 The geometries, energies and harmonic frequencies of all stationary points of the reactants, products and TSs were obtained at the MRCI-F12+Q levels. [28][29][30] The method is proved to give reliable potential energy surfaces since the higher energy accuracy.
The initial data set is initially sampled at three bond lengths R HCl 1 , R HCl 2 , and R Cl 1 Cl 2 in the range of 0.8-20 Å. Further improvement of the potential energy surface by adding points to the area around the stationary points and to the reaction path. Therefore, we obtain a primitive potential energy surface.
Based on this primitive potential energy surface, RPMD calculates from 300-1500 K for different temperatures. Exploration of dynamically relevant regions verifies the performance of the potential energy surface, which behaves unreliably in regions with lacking points. Therefore, data points in these regions are sampled to repair the potential energy surface. This process is repeated over and over again to improve the potential energy surface until all relevant dynamical results converge. To improve sampling efficiency, only those points that are not close to the existing dataset are added, using the generalized Euclidean distance

II.B Embedded Atom Neural Network Potentials
Although the PIPNN method works well for constructing potential energy surfaces, it is difficult to extend to the large system containing many atoms, which is caused by too many PIPs. Thus, Using the embedded atom neural network (EANN) method, 17 one can construct a surface of high-dimensional potential energy in which the total energy is calculated based on atomic energies. Specifically, atoms are embedded in an environment of other atoms, and their atomic energy is derived from an atomic neural network based on nonlinear transformations of electron densities, In Eq. 2, i ρ is a density-like structural descriptor that may be constructed simply by atomic orbitals of Gaussian type centered around neighboring atoms, where ij r is the Cartesian coordinate vector of embedded atom i with respect to neighbouring atom j, and r is its norm;

II.C RPMD reaction rate theory
All calculations are performed using the RPMD rate theory implemented in the RPMDrate code 38 . Since there have been lots of review articles about it, 39, 40 here we only give a brief summary related closely to the current work. For the title reaction, this Hamiltonian can be written as follows: Where ˆi p , ˆi q and i m are the momentum operator, the position operator and the 7 mass of the ith atom, respectively. Taking advantage of classical isomorphism between quantum systems and the ring polymer, each quantum particle is represented by a necklace formed by n classical beads connected by a harmonic potential: 15, 41 Where for each atom, Here QSTS ( ; ) kT   is the centroid-density QTST rate coefficient, 16, 44 evaluated at the maximum of the free energy barrier,   , along the reaction coordinate ()  q . In practice, it is calculated from the centroid potential of the mean force (PMF): 38-40 Where R  is the reduced mass between the two reactants, is the freeenergy difference which is obtained via umbrella integration along the reaction coordinate. 38,45 The second factor ( ; ) t   → is named the transmission coefficient, which provides dynamical correction and is calculated by the ratio between long-time limit and zero-time limit of the flux-side correlation function: Which captures the recrossing of the TS region and ensures that the obtained RPMD rate coefficient results do not depend on the choice of the dividing surface. 16 It should be noted that the final RPMD rate coefficients are corrected by an electronic partition function ratio of the following form: to account for the spin-orbit splitting of In addition, when only one bead is used, the results from RPMD will reduce to the classical limit. In this limit, the static and dynamic components become the same as the classical transition state theory (TST) rate coefficient and the classical transmission coefficient, respectively. Therefore, these quantities determine the limits at which quantum effects, such as ZPE and tunneling effects, can be evaluated by using more beads. The minimum number of beads considering the quantum effect can be given by the following formula: 47 min max n  = (11) Where max  is the largest vibration frequency in the system. In this work, the convergence is tested with increasing number of beads, and the numbers that yield converged results of PMF are chosen at different temperatures. In the Supporting Information, Figure S2 shows the convergence of PMF curves at 312.5 K obtained from different numbers of beads.
Additionally, there is a critical temperature named cross-over temperature 48 , Tc, which also needs to be considered: Where b i is the imaginary frequency of the reaction system in the TS. The reaction system temperature is lower than T c , which is considered as deep-tunneling region, the error of RPMD results would become large. Enough beads are needed to obtain accurate results. The cross-over temperature for the title reaction is T c =348 K.

II.D Quantum dynamics
The time-dependent quantum wave packet method is used to calculate the thermal rate coefficients as well. The Hamiltonian of the system in the reactant Jacobi coordinates can be written as 53 where R  is the reduced mass between Cl and HCl; r  is the reduced mass of the reactant HCl; R is the distance from the attacking Cl atom to the center of mass of the reactant HCl; r is the bond distance of the reactant HCl; θ is the bending angle between the vectors R and r;ˆt ot J is the total angular momentum operator of the system; ĵ is the rotational angular momentum operator of the reactant HCl; V is the potential energy operator.
The time-dependent wavefunction is expanded as where () The initial state-specific rate constant is obtained by thermal averaging the collision energy of the corresponding ICS as where Ec is the collision energy and kB is the Boltzmann constant. The electronic partition function Qe is given by 2 +

11
An L-shaped grid is used in this work. 55 The numerical parameters are listed in Table I.   Figure 2 shows contour plots as functions of the breaking (RHCl 1 ) and forming (RHCl 2 ) bonds with the bond angle ∠ClHCl relaxed.

III.A Properties of the NN PES
It is clear to know that this is a typical symmetric reaction. However, the real frequency difference of DCBKS PES is not much, but its virtual frequency difference is up to 10%.   [3] f The DIM-3C PES, see detail in Ref. [6] g The PK3 PES, see detail in Ref. [5] h The POLCl PES, see detail in Ref. [1] i The DCBKS PES, see detail in Ref. [8]

III.B RPMD rate coefficients
In this work, the thermal rate coefficients for the Cl( 2 P) + XCl (X=H, D, Mu) reactions were calculated at the temperature range of 200-1000 K. This calculation is first performed with one bead, which provides the classical limit, and then the number of beads increases until converged.
In Cayley-RPMD calculations the number of beads needed depends on different temperatures and isotopes. The number of beads should not be less than the minimum suggested by Eq. (11). The rate coefficients calculated by the RPMD method and other theoretical methods as well as the rate coefficients measured experimentally are listed in Table II.
As described in the section II C, the reaction is in deep-tunneling region at 312.5 K, since the Tc is 348 K for the title reaction. From previous discussions, the results below Tc would underestimate the rate coefficients since the title reaction is with symmetric barrier. 56 But in this work, the Cayley-RPMD results are still in good agreement with experimental values.
The left panel of Figure 3 shows the PMF of Cl( 2 P) + XCl (X=H, D, Mu) reaction at 312.5 K. The convergent RPMD barriers (with the optimal number of beads) for all three reactions are lower than the classical (single-bead) results. This is due to tunneling effects that make it easier for the three isotopes to penetrate the potential energy barrier.
The free energy barriers from different isotopic reactions decrease as ∆GD > ∆GH > ∆GMu according to the decrease in the mass of these isotopes. This order comes from the fact that the smaller the mass, the greater the tunneling capacity. The right panel of Figure 3 shows  Table II shows that the converged RPMD transmission coefficients also increase with the increasing of temperature, indicating that the recrossing is more at low temperature, which is consistent with previous studies 39,40 . As can be seen from the table II, the potential of mean force (PMF) of RPMD rates at temperatures 200 K to 1000 K, the free energy barrier increases with the increase in temperature, which is  The left panel of Figure 4 shows  Therefore, the states with simultaneous vibration-rotation excitation are not considered in the quantum wave packet calculations, which is believed to have negligible effect on the thermal rate constants. As shown in Figure 4, the quantum thermal rate constants agree reasonably well with the Cayley-RPMD results although the former are slightly higher at low temperatures. QD calculations show that, on one hand, the reaction energy threshold is visibly lower than the classical barrier height, indicating the existence of significant quantum tunneling effect for the reactions, which consistent with the PMF curves from RPMD, as shown in Figure   The right panels of Figure 4 shows Cl( 2 P) + DCl rate coefficients in the Arrhenius

IV. Summary and conclusions
In this work, we have investigate the reaction dynamics of Cl( 2 P) + HCl → HCl + Cl( 2 P). Firstly, we have developed two new full-dimensional neural network PESs for ground state of the title reaction, the permutation invariant polynomial neural network (PIP-NN) method and embedded atomic neural network (EANN) method, based on 5986 points and 6515 points at MRCI-F12+Q/AVTZ level, and with the total root mean square error was 0.043 kcal/mol and 0.056 kcal/mol, respectively. In particular, the EANN method is the first application of gas phase bimolecular reactions.
From the comparison, the performance is well, and is promising for being applied into multiple atomic reactions with active learning technique. We have also confirmed the