Rate coefficients for the O + H₂ and O + D₂ reactions: how well ring polymer molecular dynamics accounts for tunelling

Abstract

We present here extensive calculations of the O(³P) + H₂ and O(³P) + D₂ reaction dynamics spanning the temperature range from 200 K to 2500 K. The calculations have been carried out using fully converged time-independent quantum mechanics (TI QM), quasiclassical trajectories (QCT) and ring polymer molecular dynamics (RPMD) on the two lowest lying adiabatic potential energy surfaces (PESs), 1³A′ and 1³A′′, calculated by Zanchet et al. [J. Chem. Phys., 2019, 151, 094307]. TI QM rate coefficients were determined using the cumulative reaction probability formalism on each PES including all of the total angular momenta and the Coriolis coupling and can be considered to be essentially exact within the Born–Oppenheimer approximation. The agreement between the rate coefficients calculated by using QM and RPMD is excellent for the reaction with D₂ in almost the whole temperature range. For the reaction with H₂, although the agreement is very good above 500 K, the deviations are significant at lower temperatures. In contrast, the QCT calculations largely underestimate the rate coefficients for the two isotopic variants due to their inability to account for tunelling. The differences found in the disagreements between RPMD and QM rate coefficients for the reactions for both the isotopologues are indicative of the ability of the RPMD method to accurately describe systems where tunelling plays a relevant role. Considering that both reactions are dominated by tunelling below 500 K, the present results show that RPMD is a very powerful tool for determining rate coefficients. The present QM rate coefficients calculated on adiabatic PESs slightly underestimate the best global fits of the experimental measurements, which we attribute to the intersystem crossing with the singlet 1¹A′ PES.

---

1 Introduction

The gas-phase reaction between atomic oxygen in its electronic ground state and molecular hydrogen leading to a hydroxyl radical and a hydrogen atom is of paramount importance in the chain reaction of the H₂ + O₂ combustion. It is also one of the main sources of OH in the mesosphere in collisions with vibrationally excited H₂. Similarly, it seems to be important in the production of OH in photon-dominated regions where H₂ can be excited to sufficiently high rovibrational states.^1–5 For all these reasons, it is one of the benchmark systems in experimental and theoretical kinetics, and their thermal rate coefficients have been measured in a wide range of temperatures using different techniques such as flow methods, shock tubes or flash photolysis^6–20 covering a temperature range from 300 K to 5000 K. For a relatively recent compilation of rate coefficients, see ref. 21 and 22. The kinetics of the O(³P) + D₂ reaction has been also the subject of different measurements,^{8,15,19,23,24} making it possible to obtain the kinetic isotope effect (KIE), and to evaluate the importance of tunelling for the present reaction.

The reaction is slightly endothermic, , but the reactivity is limited by a large electronic barrier of 0.59 eV (≈0.53 eV for O + H₂ when the zero point energies of reactants and transition state are included). Experimental crossed-beam experiments were hampered by the barrier of the reaction,^25–27 requiring the use of the hyperthermal atomic-oxygen beam source but, nevertheless, it was possible to measure integral cross sections,²⁸ differential cross sections,²⁹ and even Λ-doublet propensities of OD(²Π) produced in the O + D₂ reaction.^30,31 In most crossed-beam experiments, D₂ is preferred over H₂ since for the same velocity of atomic oxygen in the laboratory-frame it leads to higher collision energies, making it somewhat easier to overcome the experimental difficulties arising from the high electronic energy barrier. Motivated by these and previous experiments, a myriad of dynamical calculations have been carried out intended to simulate the experimental findings and calculate the reaction rate coefficients.^{25,27–29,32–37}

The present work is aimed to compare the experimental rate coefficients for O(³P) + H₂ and O(³P) + D₂ with those calculated using rigorous quantum mechanical (QM), quasi-classical (QCT) and ring polymer molecular dynamics (RPMD) calculations, which have been carried out on the two adiabatic 1³A′ and 1³A′′ PESs, and the respective results are combined to get the total reaction rate covering a range of temperatures between 200 K and 2500 K.

Over the last decade, RPMD has proven to be a very efficient method for the determination of thermal rate coefficients for a wide variety of chemical reactions with remarkable accuracy.^38–45 In spite of the use of classical mechanics, given the isomorphism between quantum and classical statistical mechanics of harmonic ring polymers, RPMD has been shown to account for the zero point energy (ZPE) problem almost exactly and with tunnelling, if not completely, then at least to a considerable degree.^46–48 The present calculations of O + H₂ and O + D₂ on the two triplet PESs provide an excellent opportunity to test the performance of RPMD in the treatment of tunelling by comparison with fully converged, essentially exact, QM calculations on the same PESs. As it happens, RPMD almost exactly reproduces the QM results for the reaction with D₂. For the reaction with H₂ the agreement is also very good at temperatures above 500 K, but at lower temperatures, the agreement is not as good as that for the O + D₂ reaction.

The article is laid out as follows: Section 2 presents a brief review of the theoretical methods used in this work; Section 3 covers the computational details of the QM, RPMD and QCT calculations; the results and their discussion are shown in Section 4. The conclusions are presented in Section 5.

2 Theoretical background

2.1 Potential energy and minimum energy paths

The nine states of O(³P_Ja), J_a = 2, 1, 0, (where J_a is the atomic total electronic momentum) upon interaction with H₂/D₂ in their ground electronic state give rise to three triplet PESs, 1³A′′, 1³A′ and 2³A′′. Of the five degenerate sates of the ³P₂ manifold, three correlate with the 1³A′′ PES, and the remaining two correlate with the 1³A′ PES, which also correlates with one of the three degenerate states of the ³P₁ manifold. The other two spin–orbit states of ³P₁ as well as the ³P₀ state correlate with 2³A′′. Whilst both the 1³A′′ and 1³A′ PESs also correlate with the products in their ground state, OH(²Π) + H(²S_g), 2³A′′ is highly repulsive and correlates with the products in an excited state. Therefore, the reaction can take place adiabatically on the 1³A′′ and 1³A′ PESs (henceforth ³A′′ and ³A′), whose saddle points correspond to a collinear arrangement. Moreover, except for large O–H₂ distances, the two PESs are degenerate in the linear configuration and exhibit the same barrier. However, as the system draws away from the linear configuration, the degeneracy is broken and the ³A′ PES shows a steeper bending potential leading to a narrower cone of acceptance. As a result of the steeper bending potential, the bending frequency is 70% higher on the ³A′ PES,²⁷ leading to a vibrationally adiabatic barrier (including the ZPE) that is 40 meV higher than on the ³A′′ PES for the O + H₂ reaction.

As given in previous works,^{27,33,37,49,50} calculations have been carried out separately on the two adiabatic PESs, of symmetries ³A′′ and ³A′ calculated by Zanchet et al.,²⁷ which accurately reproduces the degeneracy between both PESs for collinear geometries, including the saddle point. The PESs were fitted using a many-body expansion based on ab initio energies calculated using the internally contracted multireference configuration interaction method including Davidson correction (icMCRCI + Q) based on a wave function calculated at a state-average complete active space (SA-CASSCF) level of theory including in the active site all valence electrons of the three atoms. Calculations on these PESs were able to reproduce the experimental Λ-doublet propensities of OD produced by collisions between atomic oxygen and a D₂ molecule.³⁷ These global analytical PESs overestimate the ab initio barrier height by only 6 meV for both ³A′′ and ³A′ states. In contrast, the GLDP fit of the RWKW PES by Rogers et al.²⁵ underestimates the ab initio barrier by 26 and 18 meV, respectively, and hence they are not degenerate in the collinear approach.

Fig. 1 shows the comparison of the vibrationally adiabatic minimum energy paths (MEPs) for the two isotopic variants on the ³A′ and ³A′′ PESs, related to the respective zero point energies of H₂ and D₂. The MEPs are plotted as a function of the mass-scaled reaction coordinate, s, which measures the progress of the reaction, and is given by

with

Q₂ = μ_BC^1/2R_BC, (3)

