Initial state-selected scattering for the reactions H+CH 4 /CHD 3 and F+CHD 3 employing ring polymer molecular dynamics.

The inclusion of nuclear quantum eﬀects (NQEs) in molecular dynamics simulations is one of the major obstacles for an accurate modeling of molecular scattering processes involving more than a couple of atoms. An eﬃcient method to incorporate these eﬀects is ring polymer molecular dynamics (RPMD). Here we extend the scope of our recently developed method based on non-equilibrium RPMD (NE-RPMD) from triatomic chemical reactions, to reactions involving more atoms. We test the robustness and accuracy of the method by computing the integral cross sections (ICS) for the H/F+CH 4 /CHD 3 reactions where the methane molecule is either initially in its vibrational ground or excited state (C-H stretch). Furthermore, we analyze the extent to which NQEs are described by NE-RPMD. The method shows signiﬁcant improvement over the quasiclassical trajectory (QCT) approach while remaining computationally eﬃcient.


I. INTRODUCTION
One of the challenges of physical chemistry is to study in details the dynamics of chemical reactions. The detailed understanding of the underlying mechanisms is important for many domains of research such as atmospheric 1-3 and interstellar chemistry 4-6 , combustion 7 and catalysis 8 . A key step in this direction of study is to develop accurate reactive scattering methods in order to qualitatively and quantitatively assess the influence of the reactants' rovibrational state on the reaction dynamics. 9,10 In this context, nuclear quantum effects (NQEs) can have great impact on the reaction dynamics. [11][12][13][14][15][16] Unfortunately, initial state-selected full-dimensional quantum scattering simulations that incorporate NQEs are only computationally feasible for systems up to 6 atoms with a vanishing angular momentum, although significant development in this direction has been achieved. [17][18][19][20][21][22][23][24][25] Consequently, there is a demand to develop computationally efficient simulation methods which take NQEs into account, at least approximately.
One possibility is to employ quantum simulations with constraints on the dynamics 20,21,[26][27][28][29][30][31][32][33][34][35][36][37][38][39] to reduce the number of degrees of freedom and make the simulations amenable. Those methods are called reduced dimensional quantum dynamics simulations (RDQD). However, dynamical constraints and their potential artifacts have to be elaborated for each system specifically and it is rather challenging to assess their accuracy in describing the dynamics. This renders the RDQD approach non-systematic in contrast to full-dimensional simulations. Moreover, employing effective RDQD models to a) Electronic mail: adrien.marjollet@desy.de b) Electronic mail: ludger.inhester@desy.de c) Electronic mail: ralph.welsch@gmx.de simulate systems such as the title reactions are nonetheless very expensive and can take many months to be completed. 21,40 Another very popular alternative is the so-called quasi classical trajectory method (QCT) 9, [41][42][43][44] . It circumvents the very high computational cost of quantum simulations by rendering the simulations classical, thus neglecting crucial aspects of NQEs in the dynamics. The atoms are propagated with the classical equations of motion on a potential energy surface (PES) employing initial conditions that mimic quantum rotational-vibrational states. This approach has successfully revealed insights into microscopic details of chemical reaction dynamics [45][46][47][48][49] . QCT trajectories are often easier to interpret compared to RDQD trajectories. One of the main shortcomings of the QCT approach is the ZPE leakage problem. It originates from the fact that in classical dynamics the zero-point energy (ZPE) contained in the vibrational modes of the molecules typically flows from high to low frequency modes during the simulation. Consequently, parts of the vibrational energy can flow into the reaction coordinates, leading to an artificial enhancement of the reactivity. [50][51][52] Moreover, the formation of products having less than their quantum ZPE value can occur, violating the ZPE constraints. 15,[52][53][54][55][56] In addition, tunneling effects are also absent in classical simulations. The neglect of these NQEs can lead to several inaccuracies when computing energy thresholds, reaction probabilities, and integral cross sections.
A relatively recent approach that combines the advantage of being computationally efficient and at the same time takes NQEs into account is ring polymer molecular dynamics (RPMD) [57][58][59][60][61][62][63][64][65][66] . It has been shown that RPMD is able to accurately compute thermal rate constants for many reactions 64,[67][68][69][70] . More recently, RPMD was also successfully applied to the calculations of Kubotransformed correlation functions for non-equilibrium initial conditions (NE-RPMD) 66 . RPMD and NE-RPMD have been theoretically justified as approximations of Matsubara dynamics. 65,66,71 Various results for applications of RPMD beyond the thermal equilibrium ensemble followed. [72][73][74][75] In our earlier works 76,77 , two of us developed a variant of NE-RPMD and applied it to reactive collisions of various atoms with the hydrogen molecule H 2 . Here, we extend the application of this method to more complex reactions.
To test newly developed methods, the computation of integral cross section (ICS) constitute reliable and precise indicators for the accuracy of the simulations. The ICS can be measured via challenging experiments and its variations with different isotopes and initial rovibrational states reveal important details about the reaction dynamics. 49,[78][79][80] The elementary reaction H + CH 4 → H 2 + CH 3 together with its isotopic variants play an important role in combustion chemistry, interstellar chemistry and has been the subject of extensive experimental [80][81][82] and theoretical studies 20,25,36,37,[83][84][85][86][87][88][89][90][91][92][93][94] . This reaction involves a small number of electrons and potential energy surfaces (PES) with quantitative level of accuracy were constructed via extensive high level ab initio calculations. [90][91][92][93] QCT simulations typically overestimates the reactivity for these reactions as classical dynamics ignore ZPE constraints and allow the ZPE to flow to the reaction coordinates. 15,55 Also, the presence of tunneling effects due to the light H atoms are not described in QCT. Various RDQD methods applied to tackle six-atom systems have been developed to study the reaction. [29][30][31][32][33] In particular, it was shown that restricting the non-reacting CH 3 group under C 3v symmetry yields a quantitative level of accuracy for the computation of ICS. 27,28 Also, the strong reactivity enhancing effects of the vibrational excitation of the C-H stretch in CH 4 and CHD 3 have been studied theoretically 21,24,29,31,92,[95][96][97][98][99][100] and experimentally 80,82 . Therefore, H + CH 4 and its isotopic variants constitute valuable benchmark hexatomic reactions for NE-RPMD.
Another extensively studied reaction is F + CHD 3 → FH + CD 3 . 14, 34,46,47,[101][102][103][104][105][106][107][108][109][110] The most accurate PES (PWEM-SO) for this system to date has been developed by Palma, Westermann, Eisfeld, and Manthe by considering the spin-orbit coupling in the reactant channel. 109 This intricate reaction has a very low barrier ( 0.8 kcal/mol) and its dynamics involve reactive resonances and are strongly affected by stereodynamics forces at low collision energies. 34,109 Morever, mode specificity of the reaction dynamics sparked controversies between theoretical models and experimental results regarding the effects of the C-H excited stretch on the reactivity. 103,105,111 More precisely, recent experimental studies 105,111 of the F + CHD 3 reaction found that the excitation of the C-H vibration hinders the overall reactivity while both recent QCT 47,109 and RDQD 34 simulations predict an overall enhancement of the reactivity. These discrepancies are still unexplained and several speculations exist that ascribe them to shortcomings of the experiments 34,112 or the theoretical models 34,35,109 (QCT, RDQD or the PES). NE-RPMD is capable of approximating NQEs with full-dimensional calculations, thus having the advantages of both QCT and RDQD to some extent. As such, comparing its ICS results with those of QCT and RDQD for this reaction constitutes a good test to explore the capabilities of NE-RPMD. Therefore, this reaction is well suited to benchmark methods capable of approximating NQEs in full-dimensional calculations such as NE-RPMD.
In this paper we compare the accuracy of ICS computed with the NE-RPMD, QCT, and RDQD methods for the aforementioned reactions. We compute initial state-resolved ICS for the reactive collisions H + CH 4 , H + CHD 3 , and F + CHD 3 . We analyze the impact of NQEs for the reactivity and report to what degree they are incorporated in the respective simulations. We show that the NE-RPMD approach is a promising and computationally efficient method to obtain stateresolved insights for reactive collisions involving complex polyatomic molecules.
The structure of the paper is as follows: In section II, we briefly outline the ideas behind QCT and RPMD and give references for further details. Section III explains in details the NE-RPMD approach to obtain the presented results. Section IV consists of an analysis of the obtained ICS results for the title reactions with a discussion on the NQEs described by the approach. Finally, we summarize and give our conclusions in Sec. V.

