Atom tunnelling in the reaction NH₃⁺ + H₂ → NH₄⁺ + H and its astrochemical relevance

Abstract

The title reaction is involved in the formation of ammonia in the interstellar medium. We have calculated thermal rates including atom tunnelling using different rate theories. Canonical variational theory with microcanonically optimised multidimensional tunnelling was used for bimolecular rates, modelling the gas-phase reaction and also a surface-catalysed reaction of the Eley–Rideal type. Instanton theory provided unimolecular rates, which model the Langmuir–Hinshelwood type surface reaction. The potential energy was calculated on the CCSD(T)-F12 level of theory on the fly. We report thermal rates and H/D kinetic isotope effects. The latter have implications for observed H/D fractionation in molecular clouds. Tunnelling causes rate constants to be sufficient for the reaction to play a role in interstellar chemistry even at cryogenic temperature. We also discuss intricacies and limitations of the different tunnelling approximations to treat this reaction, including its pre-reactive minimum.

---

1 Introduction

The possible routes to NH₃ formation in the interstellar medium are of significant interest in astrochemistry.¹ One promising exothermic reaction is:

NH₃⁺ + H₂ → NH₄⁺ + H. (R1)

The resulting NH₄⁺ species can generate ammonia and a hydrogen atom by recombining with an electron and subsequently dissociating a hydrogen atom.¹(R1) is hindered by a barrier of 2.1 kcal mol⁻¹ (8.8 kJ mol⁻¹).² The competitive dissociative recombination of NH₃⁺ with an electron (formation of NH + H₂) must be taken into account.¹ This means that (R1) is most relevant in dark clouds where the concentration of free electrons is low.^3,4

Early experiments reported a low rate constant of k < 5 × 10⁻¹³ cm³ s⁻¹ for (R1) at room temperature.^2,5 Smith and Adams measured the rate constant of (R1) down to approximately 80 K and found it to decrease with T.⁶ However, it is known that some exothermic ion–molecule reactions show an increase of reactivity by reducing the temperature.⁷ Experiments carried out at very low temperatures (20 to 11 K) by Luine and Dunn showed an increase in the reaction rate of (R1).⁸ They proposed a pre-reactive complex (NH₃⁺ − H₂) with an increased lifetime at lower temperature.⁸ The pressure-dependent results obtained by Böhringer supported this idea.⁹ In addition, the isotope effects detected by Smith and Adams,⁶ as well as Barlow and Dunn,¹⁰ also confirmed the presence of an entrance channel complex, from which tunnelling through the potential energy barrier could take place.^10,11 Nevertheless, for a full understanding of reaction (R1) quantum chemical calculations are needed.

In 1991 Herbst et al.¹² mapped the potential energy surface (PES) for reaction (R1) at the MP4/6-311++G(3df,3dp)//MP2/6-31G(d,p) level.^13–17 They found a weakly bound entrance channel complex (RS), as well as an exit channel complex (PS). Rate constants were calculated including a one-dimensional tunnelling approximation via an Eckart barrier. Their results qualitatively reproduced the variation of the rate coefficients with the temperature quite well, but they were quantitatively over 1.5 orders of magnitude lower than the experimental ones. Only by artificially reducing the potential energy barrier by 1 kcal mol⁻¹ and increasing the curvature at the transition state, could they improve the agreement. Their predicted rate constants still remained lower than the experimental values; however, their calculations confirmed that tunnelling plays a significant role in reaction (R1) at low temperatures.

The main goal of this study is to calculate tunnelling rates with more accurate methods than the ones found in the literature.^12,18 Firstly, a careful study of the PES comparing high-level methods has been done. Secondly, bimolecular rate constants were calculated with transition state theory including tunnelling through several approaches,¹⁹ and unimolecular rate constants were additionally computed with the instanton method.^20–37 Finally, calculations regarding the kinetic isotope effect were carried out.

2 Computational details

2.1 Optimisation procedure

All the geometry optimisations were performed with DL-FIND³⁸ in ChemShell.^39,40 For the electronic structure computations (energy and gradient) we have employed Molpro 2012.⁴¹ Starting points for the geometries were obtained from Herbst et al.¹²