where m_A, m_B, and m_C are the masses of A, B, and C, respectively. μ_BC and μ_A–BC are the reduced masses of BC and A–BC, and R_AB and R_BC are the internuclear distances A–B and B–C. In the present case, A refers to O(³P) and BC to H₂ or D₂. Q^SP₁ and Q^SP₂ are the values of Q₁ and Q₂ at the saddle point. Eqn (1) only permits evaluation the absolute value of s, such that the negative/positive values of s are chosen to correspond to the reactants/products valleys.^51,52 Note that s = 0 corresponds to the geometry of the saddle point of the electronic potential without ZPE (see Fig. S1, ESI†). When the contribution of ZPE is added to the potential, the saddle point is found at slightly negative values of s. Fig. S1 (ESI†) depicts the electronic MEPs (dashed lines) for the two reactions, which are the same on the ³A′ and ³A′′ PESs. In this figure the zero energy is that of the minimum of the asymptotic potential. Because the vibrational ZPE of D₂ is 0.192 eV whilst that of H₂ is 0.270 eV, the adiabatic barriers appear to be higher for the O + H₂ reaction. When the potentials are referred to their respective asymptotic vibrational energies, the resulting MEPs are those of Fig. 1. As shown in this figure, the vibrationally adiabatic minimum energy paths for the O + D₂ reaction are considerably broader and have slightly higher barriers than for the reaction with H₂. The respective heights of the vibrationally adiabatic barriers for O + H₂ calculated on the ³A′ and ³A′′ PESs are 0.549 eV and 0.507 eV, respectively. The difference between the barrier heights of the ³A′ (0.558 eV) and the ³A′′ (0.531 eV) are smaller for the O + D₂ due to the larger reduced mass associated with the bending vibrational mode.

Fig. 1 Vibrationally adiabatic minimum energy paths for O(³P) + H₂ (blue) and O(³P) + D₂ (red) reactions as a function of the mass-scaled reaction coordinate (see text for definition) on the ³A′ (upper panel) and ³A′′ (bottom panel) PESs calculated by Zanchet et al.²⁷ In both cases, the minimum energy paths for the two isotopes are referred to as the respective ZPE of H₂ and D₂.

Given the differences in the MEPs of the two isotopic variants, especially in their widths, the expected tunelling contribution to the reaction will be more important for O + H₂. This is underpinned by the crossover temperatures, T_c, for the two isotopic variants, which is given by

where iν_TS is the imaginary frequency at the transition state. Roughly speaking, T_c indicates the maximum temperature for which the reaction is governed by tunelling; i.e., when the temperature is lower than T_c, the system enters the deep tunelling regime and the RPMD k(T) is expected to deviate from the accurate QM results. A more rigorous definition can be found in ref. 44, 48 and 53. The crossover temperatures (imaginary frequencies of the saddle point) for the O + H₂ reaction are 417.7 K (1824 cm⁻¹) and 414.7 K (1810.85 cm⁻¹) on the analytical ³A′ and ³A′′ PESs, respectively. For the O + D₂ reaction, the respective values are 302.5 K (1321.07 cm⁻¹) on the ³A′ PES, and 312.09 K (1312.09 cm⁻¹) on the ³A′′ PES. As can be seen, there is a difference of 100 K in the crossover temperatures of the two isotopic variants. Therefore, tunelling is paramount to the O + H₂ reaction below 400 K whilst for the reaction with D₂ one would expect no deviations from RPMD above 300 K.

In this work, we will assume that the two triplet states ³A′ and ³A′′ are uncoupled. Our calculations do not include the intersystem crossing (ISC) between the singlet 1¹A′ and the triplet 1³A′′ and 1³A′ PESs. Schatz and coworkers were the first to study the influence of ICS on the O + H₂ reaction.^36,54 They determine the spin–orbit coupling matrix using a four-state model. It was found that the singlet state crosses the two triplet states near the barrier (see Fig. 1 of ref. 36). They performed trajectory surface hopping calculations and found that the effect of ICS is relatively minor and is only noticeable at collision energies above 1 eV. So even if some collisions sample the singlet PES, the spin–orbit coupling is expected to have a relatively small effect on the dynamics of the system. Similar conclusions have been drawn from a QM non-adiabatic study for the O(³P) + D₂ reaction, where, interestingly, the influence of the ISC was found to be larger than for the H₂ reaction.⁵⁵ The Renner–Teller coupling between the ³A′ and ³A′′ PESs, neither is considered. It is expected to have some effect,⁵⁶ which would be more pronounced at high orbital angular momentum (impact parameter) values. Although non-adiabatic effects are expected to be small, given the general good agreement with adiabatic results,^28,32,33,49 the comparison of the current accurate adiabatic calculations on the two separate recent PESs with the experimental rate coefficients may shed light on their importance.

2.2 QM and QCT cumulative reaction probabilities

Calculation of QM and QCT rate coefficients have been carried out by means of the cumulative reaction probabilities (CRP), which will be denoted by C_r(E).^57–62 The QM CRP at a given total angular momentum, J, and total energy, E_tot, measured from the asymptotic minimum of the reactant's potential, is given bywhere v, j and v′, j′ are the vibrational and rotational quantum numbers for reagents and products, and Ω and Ω′ are the respective projections of J onto the incoming and outgoing asymptotic directions. is the S matrix element connecting reactant's and product's states, and P^J_v,j,Ω(E_tot) is the reaction probability for initial v, j, Ω summed over all final states. The summations run over all possible states of reactants and products compatible with E_tot. The total CRP is given by

In the case of homonuclear molecules, the diatomic parity, p, has to be taken into account. The parity dependent CRP can be written as

where the summation runs over even (p = e) or odd (p = o) rotational quantum numbers j_p.

The reaction rate coefficients can be written in terms of the total CRPs

where β = 1/k_BT. The translational partition function is given byand Q_B2(T) (B₂ = H₂/D₂) is the coupled nuclear-rovibrational partition function given bywhere g^(o)_n is the degeneracy of the nuclear spin for odd rotational states, (I + 1)(2I + 1) and I(2I + 1), and g^(e)_n for even rotational states, I(2I + 1) and (I + 1)(2I + 1) for H₂ and D₂, respectively.

Of particular interest is the thermal reaction probability, which can be defined as

The thermal-CRP represents the contribution of the Boltzmann weighted total CRP to the rate coefficient in an interval of energies between E_tot and E_tot + dE_tot. Its integration over the total energy range yields k(T).

The above equations assume that there is ortho and para equilibrium at all temperatures. However, in most experiments, since the ortho ↔ para conversion is very slow, the fraction of ortho/para usually is the same as that at room temperature. At 300 K and above, the respective weights are w^(o/e) ≈ g^(o/e)_n/(g^(o)_n + g^(e)_n). The rate coefficients calculated under these circumstances can be written as⁶³

where Q^(p)_vj(T) are the rovibrational partition functions for even and odd rotational states. The corresponding thermal reaction probability is given by

Except at temperatures below 200 K, the rate coefficients obtained using eqn (8) and (12) are the same within the accuracy of the various dynamical methods. The details of the methodology used to calculate the QCT C^J_r(E) are given elsewhere.^62,64,65

In eqn (8) and (12) the calculation of rate coefficients only considers the contributions from one of the concurrent PESs. To simulate the experimental conditions, it is necessary to average the k(T) obtained for the ³A′ and ³A′′ PESs. When the correlation between the oxygen atomic states and the two triplet states is taken into account, as mentioned above, the ³A′ PES correlates with two of the five components of O(³P₂), and one component of O(³P₁), whilst the ³A′′ PES correlates with the remaining three components of O(³P₂).⁴ Under these conditions, k(T) is given by

where D is the electronic partition function of the ³P_Ja states:

D = 5 + 3 exp(−βΔE₁) + exp(−βΔE₀), (15)

and ΔE₁ and ΔE₀ are the energy differences between the ground state, O(³P₂), and O(³P₁) and O(³P₀), respectively. In this case, ΔE₁ = 227.708 K and ΔE₀ = 326.569 K.

At sufficiently high temperatures k_BT ≫ ΔE₁ and ΔE₀, and then

2.3 Ring polymer molecular dynamics