II. THEORY
In this section, we discuss the QCT method, RPMD (alongside NE-RPMD) and the NE-RPMD approach we developed, which includes aspects of QCT used in the context of NE-RPMD. Both the QCT and RPMD methods have been widely reviewed 9,41,42,57,58,65,66 , so we will only discuss the main aspects related to the NE-RPMD approach. Details of NE-RPMD applied to triatomic systems are given in Ref. 77, thus we only summarize the essential steps here.

A. QCT
In QCT, the atoms are modelled as classical pointlike masses. A standard normal-mode sampling 113 with semiclassical quantization conditions for each vibrational mode is performed to mimic a given initial vibrational state for the reactants. Specifically, positions and momenta are initialized such that the energy in each vibrational mode of the molecule is equal to the corresponding energy in the quantum state employing the harmonic approximation of the potential energy surface. The molecular vibrations are then propagated either by directly or adiabatically switching to the full (i.e. non-harmonic) potential. The rotational motion of each reactant is added using an iterative procedure as to limit spurious ro-vibrational couplings. 43,114 . The relative motion between the reactants is then set according to a fixed collision energy E col in the center-of-mass frame. The atoms are then propagated using Newton's equations of motion. This approach is very efficient, systematic, and enables feasible simulations for systems with many atoms.
There exist various effective techniques to incorporate quantum effects in the dynamics, such as discarding trajectories violating ZPE constraints or partially preventing ZPE leakage. However, these techniques remain intrinsically ad-hoc. 55,[115][116][117]