We carried out several single point calculations using the geometries obtained at the MCSCF/cc-pVTZ level.^42,43 We used the MRCISD(Q)-F12/cc-pVQZ-F12//MCSCF/cc-pVTZ level^42–45 as a reference and compared the relative energies of the RS, TS and PS of reaction (R1) (see the ESI).† Table S1 and Fig. S1† show that the relative energies at the UCCSD(T)-F12a/cc-pVDZ-F12//MCSCF/cc-pVTZ level^{42,43,45–47} agree well with the reference results. In addition, the analysis of the T₁-diagnostic values for the coupled cluster calculations allows us to conclude that a single-reference wave function can be used to describe this system, since the T₁ values were always lower than 0.02.^48,49 For all further calculations we have used the UCCSD(T)-F12a/cc-pVDZ-F12 level,^45–47 which represents a reasonable compromise between accuracy and computational efficiency.

The default thresholds of the programs were used in general. However, for the optimisation of the TS, we employed stricter convergence thresholds in order to reproduce the C_3v symmetry for the TS and the correct eigenvalue spectrum of the Hessian. We changed the maximum gradient of the geometry optimisation in DL-FIND to 5 × 10⁻⁶ atomic units. Additionally, the threshold of the SCF and UCCSD energy convergence was set to 10⁻¹⁰ Hartree for the optimisation of the stationary points.

2.2 Variational transition state theory calculations

We employed TST by using the Polyrate code version 2010.⁵⁰ We used conventional TST as well as the canonical variational TST (CVT) (the position of the TS is located variationally for every temperature). Furthermore, the transmission factors were tested with several tunnelling methods (in combination with the TST types): zero-curvature tunnelling method (ZCT), small-curvature tunnelling method (SCT), large-curvature tunnelling method (LCT), canonically-optimised multidimensional tunnelling method (COMT) and the microcanonically-optimised multidimensional tunnelling method (μOMT).^19,51 In the COMT method, the transmission factor corresponds to the largest thermally averaged value between the SCT and LCT for a given T. In the μOMT method, the transmission factor is computed from LCT and SCT results by taking the larger of the two for each energy. For (R1), SCT always results in larger transmission factors than LCT.

The minimum energy path (MEP) was calculated at the UCCSD(T)-F12/cc-pVDZ-F12 level. The Page–McIver algorithm⁵² (implemented in Polyrate) was used to follow the reaction path, with the Hessian recalculated every 9 steps. To be able to use energy and gradients in Polyrate directly from Molpro, we interfaced Polyrate to ChemShell. We improved the accuracy by using the ISPE approach (interpolated single-point energies):^50,53 the energy of the stationary points and some points along the MEP were recalculated with a larger basis set (UCCSD(T)-F12/cc-pVTZ-F12, see Fig. S1†). Note that the calculations have been done considering both the wells in the reactant region (RS) and the product side (PS).

Therefore, for the rate constant calculations including tunnelling by multidimensional transmission (MT) factors, we considered a pre-reactive complex called RS (see details in Sections 3.1 and 3.2).

An important detail in the rate calculations is to include the symmetry factor (σ) which is related to the symmetry index of the reactants and saddle point:⁵⁴

σ_rot represents the rotational symmetry number and n_* the number of optical isomers of the molecule. For (R1): σ_rot(NH₃⁺, D_3h) = 6; σ_rot(H₂, D_∞h) = 2 and σ_rot(TS, C_3v) = 3. Consequently, σ(R1) = 4. For the unimolecular reaction from the pre-reactive complex RS, σ(RS, C_s), the overall symmetry factor is σ = 1/3.

2.3 Instanton calculations

We also used the instanton method^20–37 based on Feynman path integral theory using the semiclassical approximation to obtain the rate constants. For a given temperature, it provides the most probable tunnelling path, the instanton. Instanton theory is applicable whenever the temperature is low enough for the instanton to spread out. At higher temperatures, the instanton collapses to a point which renders the theory inapplicable. For many barrier shapes⁵⁵ this collapse happens at the crossover temperature T_c,⁵⁶with Ω being the barrier frequency (the absolute value of the imaginary frequency corresponding to the transition mode) and k_B corresponding to Boltzmann's constant.

Instanton paths were optimised via a quasi Newton–Raphson method as described previously.^34,35 Energies and gradients were provided by Molpro, but instanton optimisations and Hessian calculations via finite differences were done in DL-FIND. The instanton path was discretised using 154 images.

3 Results

3.1 PES calculations