The RPMD approach exploits the isomorphism between a statistical quantum system and a fictitious classical ring polymer consisting of harmonically connected beads, allowing the real-time evolution of the quantum system to be approximated by classical trajectories. The RPMD exhibits some important features: (a) it becomes exact in the high temperature limit, where the ring polymer collapses to a single bead; (b) it is independent of the choice of the transition state dividing surface; (c) it preserves the zero-point energy (ZPE) along the reaction path; (d) it has a well-defined short-time limit, which provides an upper limit of the RPMD rate coefficient. In addition to these features, RPMD provides the exact solution for tunnelling through a parabolic barrier. Thereby, it can be expected that more realistic barriers will account for most of the tunelling contribution to reaction, even in the deep tunnelling regime.^38,39,44

A detailed description of the RPMD rate theory can be found in ref. 40, 41, and 44. The technical aspects of the computational procedure are well documented in the manual of the general RPMDrate code.⁴¹

The ring polymer Hamiltonian of a system consisting of N atoms with fictitious ring polymers of n_b beads can be written in atomic Cartesian coordinates as

where q_i,j and p_i,j are the position and momentum of the jth bead of the ith atom of the system, and q_i,j+nb = q_i,j for the polymer to be closed. The angular frequency of the harmonic springs is ω_n = [βħ/n_b]⁻¹ and β = 1/k_BT, where T is the temperature of the system.

The method begins by introducing two dividing surfaces: one located in the asymptotic reactant valley, s₀(q̄) = 0, and the other located in the transition state region, s₁(q̄) = 0.^41,44 The reaction coordinate ξ is taken to be an interpolating function that connects these dividing surfaces:

such that ξ → 0 as s₀ → 0 and ξ → 1 as s₁ → 0. These surfaces are defined in terms of the centroids of the atoms, which can be calculated by simply averaging the positions of all corresponding beads,

The correlation function formalism used in the computational procedure for the RPMD rate coefficient calculation is based on the t → +∞ limit of the ring polymer flux-side correlation function C_fs.⁴¹ The rate coefficient is then expressed using the Bennett–Chandler factorization^66,67 as

k_RPMD(T) = k_QTST(T;s₁)κ(s₁) (20)

where k_QTST(s₁) is the centroid-density quantum TST rate coefficient and is evaluated aswhere μ_R is the reduced mass of the reactants and R_∞ is the distance between the center of mass of the reactant molecules, such that the interaction potential becomes negligible. The factor p(s₁,s₀) is the quotient of the correlations functions for two different dividing surfaces which can also be expressed in terms of the centroid potential of mean force (PMF), or free energy along ξ, W(ξ):

The second factor κ(s₁) in eqn (20) is the long time limit of the time-dependent ring polymer transmission coefficient, given by

It represents a dynamical correction to k_QTST(T) that accounts for recrossing of the transition-state dividing surface at t → t_p, where t_p is the time in which κ remains invariant.

In practice κ and p are calculated at the maximum free energy value, W(ξ^‡), along the reaction coordinate ξ, which is a function used to connect the two dividing surfaces, and varies from ξ → 0 as s₀ → 0 to ξ → 1 as s₁ → 0.^40,44

3 Computational details

As mentioned above, all calculations, QCT, RPMD and QM, were performed on the triplet PESs of Zanchet et al.²⁷ Extensive TI QM scattering calculations were carried out using the coupled-channel hyperspherical coordinate formalism as implemented in the ABC code.⁶⁸ Calculations were performed for the two triatomic and diatomic parities and J ∈ [0,62] for the O + H₂ reaction and J ∈ [0,80] for O + D₂, and 60 total energies in the 0.37–2.50 eV range. The propagation was performed using 300 sectors with a maximum hyperradius of 15a₀, including on the basis of all the rovibrational states whose internal energy is below 3.25 eV, and a maximum value of the helicity quantum number of 25.

RPMD calculations have been performed using the RPMDrate code⁴¹ and the simulation parameters are similar to those used in previous works.^44,45,48 For the O + H₂ reaction, calculations were carried out using n_b = 32 in the temperature range 600–2500 K. As shown in Table S1 (ESI†), the n_b required for convergence increase at low temperatures (128 at 200 K and 64 from 250–500 K). For the O + D₂ reaction, convergence was achieved with n_b = 64 in the whole temperature range. Purely classical calculations (not QCT) correspond to setting the number of beads n_b = 1.^39,44

To calculate the PMFs, W(ξ), the path was divided in 115 windows of 0.01 width with the reaction coordinate ξ defined between −0.05 and 1.10. In each of the windows, 100 RPMD trajectories were run, restraining the value of ξ_i to the center of the window adding a harmonic potential to the Hamilton function of the system. Each trajectory was thermalized for 20 ps after which it was run for 100 ps, using a time step of 0.1 fs. Based on the PMF profiles, the RPMD rate coefficients were obtained combining the quantum transition state theory (QTST) with the RPMD recrossing factor, κ (see eqn (20)). For the latter calculations, 10⁵ child trajectories were run for 0.1 ps starting from an initial parent trajectory of 20 ps.

QCT values were computed following the cumulative reaction formalism described in ref. 62, 64 and 65. First, we calculated the cumulative reaction probabilities as a function of the total energy, E_tot, for specific values of J, with all the internal states of H₂ and D₂ microcanonically sampled and calculated using the asymptotic diatomic potential energy of the PESs. Batches of 10⁶ trajectories were calculated for each J at a variable total energy, which was sampled uniformly in the 0.6–4.0 eV range. In addition, the total CRPs, summed over all J, were calculated as a single batch (10⁸ trajectories) for each PES and isotopic variant by randomly sampling the total energy and J values proportionally to the number of projections 2J + 1. Using the data thus obtained, the rate coefficients can be calculated directly without the need to first determine the CRP and avoiding numerical integration. The details of the method can be found in ref. 69. The trajectories were run from an initial atom–diatom distance of 10 Å, and an integration step of 5 × 10⁻² fs, which ensures an energy conservation better than 1 in 10⁵. The total CRP, summed over all J, were calculated for each PES and isotopic variant.

4 Results and discussion

The comparison of QM and QCT CRPs are shown in Fig. 2 and 3 for the O + H₂ and O + D₂ reactions for J = 0 and J = 20 on the ³A′ and ³A′′ PESs, respectively. Analogous results for J = 10, 40 and 50 are shown in the ESI† (Fig. S2 for O + H₂ and Fig. S3 for O + D₂). At energies above the classical threshold, the agreement between QM and QCT is very good for all the represented C^J_r(E_tot). However, the QCT CRP drops off suddenly below the classical threshold (see inset) whilst the QM curves die out very slowly, although with small values. The classical total energy thresholds are slightly higher for the reaction with H₂ than for O + D₂. For J = 0 the respective classical thresholds are below the vibrationally adiabatic barrier measured from the asymptotic reactant minimum (see Fig. S1, ESI†), which for O + H₂ are 0.82 eV and 0.78 eV, and 0.75 eV and 0.72 eV for O + D₂, for the ³A′ and ³A′′ PESs, respectively. This shows that the ZPE of the transition state is not conserved in the QCT calculations. However, the respective classical thresholds (≈0.70 eV and 0.67 eV for O + H₂ and O + D₂, respectively) are larger than the electronic barrier (see Fig. S1, ESI†), indicating that not all the reagent ZPE is employed to overcome the barrier. This effect was also observed in the isotopic variants of the H + H₂ exchange reaction, and in particular in Mu + H₂.^48,70,71

Fig. 2 Comparison of QM (black solid line with open circles) and QCT cumulative reaction probabilities for the indicated total angular momenta, J, as a function of the total energy for the O(³P) + H₂ (blue line, left panels) and O(³P) + D₂ (red line, right panels) reactions on the ³A′ PES at the indicated values of J. In each case, the inset depicts the low energy region in the logarithmic scale. The sudden drop of the QCT calculations corresponds to the classical total energy threshold for the reaction.

Fig. 3 Same as Fig. 2 for the ³A′′ PES.