B. RPMD
The RPMD approach originally emerged from the isomorphism bridging the equilibrium quantum partition function Z to the classical partition function of fictitious ring polymers Z n . 57 Each ring polymer consists of n classical replicas (in the following called beads) joined by harmonic springs. The corresponding ring polymer Hamiltonian H n is derived from the expression of Z n as where N is the number of atoms, q ), V is the molecular potential, m i is the ith atomic mass that is used for each bead of the corresponding atom, and ω n = n β is the common frequency of each harmonic springs joining neighboring beads. 57 For equilibrium conditions, the ring polymer phase space is sampled according to the ring polymer Boltzmann distribution associated to H n at temperature n β . This is done either via thermostatted propagation in imaginary time or by employing Monte Carlo techniques. 118 The properly equilibrated ring polymers then allow for computation of exact thermal quantum observables. For example, the thermal averaged energy can be computed by where . refers to averaging over the canonical ensemble and E n is referred as the primitive energy estimator for a finite number of beads. Note the different sign in front of the harmonic spring terms of the ring polymer in Eq. 2 in contrast to the ring-polymer Hamiltonian (Eq. 1). ZPE effects and tunneling effects are incorporated through the ring polymer internal structure, which is parametrized by β and n. 68 ZPE leakage is therefore avoided, and tunneling effects are also mimicked due to the finite spatial extension of the ring polymers. 119 In 2014, RPMD was derived as an approximation to Matsubara dynamics in the context of time-dependent Kubo-transformed correlation functions. 65,71 A bit later in 2016, encouraging results for applications of RPMD to non-equilibrium time-dependent correlation functions, for the specific NE cases of a momentum impulse (or "kick") and a sudden mild shift of potential, were successfully carried out. Non-equilibrium initial conditions employed in RPMD (NE-RPMD) can also be derived from Matsubara dynamics for the case of a sudden change of potential or an initial centroid momentum "kick". 66 The extent and the conditions for which RPMD or NE-RPMD can describe real-time dynamics is not known precisely. 65,66,71 Here we employ non-equilibrium initial conditions by "kicking" the reactants centroid momenta to initialize the collision. This approach is also employed for preparing initial conditions with specific vibrational mode excitation, as described in the following section in details.