Optimisation of the geometries for the stationary points in reaction (R1) was performed at the UCCSD(T)-F12/cc-pVDZ-F12 level. We found a pre-reactive complex (RS) with the dihydrogen molecule almost parallel to the plane of the NH₃ fragment with N–H distances of 2.619 Å and 2.605 Å, see Fig. 1. The RS is 8.57 kJ mol⁻¹ below the energy of the separated reactants (3.25 kJ mol⁻¹ when the ZPE is included). In the C_3v-symmetric transition state structure, H₂ attacks the nitrogen atom head-on, with an N–H distance of 1.647 Å and a H–H distance of 0.786 Å. Compared to the separated reactants, the barrier is 6.03 kJ mol⁻¹ (12.13 kJ mol⁻¹ when the ZPE is included). Cartesian coordinates of all stationary points are available in the ESI.† A comparison to the previous potential energy surface on the MP4/6-311++G(3df,3dp)//MP2/6-31G(d,p) level¹² is given in Fig. 1. It is clear that there is a significant pre-reactive minimum, which can be expected to be populated at low temperature and high-enough pressure to allow for thermal equilibration within the lifetime of the complex.

Fig. 1 Energy profile of reaction (R1) computed here (blue) compared to a previous one¹² (red). The relative energies (kJ mol⁻¹) including ΔZPEs of the RS, TS, PS and products were calculated with respect to the reactants.

3.2 Bimolecular rate constants

To obtain bimolecular reaction rates, we have employed TST (also in its canonical version) in combination with more or less sophisticated tunnelling methods. Note that the variational determination of the transition state (model CVT vs. TST) has a negligible impact on the reaction rates. As has been said, the small curvature results, SCT, always give higher transmission factors than their LCT counterparts, so the CVT/COMT (and μOMT) results are equal to the CVT/SCT values; furthermore, the μOMT and COMT results are coincident. A full list of results is given in Table S2 of the ESI.†

In addition, we have employed the ISPE correction to improve our results, see Table S3 of the ESI.† The energies of the stationary points and some points along the path have been recalculated with the larger cc-pVTZ-F12 basis set. This lowers the barrier by 0.9 kJ mol⁻¹. ISPE is not available for CVT/LCT, but we assume that the rates remain smaller than those of the CVT/SCT approach.

The comparison of the TST and CVT/SCT rate constant values (including the ISPE correction) (see Fig. 2) shows that the reaction proceeds mostly classically above approximately 200 K. We can observe a monotonic increase in rate with temperature at T > 200 K. Therefore, in that T-range, the CVT/SCT values are practically coincident with the TST values. However, tunnelling dominates below the crossover temperature of T_c = 148.5 K. The rate constant is minimal at 70 K and slightly increases at even lower temperature. In full thermal equilibrium, this could be caused by the population of a bound pre-reactive complex (RS), which increases the attempt frequency for the transition.⁵⁷ At low pressure, however, one would expect the RS to decay to the products or back to the reactants before thermal equilibration. Note that the Polyrate code⁵⁰ does not compute the capture rates (formation of the RS complex from NH₃⁺ + H₂). However, the presence of the RS pre-reactive complex is taken into account (and it has a huge impact) during the MEP calculation. This MEP is required for the CVT rate constant and MT tunnelling factor calculations. A good description of the path is necessary because it has a large impact in the selection of the turning points during the tunnelling calculation.

Fig. 2 Comparison of calculated bimolecular rate constants to experimental values and ref. 9–11. The inset shows the same graphs at different scaling and demonstrates the effect of tunnelling.

3.2.1 Comparison to experiments. The rate constant of the title reaction was measured several times in different temperature ranges.^2,5,6,8–11 We have summarised these results in Fig. 2 and compared them to our results and the values used in the astrochemical databases UMIST and KIDA.^58,59

Before comparing to theory, it makes sense to compare the experimental results among each other. At 20 K, the rate constant obtained by Böhringer⁹ is 2.4 times larger than the one by Luine and Dunn,⁸ and 3.2 times larger than the one by Barlow and Dunn.¹⁰ Our results, extrapolated to 20 K, can be expected to be close to those by Luine and Dunn, as well as those by Barlow and Dunn, see Fig. 2. Consequently, one may speculate that the higher-temperature values by Böhringer somewhat overestimate the rate constant as well. Our high-T rate constants agree quite well with the experimental ones. There is also agreement that the rate constant has a minimum in the temperature range between about 100 K and 30 K. Below that, the rate constant increases again, which is found by both theory and experiment.