With increasing J, due to the centrifugal barrier, the classical threshold increases rapidly (Fig. S2 and S3, ESI† for J = 40 and 50). Above the classical threshold, the agreement between QCT and QM CRPs is excellent in the whole range of J. The disagreement in the low energy regime, below the classical thresholds, leads to significantly different rate coefficients at low temperatures, where QCT should underestimate the QM rate coefficients. According to the behaviour of the C^J_r(E_tot), it can be expected that at higher temperatures, where the contributions from energies below the classical threshold become less relevant, QM and QCT results should converge. Regardless of J, the PES and the isotopic variant, the CRP increases monotonically with the energy, and the density of reactive states⁷² only give rise to some oscillations, which could be attributed to quantized transition states for J = 0.

The computational burden of QM calculations of thermal rate coefficients increases considerably with T. This is because calculations at higher T involve scattering calculations at higher E_tot, which in turn require higher values of J and Ω, and also to include more rovibrational states in the basis. The excellent agreement between QM and QCT C^J_r(E_tot) at high E_tot supports the idea that QM calculations at low E_tot can be combined with QCT calculations at higher E_tot, and still yield accurate rate coefficients in the whole range of energies almost indistinguishable from the rate coefficient obtained using only QM calculations.

The thermal-CRP at 500 K and 1000 K are shown in Fig. 4 on the ³A′ and ³A′′ PESs for the two isotopic variants. As commented on previously, the thermal-CRP, C_r(E_tot,T), represents the evolution of the reactivity as a function of the total energy for a given T, and facilitates the interpretation of the k(T), which is the integral of C_r(E_tot,T). Below 500 K, most of the reactivity is associated with energies around or smaller the classical threshold, leading to a large difference between QCT and QM C_r(E_tot,T). At these energies, differences of more than a factor of 2–3 are observed between the heights of the QM and QCT thermal-CRP for O + H₂ on both PESs, and somewhat less for O + D₂. With increasing temperature, the thermal-CRP shifts towards higher energies, where QCT and QM C_r(E_tot,T) are more similar above 1.25 eV (QM and QCT thermal-CRP at 1500, 2000 and 2500 K are shown in Fig. S4 of the ESI† on the two PESs and for both isotopic variants). Overall, the low energy tail is very similar for the two approaches, indicating that the contribution from energies near the classical threshold is less relevant (see Fig. S4, ESI†). Also with increasing T, the high energy tail of the distribution extends to higher energies. In fact, whilst at 500 K only energies up to 1.1 eV are required to converge the k(T), at 2000 K energies up to 2.5 eV have to be included in the calculations.

Fig. 4 Comparison of QM (solid lines) and QCT (dashed–dotted lines) thermal cumulative reaction probabilities (see text for definition). Upper panels: O(³P) + H₂ reaction on the ³A′ (a) and ³A′′ (b) PESs at 500 K (black solid and black dashed–dotted curves) and 1000 K (blue solid and blue dashed–dotted curves). Bottom panels: O(³P) + D₂ reaction on the ³A′ (c) and ³A′′ (d) PESs at 500 K (red solid and red dashed–dotted) and 1000 K (green solid and green dashed–dotted curves).

For the calculation of rate coefficients at higher temperatures, E_tot values of above 3.0 eV are required (see Fig. S4, ESI†). Covering such a wide range of energies is not a problem for QCT, as the computational cost is only slightly dependent on the total energy. However, the computational effort of QM scattering calculations increases rapidly with E_tot. Due to the good agreement between QCT and QM CRP at high energies, one way to determine QM k(T) s at higher T values is to combine QM and QCT thermal-CRPs. Basically, what is needed is to infer the high energy tail of the QM thermal-CRP using QCT calculations. The procedure is as follows: the QCT thermal-CRP is scaled to match the QM thermal-CRP in the highest energy range where the latter has been calculated. The scaled QCT is used to extrapolate the QM thermal-CRP at energies for which there are no QM calculations. In the present work, the highest QM total energy calculated is 2.5 eV, whereas in QCT calculations, energies up to 4.0 eV were considered. To estimate the error resulting from this approximation, we calculated the k(T) at 2000 K using pure QM calculations and the QM-QCT combination. The discrepancy between the two calculations was found to be approximately 2%.

The potential mean force (PMF) profiles and transmission coefficients from the RPMD on the ³A′ and ³A′′ PESs are shown in Fig. S5 and S6 (ESI†) for the O + H₂ and O + D₂ reactions, respectively. The PMF profiles, W(ξ), for 300, 600, 1000 and 2000 K are plotted in the respective top panels. The maximum, which appear at a value of the reaction coordinate very close to ξ = 1 is slightly higher on the ³A′. Fig. S7 (ESI†) shows the comparison of the potential mean force profiles at T = 300 K and T = 2000 K obtained in the calculations on the two PESs and the two isotopic variants. At 300 K the free energy barriers (maximum of the PMF) follow a sequence dominated by the isotopic reaction (H₂(A′′) < H₂(A′) < D₂(A′′) < D₂(A′)), whereas at 2000 K the free energy barriers follow the sequence dictated by the PES (H₂(A′′) < D₂(A′′) < H₂(A′) < D₂(A′)).

The long limit of the time-dependent ring polymer transmission coefficients, κ(t), reach their maximum value at an early time, ≈30 fs. At high temperatures they converge to ≈0.78 for O + H₂ and to 0.80 for O + D₂. The results with n_b = 1 are equivalent to the purely classical ones, in which the initial quantization (ZPE) is not considered. This is the key distinction with the QCT rate coefficients, which take into account the quantization of the reactants (the detailed dependence of k^QTST and κ on the number of beads can be found in Tables S1 and S2 of the ESI†).

Thermal rate coefficients for the O + H₂ and O + D₂ reactions on the ³A′ and ³A′′ PESs are shown in Fig. 5, where the results from QCT, RPMD and QM calculations are compared. The RPMD k(T) are plotted for n_b = 1 and for the number of beads necessary for convergence. In all cases, the differences between the QCT k(T) with RPMD and QM indicate the importance of tunelling and, to a lesser extent, of the ZPE. The agreement between accurate adiabatic QM (with no approximations) and RPMD rate coefficients is remarkably good for O + D₂ on both PESs. At 200 K the QM rate coefficients are about 63% higher, but at T = 300 K the difference reduces to 16%, becoming negligible at higher temperatures. For O + H₂, the agreement is also good, but at T ≲ 400 K deviations from the QM results become apparent. At 200 K there is a factor of 2.2 between the k(T) calculated by the two methods; in any case, the difference can be considered relatively small, since k(T) changes by a factor of 10⁵ in the 200 K and 500 K interval. The temperatures in which the RPMD results start to deviate from the QM k(T) closely correspond to the respective T_c for the two isotopic variants. Nevertheless, RPMD appears to be a very robust approximation to accurate QM results, and certainly at a much lower computational cost, even when a considerable number of beads are required. This is especially the case at the temperatures above 1500 K, where the QM calculations become computationally demanding, and the RPMD requires even fewer beads for convergence.

Fig. 5 Comparison of calculated rate coefficients from QM (solid line), RPMD (solid line with points) and QCT (dashed–dotted line). Upper panels: O(³P) + H₂ reaction on ³A′ (left) and ³A′′ (right) PESs. Lower panels: O(³P) + D₂ reaction on ³A′ (left) and ³A′′ (right) PESs. For O(³P) + H₂ the RPMD k(T) were calculated with n_b = 32 for T > 600 K, n_b = 128 at 200 K and n_b = 64 in the 250–500 K range. For the O(³P) + D₂ reaction, n_b = 64 at all temperatures. Purely classical calculations (with no reactant quantization) are those obtained with n_b = 1. No electronic partition function is included.

Fig. 6 is presented to further demonstrate the relevance of tunelling on the two PESs for the two isotopic variants in order to assess the ability of the RPMD method to deal with it. The comparison between the QM and RPMD k(T) on the ³A′ and ³A′′ PESs is shown in the figure. At T = 200 K, the QM and RPMD KIEs, k^OH₂(T)/k^OD₂(T), are ≈230 and 170 on both PESs. At 300 K, the respective QM and RPMD KIEs are 27 and 18, and at 500 K they drop to 5 in both calculations. The relevant issue is that RMPD accounts for the kinetic isotope effect very accurately on both PESs above 300 K. We will explore the comparison with measurable KIE averaged over the electronic states later in this section.