A. Reactant initialization
The two reactants in the following are labelled as "T" for the target molecule (i.e. CH 4 or CHD 3 ) and "X" for the colliding attacker atom (H or F). The entire system consists of N atoms, and the target possesses 3(N − 1) degrees of freedom. Since T is a non-linear molecule, its number of internal vibrational modes is N v = 3(N − 1) − 6. The PES of the target molecule alone placed far away from X is referred to as V T and is a scalar function of the molecule's N v internal coordinates. The initial momentum and position of each bead are first sampled according to the harmonic approximation of V T . To that end, the target bead coordinates are transformed into the N v mass-weighted vibrational mode coordinates of the total potential (including the molecular potential and the ring polymer harmonic springs potential). The resulting Hamiltonian of the ring polymer for the target molecule H where are the ring polymer mass-weighted position and momentum vibrational mode coordinates and where C and L −1 are the transformation matrices from the bead representation to the ring polymer (internal) normal mode representation and the vibrational mode coordinates, respectively. The coefficients of the matrix C are chosen such that the bead coordinates corresponding to the k = 0 index are computed with the coefficient C a 0 = 1 √ n , i.e., the bead coordinatesQ are proportional to the normal mode coordinate averaged over all beads (centroid normal mode coordinate). 120 The N v internal vibrational modes are sampled with the ring polymer Boltzmann distribution corresponding to H (T) n . Provided that n is large enough compared to β max i (Ω i ) and that β is not too low (this point is discussed in details later in the manuscript), the ring polymer ensemble is equivalent to the quantum partition function for the vibrational ground state of T. To model the excited state for the ith vibrational mode associated to vibrational quantum number ν i , we add the right amount of vibrational energy by kicking the ring polymer and shifting the coordinates along a specific vibrational mode coordinate, analogously to the procedure in QCT. To that end, we shift the centroid bead coordinates and momenta by where The resulting thermal averaged energy of T becomes where is an extra contribution to the internal energy above the zero-temperature vibrational energy. Those contributions stem from the finite values of β. Since we are aiming at pure vibrational quantum states, β will be chosen sufficiently large to keep δE th small. The ring polymer attacker associated with atom X is initially placed far away from T and its internal ring polymer normal modes are sampled according to the free ring polymer Boltzmann distribution.

B. Collision initialization
The average momentum of the ring polymer atom X is and the center of mass momentum of the ring polymer molecule T is To prepare the initial conditions for the collision, we shift the bead's momenta of X and T such that and the relative collision momentum is where µ X,T = MXMT MX+MT is the reduced mass of the system and M X and M T are the total mass of X and T, respectively. Moreover, the position of the two reaction partners are shifted orthogonal to the axis connecting X and T by a given impact parameter b. The entire system is then propagated employing a modified velocity Verlet integration scheme. 120 No thermostat is attached to the system during the propagation. The origin of time t = 0 coincides in the simulations with the momentum "kick" of the reactants at collision initialization.

C. Choice of β
The spatial extension and dynamics of the ring polymers strongly depend on the spring constant ω n = n β .
While β is naturally defined as the reciprocal temperature for thermal equilibrium applications of RPMD, the status of β in the present application of NE-RPMD is unclear. 66 Previously, we introduced an ansatz for β in which β was chosen based on the average of the reciprocal kinetic energies of the two reactants. 77 We also tweaked the ansatz so that β would be higher in the presence of strong ZPE constraints for the products and, furthermore, we introduced the lower limit β − to keep the thermal energy contributions δE th in Eq. (8) low.
Here, we employ a somewhat simpler approach and set β equal to the previously employed lower limit β − . More specifically, β is given by constraining δE th to be 5% of the total energy of the system, i.e., β solves the equation Consequently, β only depends on the polyatomic reactants' initial vibrational energy and the collision energy E col . The values for β for the different target molecules used here are shown in Fig. 1. As can be seen, larger collision energies result in lower values for β. At larger collision energies, the resulting ring polymer becomes therefore spatially more confined. This behavior reflects the expected tendency that in general, for higher collision energies, tunneling effects become relatively less important. We note that other NQEs such as ZPE constraints can still be strongly present in the dynamics for certain reactions also at high collision energies and that these effects are still captured with NE-RPMD also at low β. 12,77 FIG. 1. β values for each reactive system treated.

D. Calculation of integral cross sections
In our simulations, we propagate the system for a sufficiently long time after the collision such that the reaction products can be distinguished. N trajectories are simulated with different impact parameters b that has been sampled such that the squared impact parameter b 2 is uniformly distributed between 0 and b max , where b max is the maximal value beyond which no reaction occurs.
From the set of QCT trajectories, we calculate the integral cross section as where the vector ν T refers to the initial vibrational quantum numbers of T, and I R is a function taking as argument the positions q (i) (t) of the atoms of the trajectory i at time t and returns 1 if a reaction occurred or 0 otherwise. 41 It is important to note that QCT trajectories can lead to reaction products with internal energies below their ZPE. Here, the QCT ICS results are computed by considering all product outcomes after collision regardless of any ZPE violations. This approach, namely the histogram binning method, is employed as we avoid using ad-hoc techniques 115,116 to cure ZPE violations. This way, we are able to pinpoint which NQEs are described by NE-RPMD when comparing with the QCT.
Analogously, for the NE-RPMD simulations, the ICS are computed as follows where q (i) (t) are the centroid coordinates for trajectory i at time t. Notably, in QCT all trajectories have a fixed energy as defined by ν T and E col , whereas in NE-RPMD, the initial sampling of the ring polymers lead to energy estimator values that follow a distribution with finite width. 121 To compute ICS with NE-RPMD, we proceed by associating the outcome of each trajectory (reaction or not) to the vibrational state associated to the ν i 's in Eq. (8), which define the sampling of the ring polymer's internal vibrational modes.