The rate expressions used in the kinetic databases KIDA and UMIST agree remarkably well with our results at low temperature. Also at high temperature, the expression used in KIDA agrees well with our data, while UMIST somewhat underestimates the rate constant. At intermediate temperatures, where the rate constant is minimal, KIDA, UMIST and our results are always within a factor of 3.

3.2.2 Instanton calculations: intricacies and limitations. The crossover temperature (T_c) for this reaction is 148.5 K. Instantons were calculated in the range of 145 to 50 K. The PES for this reaction leads to several difficulties for semiclassical instanton theory. The first is the difference in symmetry between the TS and the RS. The instanton at 145 K shows a relatively small delocalisation and retains the C_3v symmetry of the TS, see Fig. 3. However, by lowering T to just 140 K, the delocalisation increases and the C_3v-symmetric instanton becomes unstable. A symmetry breaking takes place; first-order saddle point paths in the Euclidean action at lower temperature have the C_s symmetry of the RS. Similar temperature-dependent symmetry breakings have also been observed in other cases.⁶⁰ In the temperature range of the instability between the instantons of different symmetry, the vibrational frequency of the instanton associated with the instability becomes small, leading to unrealistically high predicted rate constants. Combined with the previous finding⁶¹ that instanton theory often overestimates the rate constant close to T_c, instanton rates are probably only reliable well below T_c.

Fig. 3 Left: instanton at 145 K (C_3v). Right: instanton at 50 K (C_s). Images are coloured from the RS-end (red) to the PS-end (blue).

The intrinsic reaction coordinate (IRC) is C_3v-symmetric at the TS and quite some way down into the reactant well, see Fig. 4. There, a dihedral angle which shows the distortion from C_3v symmetry is plotted against the distance of the attacking hydrogen atom to the nitrogen atom for the IRC and several instanton paths. The instanton paths for T ≤ 140 K differ substantially from the IRC. This phenomenon is known as corner cutting by the tunnelling motion.⁶² It can also be seen in Fig. 4 that the lower the temperature gets, the longer the instanton paths become. The fact that they end in similar areas of configuration space on the product side is due to the fact that both ends of the instanton path have the same energy. So they end when the energy of the PS side drops below the energy of the RS. Here, it should be noted that only half of each instanton path is visible in Fig. 4: the paths are closed on themselves and follow the line shown in Fig. 4 backward and forwards.

Fig. 4 Symmetry breaking in instanton calculations: the plot shows instantons, the IRC and the stationary points projected onto a dihedral angle and the N–H bond length. C_3v symmetry corresponds to a dihedral angle of 0. Equipotential lines at −3, −1, 1, … kJ mol⁻¹ are shown as thin dashed lines.

A different problem in instanton calculations concerns the evaluation of bimolecular rates. Since the barrier is rather low (almost submerged) compared to the RS, all instantons for T < 145 K have a tunnelling energy (the energy of the turning points of the instanton) below the energy of the separated reactants. Therefore, bimolecular rate constants are not available from instanton theory. It would be possible to calculate temperature-dependent bimolecular rates by calculating energy-dependent transmission coefficients and taking the Boltzmann average of these based on the energy of the separated reactants. This was outside the scope of this paper, though. Instead, we calculated unimolecular rate constants.

3.3 Unimolecular rate constants

The title reaction is expected to be relevant on the surface of interstellar dust grains, where unimolecular rate constants are relevant. Different types of surface reactions can be distinguished: if one reaction partner is adsorbed on the surface and the other one comes in from the gas phase, one speaks of Eley–Rideal-type reactions. Those are described by bimolecular rate constants. If both partners are adsorbed on the surface and move diffusively until they meet, a Langmuir–Hinshelwood mechanism takes place. The corresponding rate constants can be approximated by unimolecular gas-phase rate constants from the RS minimum. Adsorption of one or both reactants on a dust grain can, of course, change the PES and especially the barrier height. For a reasonably inert surface this effect might be neglected in a first approximation. The other effect of adsorption is the greatly facilitated thermal equilibration. Excess heat produced by formation of a pre-reactive complex or by the reaction itself is rapidly dissipated into the bulk, which justifies the use of thermal rate constants.