Fig. 6 QM and RPMD rate coefficients for O(³P) + H₂ and O(³P) + D₂ reactions on ³A′ (left) and ³A′′ (right) PESs. No electronic partition function is included. Note the excellent agreement between the QM and RPMD for the O(³P) + D₂ reaction on both PESs.

With regard to the comparison of the reactivity on the two PESs, it can be observed that the rate coefficients on the ³A′′ are greater for the two isotopic variants. For O + H₂, the ratio (T)/k^A′(T) is ≈3 at 200 K, 2.1 at 500 K and 1.8 at 1000 K. At 2500 K the rate coefficient for the reaction with H₂ is only 50% bigger than for O + D₂. The ratios obtained from the RPMD are very similar. Two factors contribute the higher reactivity observed on the ³A′′ PES: (i) as illustrated in Fig. 1, the adiabatic barrier is slightly smaller and somewhat narrower, so tunelling is likely to be more efficient in this PES; (ii) the bending potential on the ³A′ PES is steeper leading to a narrower cone of acceptance, whose effect is especially noticeable at collision energies above 1.0 eV. The higher k^A′′(T)/k^A′(T) ratio at low temperatures suggests that the former effect is more important. This result agrees with the energy dependence of the ratio between the cross section on the ³A′ and ³A′′ obtained for the O + H₂ (j > 0).^49,50 The k^A′′(T)/k^A′(T) ratio at low T is slightly smaller for O + D₂ reaction, confirming that tunelling is less important for this isotopic variant.

The comparison of the experimental data^{8–10,13–17,20} and the theoretical k(T), are shown in Fig. 7 and 8 for the O + H₂ and O + D₂ reactions, respectively. For a proper comparison, it is necessary to combine the rate coefficients calculated separately on the ³A′ and the ³A′′ PESs, as it has been discussed in Section 2.2. Eqn (14), which takes into account the correlation of the states of ³P₂ and ³P₁ and the ³A′ and ³A′′ PESs, has been used in this work. Although in the high temperature limit, eqn (14) is reduced to eqn (16), even at 1000 K there is a deviation of 7%.

Fig. 7 Comparison of experimental rate coefficients with those obtained with QM (blue solid line), RPMD (green line with solid circles), and QCT (red short-dashed line) for the O(³P) + H₂ reaction between 200 and 2500 K. Experimental results: Westenberg and de Haas,^8,9 Dubinsky and McKenney,¹⁰ Campbell and Handy^11,12 Presser and Gordon,¹³ Sutherland et al.,¹⁴ Marshall and Fontijn,¹⁵ Natarajan and Roth,¹⁶ Shin et al.,¹⁷ Yang,¹⁹ and Ryu et al.²⁰

Fig. 8 Same as Fig. 7 for the O(³P) + D₂ reaction. Experimental results: Westenberg and de Haas,⁸ Presser and Gordon,¹³ Marshall and Fontijn,¹⁵ Michael,²³ Zhu et al.²⁴ and Yang et al.¹⁹

We will first consider the comparison of the present theoretical results with the experimental data. Table 1 shows the theoretical rate coefficients together with the average of the most reliable best-fits to the experimental data for both isotopic variants, which can be found in the ESI.†^{14,15,17,19,21,24} The most common three-parameter empirical expression of k(T) used over a wide range of temperatures is

k(T) = AT^mexp(−E_a/T) (24)

Table 1 Comparison of the O(³P) + H₂ and O(³P) + D₂ total (summing the contributions from ³A′ and ³A′′ PESs, eqn (14)) rate coefficients calculated using the QCT, RPMD, and QM methods. The parentheses denote powers of ten. Units are in cm³ s⁻¹

O + H2					O + D2
T(K)	k(T)^QCT	k(T)^RPMD	k(T)^QM	k(T)^exp	k(T)^QCT	k(T)^RPMD	k(T)^QM	k(T)^exp
200	3.166(−23)	1.398(−20)	3.068(−20)	—	2.093(−24)	7.880(−23)	1.293(−22)	—
300	2.762(−19)	2.661(−18)	4.525(−18)	(9 ± 2)(−18)	4.092(−20)	1.425(−19)	1.665(−19)	(5 ± 2)(−19)
400	2.909(−17)	9.543(−17)	1.218(−16)	(2.4 ± 0.4)(−16)	6.539(−18)	1.236(−17)	1.312(−17)	(3.1 ± 0.2)(−17)
500	5.105(−16)	1.046(−15)	1.205(−15)	(1.9 ± 0.5)(−15)	1.480(−16)	2.196(−16)	2.278(−16)	(5.1 ± 0.1)(−16)
600	3.619(−15)	5.907(−15)	6.420(−15)	(9 ± 2)(−15)	1.246(−15)	1.594(−15)	1.685(−15)	(3.1 ± 0.9)(−15)
700	1.519(−14)	2.187(−14)	2.296(−14)	(2.7 ± 0.8)(−14)	5.930(−15)	7.224(−15)	7.453(−15)	(1.3 ± 0.4)(−14)
800	4.584(−14)	6.062(−14)	6.274(−14)	(7 ± 2)(−14)	1.966(−14)	2.275(−14)	2.361(−14)	(3.9 ± 0.9)(−14)
900	1.106(−13)	1.388(−13)	1.418(−13)	(1.5 ± 0.4)(−13)	5.109(−14)	5.627(−14)	5.955(−14)	(1.0 ± 0.2)(−13)
1000	2.279(−13)	2.613(−13)	2.792(−13)	(2.8 ± 0.7)(−13)	1.117(−13)	1.242(−13)	1.275(−13)	(2.0 ± 0.4)(−13)
1200	7.010(−13)	7.792(−13)	8.108(−13)	(9 ± 1)(−13)	3.749(−13)	4.052(−13)	4.169(−13)	(7.0 ± 0.2)(−13)
1500	2.300(−12)	2.455(−12)	2.544(−12)	(2.9 ± 0.2)(−12)	1.336(−12)	1.400(−12)	1.459(−12)	(2.2 ± 0.2)(−12)
1800	5.320(−12)	5.700(−12)	5.764(−12)	(6.9 ± 0.6)(−12)	3.255(−12)	3.461(−12)	3.528(−12)	(5.4 ± 0.1)(−12)
2000	8.244(−12)	8.765(−12)	8.931(−12)	(1.1 ± 0.2)(−11)	5.164(−12)	5.510(−12)	5.627(−12)	(8.5 ± 0.2)(−12)
2200	1.192(−11)	1.254(−11)	1.290(−11)	(1.4 ± 0.3)(−11)	7.605(−12)	8.172(−12)	8.321(−12)	(1.24 ± 0.02)(−11)
2500	1.881(−11)	2.017(−11)	2.040(−11)	(2.4 ± 0.5)(−11)	1.225(−11)	1.332(−11)	1.351(−11)	(2.01 ± 0.03)(−11)

---

A number of bibliographical expressions that fit the experimental measurements are listed in the ESI,† along with the temperature interval within which the expressions are valid. The most recent global expression for the k(T) of the reaction with H₂ has been provided by Baulch et al.²¹ as the sum of two Arrhenius equations with four parameters, presumably covering the entire range 300–3300 K. It is important to note that the experimental points exhibit an average deviation of 15–20% from the equation presented in ref. 21. The experimental values shown in Table 1 have been obtained by averaging the different expressions at each temperature using the interval where each fit is assumed to be reliable. Fig. 9 presents the measurements of the rate coefficients for the O + H₂ reaction together with the bibliographical best-fits in semi-log plots in order to illustrate the accuracy of the various global expressions. The inset to the figure displays the corresponding linear plots in the temperature range 1500–2500 K, where the QM k(T) is depicted as a red solid line.