A. System details and computational details
The PWEM-SO 109 and PIP-NN 94 potential energy surfaces are employed for the reactions F + CHD 3 and H + CH 4 /CHD 3 , respectively. We do not consider any non-adiabatic effects or spin-orbit couplings as we are strictly interested in testing the capability of the NE-RPMD method for describing NQEs in the context of reactive scattering for the given PES. All the simulations are carried out with no initial rotational motion for the target molecule. The vibrational ground state is referred with the notation ν = 0 which indicates that all vibrational quantum numbers are zero. The vibrational mode corresponding to the symmetric C-H excited stretch in CHD 3 has the vibrational quantum number ν 1 and a harmonic frequency of Ω 1 = 3147 cm −1 .
NE-RPMD simulation were run with n = 3 βω max where βω max refers to the smallest integer i such that i ≥ βω max and ω max is the highest harmonic frequency of the polyatomic reactant. This is the number of beads beyond which we observe little change for the ICS. The time step employed for the QCT and NE-RPMD simulations is 0.02 fs. The initial distances between the reactants' center of mass are 12 a 0 for H + CH 4 /CHD 3 and 14 a 0 for F + CHD 3 . The maximal impact parameters found for the NE-RPMD simulations are in most cases higher than in QCT. For the H + CH 4 /CHD 3 ground state reactions b max = 6 a 0 is employed, while we have b max = 8 a 0 in the presence of the excited C-H stretch in CHD 3 . For the reaction F + CHD 3 , b max is set to 10 a 0 . For each collision energy and reactant vibrational state, 20,000 trajectories were used to compute the ICS. To benchmark our NE-RPMD approach, we compare the NE-RPMD ICS results with our QCT results and with RDQD results taken from the literature. For the reaction H + CH 4 we compare our results with RDQD ICS taken from Refs. 98 and 99 and for the reaction F + CHD 3 from Ref. 34.