Unimolecular rate constants were calculated with instanton theory and CVT/SCT theory. Results are shown in Fig. 5 and Table S4 of the ESI.† At high temperatures, close to T_c, instanton theory is known to overestimate the rate constant.⁶¹ This is also seen here. Between 120 K and 80 K there is reasonable agreement between the methods. At low temperature, instanton theory predicts a temperature-independent unimolecular rate constant, while the one predicted by CVT/SCT continues to decrease. Below 80 K, one can assume instanton rate constants to be more trustworthy than the SCT data because instanton theory uses an optimised tunnelling path. As shown in Fig. 4, the tunnelling path differs significantly from the IRC at low temperature. It should be noted that the instanton optimisations for T < 60 K are not perfectly converged, which causes the slight wiggles in Fig. 5.

Fig. 5 Unimolecular rate constants calculated with instanton theory compared to CVT/SCT.

3.4 Kinetic isotope effects

We have also studied the effect of replacing H₂ by HD and D₂ in (R1), the kinetic isotope effect (KIE). This influences the deuterium fractionation of molecules observed in the interstellar medium.^4,63 The KIE is defined as the ratio between the rate constant of a reaction with a lighter isotope and the rate constant of the reaction with a heavier isotope. The following reactions were investigated:

NH₃⁺ + HD → NH₄⁺ + D, (R2)

NH₃⁺ + DH → NH₃D⁺ + H, (R3)

NH₃⁺ + D₂ → NH₃D⁺ + D. (R4)

Note that (R4) has been studied experimentally by Adams and Smith, as well as Barlow and Dunn.^10,11(R2) and (R3) both represent the reaction of NH₃⁺ with HD, but the products are different. From the resulting rate constant, one can calculate a branching ratio.

We have computed the rates of the deuterated systems with Polyrate, which includes a scaling procedure to recalculate the rates from the non-deuterated reactions.^64,65 TST in combination with the SCT method (TST/SCT), including the ISPE correction, was used. The approach is almost equivalent to the one employed in Fig. 2, with only CVT being replaced by TST. For the non-deuterated case, these two methods resulted in almost identical rate constants, so we assume that they are very similar for the deuterated cases as well. The most important results are shown in Fig. 6. Furthermore, some KIE values are listed in Table 1.

Fig. 6 Effect of deuteration on the reaction rates, calculated with TST/SCT, ISPE correction at the UCCSD(T)-F12/cc-pVTZ-F12 level. Reactions (R1) (noted as H₂, black), (R2) (noted as HD, red), (R3) (noted as DH, green) and (R4) (noted as D₂, blue) are plotted.

Table 1 KIEs values for reactions (R2)–(R4) and branching ratio (BR) for (R2)vs.(R3)

T, K	KIE (HD)	KIE (DH)	BR (%)	KIE (D2)	KIE (D2) ref. 11
25	3.8	10.0	72	11.2	—
35	3.7	9.4	72	10.3	—
45	3.5	9.0	72	9.5	—
55	3.4	8.6	72	8.8	—
65	3.3	8.3	72	8.2	—
75	3.2	7.9	71	7.7	—
85	3.1	7.6	71	7.1	11.3
100	3.0	7.2	71	6.5	—
110	3.0	6.9	70	6.1	—
125	2.9	6.5	69	5.5	—
145	2.9	6.1	68	5.0	—
170	2.8	5.7	67	4.4	—
200	2.8	5.3	65	3.9	3.3
300	2.7	4.5	63	2.9	1.1

---

4 Conclusions

We have calculated bimolecular and unimolecular rate constants for the title reaction with different rate theories. Our results agree to a reasonable degree with available experimental data. They also confirm the fitted rate expressions currently used in astrochemical models. It is clear that at low temperatures (<70 K), the bimolecular rate constant increases with decreasing temperature. This unusual behaviour is caused by the quantum mechanical tunnel effect. Without tunnelling (see the TST graph in Fig. 2), the rate constant decreases strongly at lower temperature.

It is clear that atom tunnelling plays a key role in the title reaction. The shape of the potential energy surface, including an almost submerged barrier and a tunnelling path which differs significantly from the minimum-energy path, lead to challenges in the theoretical description. For the title reaction, instanton theory is only applicable well below the crossover temperature. In that region, however, one can expect it to provide more accurate rate constants than CVT/SCT. The latter assumes a tunnelling path close to the minimum-energy path by construction. It accounts for some corner cutting, but can not adapt to qualitatively different tunnelling paths. Nevertheless, the agreement with experiment is quite satisfactory.