Fig. 9 Comparison of individual measurements^8–17,20 and global fits^14–16,19 of rate coefficients for the O + H₂ reaction. The global fit by Baulch et al.²¹ is shown as a black line. The inset shows the different data in the 1500–2500 K temperature range in the linear scale. As can be seen, the recommended global k(T) by Baulch et al. appears as an upper bound of the individual measurements in this temperature interval. The present QM rate coefficients are shown as a red solid line.

The results given by Baulch's global expression are plotted a black solid line, and as can be seen, it represents almost an upper bound of the experimental points in this interval. This is due to the fact that the parameters utilised in this equation also fit the measurements^16,18,73 at temperatures above 2500 K (see ref. 22). However, the double Arrhenius equation is not sufficiently flexible to encompass the entire range of temperatures. As illustrated in Fig. 7 and 9, the present QM results exhibit a slight underestimation of the experimental rate coefficients. The discrepancy between QM and experimental k(T) values is significant up to 500 K, although it never exceeds a factor of 2 even at the lowest temperatures (300 K). The degree of agreement between theory and calculations improves above 700 K.

The agreement between QM and experimental rate coefficients^8,15,23,24 is worse for the reaction with D₂, as can be seen in Fig. 8 and 10. These figures demonstrate that the QM results systematically underestimate the experimental determinations. In any case, according to Table 1, the difference is not greater than a factor of 3 at 300 K, decreasing to 1.5 at temperatures above 1800 K. For this isotopic variant, the deviations between the different experimental data sets is much smaller (as evidenced by the uncertainties associated with the experimental data), and the QM rate coefficients appear to represent the lower limit of the experimental results, as shown in the inset of Fig. 10.

Fig. 10 Comparison of individual measurements^{8,13,15,19,23,24} and global fits^15,16,19,23 of rate coefficients for the O + D₂ reaction. The inset shows a blow-up of the 1500–2500 K range. The red solid line represents the present QM k(T).

With regard to the comparison of the predictions of the rate coefficients by the different theoretical approaches, it can be observed that the agreement between QM and RPMD k(T) above 500 K is very good for the reaction with H₂, as shown in Fig. 7 and Table 1. In the worst case, 300–400 K, the discrepancy between the two approaches reaches up to a factor of 2. For the reaction with D₂, as shown in Fig. 8 and Table 1, the agreement is even better, demonstrating the excellent performance of the RPMD method at temperatures near or above the crossover temperature. As expected, the QCT method predicts very low rate coefficients below 500 K due to the inability of the QCT results to account for tunnelling. At 300 K and 500 K there are differences by a factor of 16, and 2.4, respectively. Only at T > 800 K the respective rate coefficients are comparable. For the O + D₂ reaction, the QCT k(T) are closer to the QM ones. At 300 K, there is a factor of 4 which decreases to 1.53 at 500 K. At T > 800 K, the discrepancy is only of a factor of 1.2 with respect to the QM k(T).

A crucial piece of information regarding the tunelling effect is revealed by the KIE. Fig. 11 displays the comparison of the present theoretical results along with the individual experimental data^13,24 and the global fits by Marshall and Fontjin¹⁵ and Michael,²³ whose expressions are given in the ESI.† The first observation is that the present theoretical KIE is above the global fits to the experimental data, which is not surprising given that the theoretically predicted rate coefficients for O + D₂ are smaller than most of the experimental measurements, whilst they are closer for the O + H₂ reaction, as can be seen in Fig. 7 and 8.

Fig. 11 Comparison of experimental and calculated QM (blue solid line), RPMD (green short-dashed line), and QCT (red dashed–dotted line) kinetic isotopic effect (KIE). The black dashed–dotted line curve is the ratio of the experimental fits from ref. 23 in the 350–2500 K range. The solid dark-yellow line in the 390–1450 K range is the KIE given by Marshall and Fontjin.¹⁵ The black symbols are experimental data from Presser and Gordon¹³ and red symbols are from Zhu et al.²⁴

The KIE predicted by the RPMD calculations is in close agreement with that obtained by the present accurate QM calculations. Only at temperatures below 500 K do the discrepancies become discernible, as shown in the inset of Fig. 11. Even at 300 K, below the crossover temperature for O + H₂, the RPMD KIE is 19 to be compared with 27 as obtained in the QM calculations. As anticipated, the QCT KIE is considerably smaller than those derived from other theoretical approaches. Only above 1000 K does it begin to converge with the QM and RPMD KIEs; at 2500 K the KIE for the three approaches is 1.5. The good agreement between the QM and the RPMD KIEs lends credence to the latter method, which is more efficient than the fully converged QM calculations, especially at relatively high temperatures. The results presented in Fig. 7 and 8, along with the data from Table 1 provide a clear insight into the ability of the RPMD method to account for tunelling in stark contrast to the results of the QCT calculations.

Balakrishnan performed converged QM scattering calculations on the two PESs for the O + H₂ reaction allowing him to calculate the rate coefficients up to 1000 K.³³ Specifically, he used the GLDP fits of the ³A′′ and ³A′ PESs by Rogers et al.²⁵ Balakrishnan assumed that the electronic partition function of the triplet oxygen atom has 9-fold degeneracy, resulting in eqn (16), which is only approximately valid above 1000 K. Once corrected by the appropriate partition function, a good agreement was found with the experimental results as well as with a semi-classical transition state theory.²² However, it should be emphasized that the barriers on the GLDP ³A′′ and ³A′ PESs are smaller than those predicted by the higher-level PESs using in the present work. Furthermore, the two PESs are not degenerate for the collinear arrangement. To the best of our knowledge, there are no previous accurate QM calculations of the rate coefficients for the O + D₂ reaction to compare with the present results.

We have therefore come across two interesting findings: (i) the present state-of-the-art QM results on adiabatic PESs slightly underestimate the experimental results, especially at low temperatures; and (ii) the experimental results agree worse with the QM calculations for O + D₂ than for O + H₂. The fact that QM and RPMD rate coefficients are in good agreement with each other, especially for O + D₂ rules out a systematic error in the QM calculations. As far as the electronic calculations are concerned, the PESs used in this work can be considered the most accurate one available in the literature. Therefore, the only plausible explanation for these findings is that calculations require the inclusion of the ISC between the triplet and singlet PESs in a non-adiabatic treatment.

As commented in the above text, there are only few studies that include the spin–orbit coupling with the singlet PES. The QCT-TSH calculations by Maiti and Schatz³⁶ seemed to lead to larger cross sections for the O + H₂ (v = 0, j = 0) reaction when the ICS was considered. More recent non-adiabatic TD QM by Zhao⁵⁵ for the O + D₂ (v = 0, j = 0) reaction concluded that the spin–orbit influence can be possibly ignored for the title reactions. An interesting aspect is that the non-adiabatic effect seemed to be slightly more important for the reaction with the D₂ isotopic variant (see ref. 55). Unfortunately, the existing calculations are restricted to the ground rovibrational state and made use of the triplet PESs by Rogers et al.²⁵ for both adiabatic and non-adiabatic calculations. Moreover, apart from their qualitative insight, the existing non-adiabatic results should be taken with some reservation since they are based on TSH or in the QM centrifugal sudden approximation. From the comparison between the bulk of the experimental rate coefficients and the current QM results, which can be considered as the most accurate ones carried out on the two separated triplet PESs, it can be concluded that the influence of ISC, although small, cannot be neglected, and that this effect is more important for the O + D₂ reaction. A plausible qualitative explanation for the larger effect of the ISC in the latter reaction can be suggested by the MEP profiles shown in Fig. 1. As discussed in Section 2, the barriers on the two triplet PESs are broader and slightly higher for the O + D₂ isotopic variant. As such, a more effective crossing with the singlet PES can be expected to result in slightly larger cross sections. To assess if ISC could be expected to be more important for O + D₂ than for O + H₂, we calculated the possible crossing of non-reactive trajectories to the singlet PES. We found that singlet–triplet crossing was important only at high temperatures (above 1000 K) and, especially for O + D₂.