B. Stability of the centroid vibrational C-H stretch excitation model
ZPE constraints, which in particular forbid the leakage of the zero-point vibrational energy to other degrees of freedom, are described in RPMD intrinsically via the internal ring polymer structure. Conversely, a potential issue when performing initial excitation of a vibrational mode on a ring polymer polyatomic reactant via "kicking" certain coordinates is the possible leakage of the excitation energy to other modes on the time period it takes for the reactants to collide This leakage occurs because we are not involving the ring polymer spring terms while performing the excitation, rather only the centroids [see Eq. (6)]. It is therefore crucial to monitor how the energy in the excited vibrational mode evolves as a function of time which requires a mapping of phase space coordinates to vibrational quantum states. For QCT, the most straightforward assignment of quantum states is the determination of the corresponding harmonic actions numbers. This is done by employing a semi-classical quantization of the vibrational modes in the harmonic regime as described in Ref. 116. Here, we employ an analogous procedure for the ring polymers to assess the stability of the C-H stretch vibrational excitation.
We perform initial sampling corresponding to the C-H stretch excited state as explained in subsection III A employing the harmonic approximation of the potential. The ring polymer is then propagated employing the full potential. Computing harmonic action numbers to each vibrational mode along the trajectory requires evaluating the displacements of the molecule in terms of vibrational mode coordinates. For this task, it is desirable to align the molecule at each time step of the simulation via translations and rotations as close as possible to an equilibrium reference geometry. Following the procedure described by Eckart 122 , we shift the molecular coordinates to the center-of-mass frame and employ a rotation defined by the rotation matrix C, such that whereq i is the ring polymer centroid position of the ith atom and q eq i is the position of the ith atom in the equilibrium reference geometry of the molecule. 123 The ringpolymer phase space coordinates are then determined as vibrational mode coordinates and momenta as described in Eq. (4) and Eq. (5). The energy for the ith vibrational mode is computed as and a real harmonic action number for each mode is obtained asν In order to quantify the extent to which the excitation energy leaks to other modes we have run 100 trajectories for the ring polymer CHD 3 (ν 1 = 1) molecule, and we have computed the harmonic action number for the C-H stretch modeν 1 at each time step employing the procedure described above. The harmonic action number is shown in Fig. 2 as function of time. For this simulation, we used β = 450 as this value is representative of the β values employed for the collision simulations in this work. We compare these results with the harmonic action numbers obtained from QCT simulations. An initial spike at t = 0 can be seen for both, QCT and NE-RPMD. This spike can be attributed to the anharmonicities in the full potential, which are not considered in the preparation of the initial vibrational state. We observe that for times t ≤ 250 fs after excitation, in which most of the reactions studied here occur, more than 80 percents of the initial energy remains in the initially excited mode and only for longer times the harmonic action number becomes lower than 0.8. For times t < 400 fs, the harmonic action numbers obtained from the QCT simulations tend to decay a bit faster than the ones obtained from NE-RPMD. Overall, the harmonic action numbers obtained from the QCT and the NE-RPMD evolve qualitatively similar. We conclude that the simulations maintain the vibrationally excited state for sufficiently long times. A propagation time of t ≥ 250 fs becomes potentially problematic for reactions at very low collision energies. As we will discuss later, leakage of excitation energy becomes a problem for the studied F+CHD 3 (ν 1 = 1) reaction, where very low collision energies are sufficient to trigger a reaction. In the following we discuss the collisions H + CH 4 and H + CHD 3 , where the target molecules are in their rovibrational ground states. Figure 3 shows the calculated QCT and NE-RPMD ICS as a function of the collision energy for the reactions H + CH 4 → H 2 + CH 3 [ Fig. 3(a)] and H + CHD 3 → H 2 + CD 3 [ Fig. 3(b)]. We compare this data with ICS calculated via RDQD from Ref. 98. These results are based on the reduced-dimensional model introduced by Palma and Clary 27 employed alongside an established wave packet method to study atom-methane reactions. 98 It has been demonstrated that the frozen coordinates in this reduced-dimensional model act largely as spectators in the reaction dynamics. The model has also been validated by demonstrating that ICS for the ground state reaction calculated via QCT are unaffected by the dimensional reduction for a wide range of collision energies and different isotopic variants. 15 The ICS results from these RDQD simulations are therefore considered to be quantitatively accurate. For both reactions we observe that RDQD predict almost the same reaction thresholds at around 0.4 eV, below which the reactivity is zero. This is because both reactions have the same ZPE corrected barrier and tunneling effects provide for both reactions the same contributions. 98 The QCT thresholds for both reactions are predicted at around 0.3 eV, which is substantially lower than in RDQD. The lower threshold for QCT is explained by the presence of ZPE leakage that allows vibrational energy to leak into the reaction coordinates and enables the reaction to happen below the barrier height. 15 On a side note, the absence of tunneling contributions in QCT, which increases the threshold, is overcompensated by the ZPE leakage.
Remarkably, NE-RPMD predicts a reaction threshold approximately around 0.4 eV for H + CH 4 and approximately 0.45 eV for H + CHD 3 in good agreement with RDQD. These results indicate that the ZPE leakage is well mitigated with NE-RPMD. As the collision energy increases, the QCT ICS for the H + CH 4 → H 2 + CH 3 reaction increases more rapidly and also to larger values than the RDQD ICS. A similar observation can be seen for H + CHD 3 → H 2 + CD 3 . Overall, the RDQD cross sections are much smaller compared to QCT, only reaching 0.3 a 2 0 at most compared to 1.3 a 2 0 for QCT for H + CH 4 at 1 eV, while for H + CHD 3 at 1 eV QCT predicts an ICS of 0.16 a 2 0 and 0.075 a 2 0 for RDQD. This suggests that NQEs are relevant for the whole range of collision energies considered above the reaction threshold. More specifically, we attribute this discrepancy to violation of ZPE constraints in QCT that leads to an overestimation of the reactivity. The ICS calculated from NE-RPMD are much closer to those calculated from RDQD. For collision energies above the threshold, NE-RPMD predicts somewhat larger ICS for H + CH 4 and marginally lower ICS for H + CHD 3 compared to RDQD.
We conclude that NE-RPMD covers the major part of the ZPE effects relevant for the two considered reactions. These encouraging results are reminiscent of the previous successes of the method for triatomic systems. We now turn to the ICS for the reaction H + CHD 3 → H 2 + CD 3 in the presence of the C-H excited stretch. Figure 4 shows a comparison of ICS as a function of collision energy calculated from NE-RPMD, QCT, and RDQD. As for the vibrational ground-state, the RDQD ICS results are taken from Ref. 98. Accordingly, they involve the same reduced dimensional model as for the ground state reactions.
From the comparison with the ground state ICS (Fig. 3  (bottom)), one can see that all simulations consistently indicate that in the presence of the C-H stretch the reactivity is increased substantially which is consistent with earlier simulations. 24,[96][97][98] The QCT simulations indicate a reaction threshold of E col 0.2 eV. In contrast to the ground state reactions discussed in Fig. 3, where the reaction threshold for QCT was lower than for RDQD, the reaction threshold for QCT is here somewhat larger than the reaction threshold for RDQD of approximately 0.12 eV. This suggests that the excitation energy deposited into the C-H vibrational stretch mode cannot transfer as efficiently to the reaction coordinate in QCT as compared to RDQD. Also, QCT simulations do not account for any tunneling contributions to the reactivity that may be responsible for the lowered reaction threshold. Interestingly, NE-RPMD predicts a reaction threshold of approximately 0.08 eV which is slightly lower than the one for RDQD. For collision energies up to 0.6 eV we observe that NE-RPMD predicts ICS in good agreement with RDQD while QCT underestimate slightly the reactivity. For collision energies from 0.6 eV to 1.4 eV, both QCT and NE-RPMD ICS remain slightly above the RDQD ICS and NE-RPMD results tend to merge with the QCT predictions at these higher collision energies.
In summary, we find that the NE-RPMD ICS follow the same qualitative trend as RDQD for low collision energies and becomes close to QCT ICS at high collision energies. Compared to the ground state reaction, all simulation methods consistently indicate a considerable reactivity enhancement (around 10 fold) caused by the C-H stretch (ν 1 = 1).
E. F+CHD3(ν1 = 0, 1) Here we compare the ICS for the reactions F + CHD 3 (ν 1 = 0, 1) → FH + CD 3 computed with QCT, RDQD, and NE-RPMD. The comparison between RDQD and QCT and between RDQD and NE-RPMD for this reaction will be somewhat semi-quantitative as an additional approximation 26 alongside with the aforementioned Palma-Clary reduced-dimensionality model (see section IV C) is employed to compute the RDQD ICS, which are taken from Ref. 34 There has been many discussion about the intriguing discrepancies between experimental results and the QCT and RDQD predictions regarding the influence of the C-H stretch excitation on the dynamics of the F + CHD 3 reaction 34,46,103,109 . It is not yet clear what is the cause of these discrepancies and whether the experiments or the simulations are lacking accuracy. One of the possible reasons for disagreement is regarding the QCT simulations the imperfect treatment of quantum effects such as ZPE and tunneling and regarding the RDQD model potential effects of dimensionality reduction. In this context, the presented results may thus provide additional information for investigating this discrepancy, since NE-RPMD describe NQEs to some extend, but on the other hand does not rely on dimensional reduction (as RDQD). Another possible cause of deviation are the limitations of the potential energy surface. In this work, we focus on assessing the performance of our NE-RPMD approach for computing ICS and describing the quantum effects on the PWEM-SO PES.
In contrast to the reactions considered in Fig. 3 and Fig. 4, ZPE constraints are expected to be less relevant, because of the higher mass of the attacker atom. Accordingly, previous works based on QCT trajectories showed that reactivity trends are unaffected by ad-hoc treatments regarding the ZPE of the products. 106 The calculated ICS are shown in Fig. 5 as a function of collision energy and compared with the earlier published RDQD results 34 . For the ground state reaction [ Fig. 5(a)], the QCT reaction threshold is at around E col = 0.03 eV, which corresponds to the classical barrier for the PWEM-SO PES. 109 The RDQD reaction threshold is considerably lower (< 0.010 eV); a complete decline of the ICS is not visible within the range of collision energies considered here. QCT predicts much lower ICS than RDQD for the entire range of collision energies. This has been attributed to the relevance of tunneling contributions that are absent in QCT. 34 NE-RPMD predicts ICS that are somewhat larger than the ICS provided by QCT. They are still significantly smaller, but slightly closer to the values provided by RDQD.
These results indicate that NE-RPMD only captures a rather limited part of the tunneling contributions involved in the reaction, which is in accordance with a similar observation for F+H 2 (ν = 0). 77 A possible cause is the relatively low employed β that severely limits the spatial extension of the ring polymer of the atom F, which in turn limits the description of the tunneling effects.
For the C-H stretch-excited case the calculated ICS are compared in Fig. 5 (b). Similar to what we have seen before for H + CHD 3 , vibrational excitation of the C-H bond stretch increases the reactivity considerably. For the collision energies considered here, we observe that the QCT ICS are considerably lower than their RDQD counter-parts. NE-RPMD predicts slightly lower ICS as compared to QCT except at higher collision energies (≥ 0.14 eV) where NE-RPMD ICS are a bit higher than QCT. Compared to RDQD, NE-RPMD underestimates the reactivity considerably. The modeling of the excited state reaction via NE-RPMD faces particular difficulties, mainly due to the fact that the reaction is nearly barrierless. As a result, the collision energies of interest involve a long simulation time and thus substantial vibrational energy leakage (up to 40 % near the threshold) of the C-H stretch vibrational energy to the other vibrational modes, which are not as efficient at promoting the reaction.
In addition to the discussed causes for the discrepancies for the ground and vibrational excited state, we note that the reaction involves resonance effects that are not described in NE-RPMD and that can strongly affect the reaction dynamics. 34 Those considerations constitute a limit of our NE-RPMD simulations for very low collision energies and contribute to NE-RPMD underestimating the reactivity.
for the corresponding reaction. Note that the range of collision energies in this figure is considerably larger than in Fig. 5. Displaying R helps to assess the efficiency of the C-H stretch excitation in promoting the reaction. For all collision energies, it is observed that the results from QCT, NE-RPMD, and RDQM are systematically larger than 1, i.e., excitation increases the reactivity (in contrast to the earlier experimental findings 103,105,111 ). The QCT ratios are significantly higher in comparison to its RDQD counterpart for the 0.5 -2.5 eV collision energy range. As mentioned before, this is due to the absence of tunneling contributions in the QCT dynamics for the ground state reaction. 34 NE-RPMD predictions of the enhancement ratio are closer to the RDQD results. Given the large discrepancy of the individual ICS, the good agreement for the enhancement ratio must be attributed to some cancellation of error and that there is some additional reactivity at low collision energies for the ground state due to tunneling contributions that are covered within NE-RPMD.