Acknowledgements

This work was financially supported by the German Research Foundation (DFG) within the Cluster of Excellence in Simulation Technology (EXC 310/2) at the University of Stuttgart as well as by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 646717, TUNNELCHEM).

References

1. A. G. G. M. Tielens, Rev. Mod. Phys., 2013, 85, 1021–1081.
2. F. C. Fehsenfeld, W. Lindinger, A. L. Schmeltekopf, D. L. Albritton and E. E. Ferguson, J. Chem. Phys., 1975, 62, 2001–2003.
3. V. G. Anicich and W. T. Huntress Jr, Astron. Astrophys., Suppl. Ser., 1986, 62, 553–672.
4. F. Daniel, C. Rist, A. Faure, E. Roueff, M. Gérin, D. C. Lis, P. Hily-Blant, A. Bacmann and L. Wiesenfeld, Mon. Not. R. Astron. Soc., 2016, 457, 1535–1549.
5. J. K. Kim, L. P. Theard and W. T. Huntress, J. Chem. Phys., 1975, 62, 45–52.
6. D. Smith and N. G. Adams, Mon. Not. R. Astron. Soc., 1981, 197, 377–384.
7. E. Ferguson, Annu. Rev. Phys. Chem., 1975, 26, 17–38.
8. J. A. Luine and G. H. Dunn, Astrophys. J., 1985, 299, L67–L70.
9. H. Böhringer, Chem. Phys. Lett., 1985, 122, 185–189.
10. S. Barlow and G. Dunn, Int. J. Mass Spectrom. Ion Processes, 1987, 80, 227–237.
11. N. G. Adams and D. Smith, Int. J. Mass Spectrom. Ion Processes, 1984, 61, 133–139.
12. E. Herbst, D. J. DeFrees, D. Talbi, F. Pauzat, W. Koch and A. D. McLean, J. Chem. Phys., 1991, 94, 7842–7849.
13. C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618–622.
14. R. Krishnan and J. A. Pople, Int. J. Quantum Chem., 1978, 14, 91–100.
15. R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654.
16. T. Clark, J. Chandrasekhar, G. W. Spitznagel and P. V. R. Schleyer, J. Comput. Chem., 1983, 4, 294–301.
17. M. J. Frisch, J. A. Pople and J. S. Binkley, J. Chem. Phys., 1984, 80, 3265–3269.
18. J. Ischtwan, B. J. Smith, M. A. Collins and L. Radom, J. Chem. Phys., 1992, 97, 1191–1210.
19. A. Fernandez-Ramos, B. A. Ellingson, B. C. Garrett and D. G. Truhlar, in Rev. Comput. Chem., ed. K. B. Lipkowitz and T. R. Cundari, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2007, vol. 23, ch. Variational Transition State Theory with Multidimensional Tunneling, pp. 125–232.
20. J. S. Langer, Ann. Phys., 1967, 41, 108.
21. J. S. Langer, Ann. Phys., 1969, 54, 258–275.
22. W. H. Miller, J. Chem. Phys., 1975, 62, 1899.
23. S. Coleman, Phys. Rev. D: Part. Fields, 1977, 15, 2929.
24. C. G. Callan Jr and S. Coleman, Phys. Rev. D: Part. Fields, 1977, 16, 1762.
25. E. Gildener and A. Patrascioiu, Phys. Rev. D: Part. Fields, 1977, 16, 423.
26. I. Affleck, Phys. Rev. Lett., 1981, 46, 388–391.
27. S. Coleman, Nucl. Phys. B, 1988, 298, 178.
28. P. Hänggi, P. Talkner and M. Borkovec, Rev. Mod. Phys., 1990, 62, 251.
29. V. A. Benderskii, D. E. Makarov and C. A. Wight, Adv. Chem. Phys., 1994, 88, 55.
30. M. Messina, G. K. Schenter and B. C. Garrett, J. Chem. Phys., 1995, 103, 3430.
31. J. O. Richardson and S. C. Althorpe, J. Chem. Phys., 2009, 131, 214106.
32. M. Kryvohuz, J. Chem. Phys., 2011, 134, 114103.
33. S. C. Althorpe, J. Chem. Phys., 2011, 134, 114104.
34. J. B. Rommel, T. P. M. Goumans and J. Kästner, J. Chem. Theory Comput., 2011, 7, 690.
35. J. B. Rommel and J. Kästner, J. Chem. Phys., 2011, 134, 184107.
36. M. Kryvohuz, J. Phys. Chem. A, 2014, 118, 535–544.
37. J. O. Richardson, J. Chem. Phys., 2016, 144, 114106.
38. J. Kästner, J. M. Carr, T. W. Keal, W. Thiel, A. Wander and P. Sherwood, J. Phys. Chem. A, 2009, 113, 11856.
39. ChemShell, A Computational Chemistry Shell, see http://www.chemshell.org.
40. S. Metz, J. Kästner, A. A. Sokol, T. W. Keal and P. Sherwood, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 101.
41. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, W. Györffy, D. Kats, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. Köppl, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, D. P. O'Neill, P. Palmieri, D. Peng, K. Pflüger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson and M. Wang, MOLPRO, version 2012, a package of ab initio programs, 2012, see https://www.molpro.net/.
42. H.-J. Werner and P. J. Knowles, J. Chem. Phys., 1985, 82, 5053–5063.
43. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007.
44. T. Shiozaki and H.-J. Werner, Mol. Phys., 2013, 111, 607–630.
45. K. A. Peterson, T. B. Adler and H.-J. Werner, J. Chem. Phys., 2008, 128, 084102.
46. T. B. Adler, G. Knizia and H.-J. Werner, J. Chem. Phys., 2007, 127, 221106.
47. G. Knizia, T. B. Adler and H.-J. Werner, J. Chem. Phys., 2009, 130, 054104.
48. T. J. Lee and P. R. Taylor, Int. J. Quantum Chem., 1989, 36, 199–207.
49. T. J. Lee, Chem. Phys. Lett., 2003, 372, 362–367.
50. J. Zheng, S. Zhang, B. J. Lynch, J. C. Corchado, Y.-Y. Chuang, P. L. Fast, W.-P. Hu, Y.-P. Liu, G. C. Lynch, K. A. Nguyen, C. F. Jackels, A. Fernandez-Ramos, B. A. Ellingson, V. S. Melissas, J. Villa, I. Rossi, E. L. Coitino, J. Pu, T. V. Albu, A. Ratkiewicz, R. Steckler, B. C. Garrett, A. D. Isaacson and D. G. Truhlar, POLYRATE-version 2010, University of Minnesota, Minneapoli, USA, 2010.
51. A. Fernandez-Ramos, J. A. Miller, S. J. Klippenstein and D. G. Truhlar, Chem. Rev., 2006, 106, 4518–4584.
52. M. Page and J. W. McIver, J. Chem. Phys., 1988, 88, 922–935.
53. Y.-Y. Chuang, J. C. Corchado and D. G. Truhlar, J. Phys. Chem. A, 1999, 103, 1140–1149.
54. A. Fernández-Ramos, B. A. Ellingson, R. Meana-Pañeda, J. M. C. Marques and D. G. Truhlar, Theor. Chem. Acc., 2007, 118, 813–826.
55. S. Álvarez-Barcia, J. R. Flores and J. Kästner, J. Phys. Chem. A, 2014, 118, 78.
56. M. J. Gillan, J. Phys. C: Solid State Phys., 1987, 20, 3621.
57. J. Meisner and J. Kästner, Angew. Chem., Int. Ed., 2016, 55, 5400–5413.
58. UMIST database, http://www.udfa.ajmarkwick.net.
59. KIDA database, http://www.kida.obs.u-bordeaux1.fr.
60. J. Meisner, J. B. Rommel and J. Kästner, J. Comput. Chem., 2011, 32, 3456.
61. T. P. M. Goumans and J. Kästner, J. Phys. Chem. A, 2011, 115, 10767.
62. R. A. Marcus and M. E. Coltrin, J. Chem. Phys., 1977, 67, 2609.
63. E. Roueff, J. C. Loison and K. M. Hickson, Astron. Astrophys., 2015, 576, A99.
64. P. L. Fast, J. C. Corchado and D. G. Truhlar, J. Chem. Phys., 1998, 109, 6237–6245.
65. J. Pu and D. G. Truhlar, J. Chem. Phys., 2002, 117, 10675–10687.

---

Footnote

† Electronic supplementary information (ESI) available: Rate constants, energies of the benchmark, coordinates of stationary points. See DOI: 10.1039/c6fd00096g

---