Accurate non-adiabatic calculations of the rate coefficients for the title reactions represent an enormous challenge. Even this type of QM calculation at fixed energies, including excited rovibrational states, appears to be a very costly computational endeavour. To the best of our knowledge, the existing non-adiabatic, accurate calculations (converged and including the Coriolis coupling between the different helicities) are limited to the initial rotational state j = 0. In any case, further non-adiabatic, more accurate calculations at fixed energies or using TD-WP for the O + H₂ and O + D₂ reactions would be valuable for assessing the effect of the ISC and would serve to explain the discrepancies between experimental rate coefficients and accurate adiabatic results.

Given the excellent general agreement between RPMD and rigorous QM results, a promising alternative could be the development of non-adiabatic RPMD. This was first presented by Shushkov et al.⁷⁴ for a model system. They used TSH for the instantaneous transitions between PESs and adiabatic evolution of the ring-polymer beads on single PES. Application of this methodology to systems like the O + H₂ reaction might shed light on the importance of the ISC at a much lower computational cost than QM calculations.

5 Conclusions

Extensive quasiclassical (QCT), quantum mechanical (QM) and ring polymer molecular dynamics (RPMD) calculations have been carried out for the reactions of O(³P) with H₂ and D₂ in a range of total energies between 0.37 and 2.5 eV. The QM calculations for the two isotopologues comprise all the total angular momenta and helicity projections necessary for convergence, allowing the determination of rate coefficients up to 2500 K. For the RPMD calculations tests were performed to select the number of beads to ensure convergence. QCT calculations were carried out by means of the cumulative reaction probabilities. The rate coefficient calculations were performed on the two adiabatic ³A′ and ³A′′ PESs, and the respective rate coefficients were combined using an electronic partition function that relates each PES to the O(³P_Ja) atomic states. No attempt has been made to include the spin–orbit crossing or the Renner–Teller coupling between the two PESs.

Converged RPMD results calculated on the two PESs for the two isotopic variants served as a probe of the extent to which this method can tackle the tunnelling effect. Whereas the agreement between the QM and RPMD rate coefficients is almost perfect for the O + D₂ reaction at temperatures above 300 K, for O + H₂ the RPMD predictions are below the QM results at T < 400 K. These temperatures correspond to the respective crossover temperatures.

As a counterpoint, the k(T) predicted by QCT calculations lie well below the experimental or the QM results up to 1000 K. This is also evident in the comparison of the QCT and QM cumulative reaction probabilities at a given total angular momentum. As expected, the discrepancies are smaller for the reaction with D₂ than for O + H₂. If we attribute the discrepancy to tunelling through the barrier, it is found that this effect is also important for the O + D₂ reaction up to 800 K although to a lesser extent than for the reaction with H₂. In contrast, RPMD captures the tunelling effect almost completely for O + D₂ above 300 K and for O + H₂ above 400 K.

The present QM k(T), calculated adiabatically on both PESs and weighted with the electronic partition functions, slightly underestimate the best global fits of the experimental measurements. The discrepancies are rather small for the O + H₂ reaction, in particular for temperatures above 500 K. Interestingly, the differences are larger for the reaction with D₂. Considering that the ab initio PES are more accurate than any previous PESs, and that the RPMD and QM results are generally in excellent agreement, the most plausible explanation for the discrepancies with the experimental results is that a more accurate treatment of the reaction requires a non-adiabatic treatment that include the spin–orbit crossing of the singlet (¹A′) and triplet (³A′ and ³A′′) PESs. Future non-adiabatic calculations for the title reactions will delimit the accuracy of the adiabatic calculations for reactions such as those studied in the present work.

Data availability

Hereunder, we let you know that all data contained in this submission are available.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

The authors gratefully acknowledge grant PID2020-113147GA-I00 and PID2021-122839NB-I00 funded by Spanish Ministry of Science and Innovation (MCIN/AEI/10.13039/MCIN/AEI/10.13039/501100011033). A. V. acknowledges grant no. EDU/1508/2020 (Junta de Castilla y León and European Social Fund). A. Z. also acknowledges grants PID2019-107115GB-C21 and PID2021-122549NB-C21, and the COST Action CA21101 (COSY).

Notes and references