V. CONCLUSIONS
The current work extends the scope of NE-RPMD to the study of state-resolved reactive scattering dynamics for reactions involving more than 3 atoms. We have compared the ICS results for the reactions H + CH 4 , H + CHD 3 , and F + CHD 3 with and without vibrational C-H stretch excitation in CHD 3 . Good agreement for the H+CH 4 /CHD 3 /CD 4 ground state reactions were found. We show that the NE-RPMD simulations describe the dominant ZPE effects present in the dynamics for a wide range of collision energies. However, we also see that NE-RPMD only performs marginally better than QCT for the reaction F+CHD 3 , where relevant tunneling effects cannot be described with NE-RPMD.
For vibrationally excited states, we find that the NE-RPMD method can also give accurate ICS that are sim-ilar to those obtained from much more involved RDQD simulations. The centroid stretch excitation scheme employed to excite the C-H stretch vibrational mode in CHD 3 was found to be robust and fairly accurate for the time period of most of the simulations (if ≤ 200 fs). When employed to study the C-H excited H + CHD 3 reaction, good agreement with the RDQD ICS is found.
In summary, except for the reactive collision F+CHD 3 , where tunneling effects are dominant, NE-RPMD simulations show good agreement with RDQD simulations, especially for collision energies around the threshold. For the discrepancies that arise at larger collision energy, one has to keep in mind that the dimensional reduction on which RDQD results rely may not be strictly valid anymore. Compared to QCT, we can conclude that the here studied NE-RPMD method is able to overcome key limitations of the previously applied QCT methodologies with little extra computational efforts. NE-RPMD can be relatively easily built in into existing MD codes and can be applied to considerably larger systems thanks to its scalability. However, the current NE-RPMD method employs an inconsistent and unclear choice for the β parameter. As such, our NE-RPMD approach as it stands and the present choice of β can be seen as a possible starting point for further refinements of the methods toward broader and consistent applications of NE-RPMD for state-selective reactive scattering simulations.