1. D. Fedele, S. Bruderer, E. F. van Dishoeck, J. Carr, G. J. Herczeg, C. Salyk, N. J. Evans, J. Bouwman, G. Meeus, T. Henning, J. Green, J. R. Najita and M. Güdel, Astron. Astrophys., 2013, 559, A77.
2. B. Godard, E. Falgarone and G. Pineau des Forêts, Astron. Astrophys., 2014, 570, A27.
3. R. A. Sultanov and N. Balakrishnan, Astrophys. J., 2005, 629, 305–310.
4. A. Veselinova, M. Agúndez, J. R. Goicoechea, M. Menéndez, A. Zanchet, E. Verdasco, P. G. Jambrina and F. J. Aoiz, Astron. Astrophys., 2021, 648, A76.
5. M. Zannese, B. Tabone, E. Habart, J. R. Goicoechea, A. Zanchet, E. F. van Dishoeck, M. C. van Hemert, J. H. Black, A. G. G. M. Tielens, A. Veselinova, P. G. Jambrina, M. Menendez, E. Verdasco, F. J. Aoiz, L. Gonzalez-Sanchez, B. Trahin, E. Dartois, O. Berne, E. Peeters, J. He, A. Sidhu, R. Chown, I. Schroetter, D. V. D. Putte, A. Canin, F. Alarcon, A. Abergel, E. A. Bergin, J. Bernard-Salas, C. Boersma, E. Bron, J. Cami, D. Dicken, M. Elyajouri, A. Fuente, K. D. Gordon, L. Issa, C. Joblin, O. Kannavou, B. Khan, O. Lacinbala, D. Languignon, R. L. Gal, A. Maragkoudakis, R. Meshaka, Y. Okada, T. Onaka, S. Pasquini, M. W. Pound, M. Robberto, M. Rollig, B. Schefter, T. Schirmer, S. Vicente and M. G. Wolfire, Nat. Astron., 2024, 8, 577–586.
6. M. A. A. Clyne, B. A. Thrush and R. G. W. Norrish, Proc. R. Soc. London, 1963, 275, 544–558.
7. A. A. Westenberg and N. de Haas, J. Chem. Phys., 1967, 46, 490–501.
8. A. A. Westenberg and N. de Haas, J. Chem. Phys., 1967, 47, 4241–4246.
9. A. A. Westenberg and N. de Haas, J. Chem. Phys., 1969, 50, 2512–2516.
10. R. N. Dubinsky and D. J. McKenney, Can. J. Chem., 1975, 53, 3531–3541.
11. I. M. Campbell and B. J. Handy, J. Chem. Soc., Faraday Trans. 1, 1975, 71, 2097–2106.
12. I. M. Campbell and B. J. Handy, J. Chem. Soc., Faraday Trans. 1, 1978, 74, 316–325.
13. N. Presser and R. J. Gordon, J. Chem. Phys., 1985, 82, 1291–1297.
14. J. W. Sutherland, J. V. Michael, A. N. Pirraglia, F. L. Nesbitt and R. B. Klemm, 21st Sym. Int. Combust. Inst., 1986, 929–941.
15. P. Marshall and A. Fontijn, J. Chem. Phys., 1987, 87, 6988–6994.
16. K. Natarajan and P. Roth, Combust. Flame, 1987, 70, 267–279.
17. K. S. Shin, N. Fujii and W. C. J. Gardiner, Chem. Phys. Lett., 1989, 161, 219.
18. D. F. Davidson and R. K. Hanson, Combust. Flame, 1990, 82, 445–447.
19. H. X. Yang, K. S. Shin and W. C. J. Gardiner, Chem. Phys. Lett., 1993, 207, 69–74.
20. S. O. Ryu, S. M. Hwang and M. J. Rabinowitz, Chem. Phys. Lett., 1995, 242, 279–284.
21. D. L. Baulch, C. T. Bowman, C. J. Cobos, R. A. Cox, T. Just, J. A. Kerr, M. J. Pilling, D. Stocker, J. Troe, W. Tsang, R. W. Walker and J. Warnatz, J. Phys. Chem. Ref. Data, 2005, 34, 757–1397.
22. T. L. Nguyen and J. F. Stanton, J. Phys. Chem. A, 2014, 118, 4918–4928.
23. J. V. Michael, J. Chem. Phys., 1989, 90, 189–198.
24. Y. F. Zhu, S. Arepalli and R. J. Gordon, J. Chem. Phys., 1989, 90, 183–188.
25. S. Rogers, D. Wang, A. Kuppermann and S. Walsh, J. Phys. Chem. A, 2000, 104, 2308.
26. J. Brandao, C. Mogo and B. C. Silva, J. Chem. Phys., 2004, 121, 8861–8868.
27. A. Zanchet, M. Menéndez, P. G. Jambrina and F. J. Aoiz, J. Chem. Phys., 2019, 151, 094307.
28. D. J. Garton, T. K. Minton, B. Maiti, D. Troya and G. C. Schatz, J. Chem. Phys., 2003, 118, 1585–1588.
29. D. J. Garton, A. L. Brunsvold, T. K. Minton, D. Troya, B. Maiti and G. C. Schatz, J. Phys. Chem. A, 2006, 110, 1327–1341.
30. S. A. Lahankar, J. Zhang, K. G. McKendrick and T. K. Minton, Nat. Chem., 2013, 5, 315.
31. S. A. Lahankar, J. Zhang, T. K. Minton and K. G. McKendrick, J. Am. Chem. Soc., 2014, 136, 12371–12384.
32. N. Balakrishnan, J. Chem. Phys., 2003, 119, 195–199.
33. N. Balakrishnan, J. Chem. Phys., 2004, 121, 6346–6352.
34. P. F. Weck, N. Balakrishnan, J. Brandao, C. Rosa and W. Wang, J. Chem. Phys., 2006, 124, 074308.
35. M. Braunstein, S. Adler-Golden, B. Maiti and G. C. Schatz, J. Chem. Phys., 2004, 120, 4316–4322.
36. B. Maiti and G. C. Schatz, J. Chem. Phys., 2003, 119, 12360.
37. P. G. Jambrina, A. Zanchet, J. Aldegunde, M. Brouard and F. J. Aoiz, Nat. Commun., 2016, 7, 13439.
38. I. R. Craig and D. E. Manolopoulos, J. Chem. Phys., 2005, 123, 034102.
39. R. Collepardo-Guevara, Y. V. Suleimanov and D. E. Manolopoulos, J. Chem. Phys., 2009, 130, 174713.
40. Y. V. Suleimanov, R. Collepardo-Guevara and D. E. Manolopoulos, J. Chem. Phys., 2011, 134, 044131.
41. Y. V. Suleimanov, J. W. Allen and W. H. Green, Comput. Phys. Commun., 2013, 184, 833–840.
42. R. Pérez de Tudela, Y. V. Suleimanov, J. O. Richardson, V. Sáez Rábanos, W. H. Green and F. J. Aoiz, J. Phys. Chem. Lett., 2014, 5, 4219–4224.
43. Y. Li, Y. V. Suleimanov, W. H. Green and H. Guo, J. Phys. Chem. A, 2014, 118, 1989–1996.
44. Y. Suleimanov, F. J. Aoiz and H. Guo, J. Phys. Chem. A, 2016, 120, 8488–8502.
45. M. Menéndez, P. G. Jambrina, A. Zanchet, E. Verdasco, Y. V. Suleimanov and F. J. Aoiz, J. Phys. Chem. A, 2019, 123, 7920–7931.
46. R. Pérez de Tudela, F. J. Aoiz, Y. V. Suleimanov and D. E. Manolopoulos, J. Phys. Chem. Lett., 2012, 3, 493–497.
47. Y. Li, Y. V. Suleimanov, M. Yang, W. H. Green and H. Guo, J. Phys. Chem. Lett., 2013, 4, 48–52.
48. Y. V. Suleimanov, R. P. de Tudela, P. G. Jambrina, J. F. Castillo, V. Sáez-Rábanos, D. E. Manolopoulos and F. J. Aoiz, Phys. Chem. Chem. Phys., 2013, 15, 3655–3665.
49. P. G. Jambrina, A. Zanchet, M. Menéndez, V. J. Herrero and F. J. Aoiz, Phys. Chem. Chem. Phys., 2019, 21, 25389–25396.
50. A. Veselinova, M. Menéndez, L. González-Sánchez, A. Zanchet, F. J. Aoiz and P. G. Jambrina, Phys. Chem. Chem. Phys., 2024, 26, 6752–6762.
51. G. A. Natanson, B. C. Garrett, T. N. Truong, T. Joseph and D. G. Truhlar, J. Chem. Phys., 1991, 94, 7875–7892.
52. J. Corchado, Y. Chuang, P. Fast, W.-P. Hu, Y.-P. Liu, G. Lynch, K. Nguyen, C. Jackels, A. Fernandez Ramos, B. Ellingson, B. Lynch, J. Zheng, V. Melissas, J. Villà, I. Rossi, E. Coitino, J. Pu, T. Albu, R. Steckler, B. Garrett, A. Isaacson and D. Truhlar, POLYRATE-v. 9.7, University of Minnesota, Minneapolis, Minnesota 55455, 2007.
53. J. O. Richardson and S. C. Althorpe, J. Chem. Phys., 2009, 131, 214106.
54. M. R. Hoffman and G. C. Schatz, J. Chem. Phys., 2000, 113, 9456.
55. J. Zhao, J. Chem. Phys., 2013, 138, 134309.
56. A. F. Wagner, J. M. Bowman and L. B. Harding, J. Chem. Phys., 1985, 82, 1866–1872.
57. W. H. Miller, Acc. Chem. Res., 1976, 9, 306–312.
58. W. H. Miller, Acc. Chem. Res., 1993, 26, 174–181.
59. S. Chapman, B. C. Garrett and W. H. Miller, J. Chem. Phys., 1975, 63, 2710–2716.
60. D. C. Chatfield, R. S. Friedman, D. W. Schwenke and D. G. Truhlar, J. Chem. Phys., 1992, 96, 2414–2421.
61. D. C. Chatfield, S. L. Mielke, T. C. Allison and D. G. Truhlar, J. Chem. Phys., 2000, 112, 8387–8408.
62. F. J. Aoiz, M. Brouard, C. J. Eyles, J. F. Castillo and V. Sáez-Rábanos, J. Chem. Phys., 2006, 125, 144105.
63. P. G. Jambrina, M. Lara, M. Menéndez, J.-M. Launay and F. J. Aoiz, J. Chem. Phys., 2012, 137, 164314.
64. F. J. Aoiz, V. J. Herrero, M. P. de Miranda and V. Sáez-Rábanos, Phys. Chem. Chem. Phys., 2007, 9, 5367–5373.
65. P. G. Jambrina, F. J. Aoiz, C. J. Eyles, V. J. Herrero and V. Sáez-Rábanos, J. Chem. Phys., 2009, 130, 184303.
66. C. H. Bennett, Molecular Dynamics and Transition State Theory: The Simulation of Infrequent Events, American Chemical Society, 1977, vol. 46, ch. 4, pp. 63–97.
67. D. Chandler, J. Chem. Phys., 1978, 68, 2959.
68. D. Skouteris, J. F. Castillo and D. E. Manolopoulos, Comput. Phys. Commun., 2000, 133, 128–135.
69. F. J. Aoiz, M. Brouard, C. J. Eyles, J. F. Castillo and V. Sáez-Rábanos, J. Chem. Phys., 2006, 125, 144105.
70. P. G. Jambrina, E. García, V. J. Herrero, V. Sáez-Rábanos and F. J. Aoiz, Phys. Chem. Chem. Phys., 2012, 14, 14596.
71. J. Aldegunde, P. G. Jambrina, E. García, V. J. Herrero, V. Sáez-Rábanos and F. J. Aoiz, Mol. Phys., 2013, 21, 3169.
72. D. Chatfield, S. Mielke, T. Allison and D. Truhlar, J. Chem. Phys., 2000, 112, 8387.
73. S. Javoy, V. Naudet, S. Abid and C. E. Paillard, Int. J. Chem. Kinet., 2000, 32, 686–695.
74. P. Shushkov, R. Li and J. C. Tully, J. Chem. Phys., 2012, 137, 22A549.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp01711k

---