Development and Validation of a Parameter-Free Model Chemistry for the Computation of Reliable Reaction Rates

A recently developed model chemistry (jun-Cheap) has been slightly modified and proposed as an effective, reliable, and parameter-free scheme for the computation of accurate reaction rates with special reference to astrochemical and atmospheric processes. Benchmarks with different sets of state-of-the-art energy barriers spanning a wide range of values show that, in the absence of strong multireference contributions, the proposed model outperforms the most well-known model chemistries, reaching a subchemical accuracy without any empirical parameter and with affordable computer times. Some test cases show that geometries, energy barriers, zero point energies, and thermal contributions computed at this level can be used in the framework of the master equation approach based on the ab initio transition-state theory for obtaining accurate reaction rates.


Introduction
For many years, scientists were skeptical about the presence of molecular systems in the interstellar space due to the harsh physical conditions (low temperature and pressure in the presence of high-energy radiation fields) characterizing this environment. However, contrary to those expectations, more than 200 molecules have been now identified in the interstellar and circumstellar medium (ISM), 1 including several so-called complex organic interstellar molecules (iCOMs), namely molecules containing carbon and a total of more than 6 atoms. 2 Most of the observed species should have a very short life-time according to earth-based standards, but the inter-molecular processes leading to thermodynamic equilibrium are not effective in the ISM due to its extreme physical parameters. 3,4 This situation calls for a strong interplay between observations, laboratory studies and computational approaches to understand the chemical evolution in these regions and to explain the observed abundances of different species.
Astrochemical models are virtual laboratories including thousands of reactions and whose main goal is to reproduce the observational data to the best possible extent. Although the available astrochemical models show widely different degrees of sophistication, 5 all of them share the same basic ingredients: 6 a set of initial conditions (total density, temperature, etc.) and a panel of chemical reactions characterized by their respective temperature-dependent rate constants and most likely exit channels. In order to improve the current predictions provided by these models, the reactions responsible for the largest uncertainties on the abundances must be studied in more detail by laboratory experiments and/or theoretical methods to provide improved rate constants and branching ratios.
Chemical kinetics plays a fundamental role also in the different, but related context of atmospheric models that try to reproduce and interpret the large number of chemical processes occurring in the troposphere. Reaction rate coefficients and product yields have been traditionally obtained either by means of suitable experimental techniques 7 or estimated using structure-activity relationships. 8 The massive number of organic compounds released in the atmosphere and the corresponding huge number of possible reactions ruling their oxidation/degradation pathways make experimental measurement of even a small fraction of key processes a daunting task. In recent years, computational chemistry has begun to contribute substantially to a better understanding of several important reaction sequences in the atmosphere. 9 These contributions have, at their heart, the use of electronic structure calculations to determine the energies and other characteristics (mainly geometries and vibrational frequencies) of stable species, reactive complexes and transition states, which are then used in theoretical frameworks to determine rate coefficients. The main factor limiting the accuracy of this process is the computation of accurate values for all the energy barriers ruling the different elementary steps. Next, zero point energies and finite temperature contributions come into play, whose contributions can become non-negligible already for medium size systems.
Several non-empirical procedures have been developed and employed for the generation of accurate thermochemical data, which for small systems come close to the full configuration interaction (FCI) complete basis set (CBS) limit. 10 Among the most successful approaches there are the Weizmann-n series (with the most accurate being W4 11 ), the focal point analysis (FPA) approach, 12,13 the Feller-Dixon-Peterson model (FDP) 14 and the extrapolated abinitio thermochemistry (HEAT) protocol. [15][16][17] A simplified version of the HEAT protocol is obtained by retaining only the extrapolation to the CBS limit at the CCSD(T) level and the incorporation of the core-valence corrections, thus leading to the model referred to in the following as CBS-CV. This approach is rather well tested in the literature and was shown to provide results with an accuracy well within 0.5 kcal mol −1 . Recently, alternative protocols have been proposed, which employ explicitly-correlated approaches: 10,18 thanks to the faster convergence to the complete basis set limit, these approaches allow some computer time saving, but the rate determining step remains the evaluation of higher level contributions.
For larger molecular systems, more approximate composite methods are unavoidable, which aim at reaching the so-called chemical accuracy (1 kcal mol −1 ). The most well known among these so-called model chemistries are the last versions of the Gn 19 (G4 20 ) and CBS 21 (CBS-QB3 22 ) families. However, all these models include some empirical parameters and employ geometries, which are not fully reliable for transition states and non-covalent complexes ruling the entrance channels of most reactions of astrochemical and atmospheric interest. As a matter of fact, the most reliable protocols (e.g. HEAT) push geometry optimizations to the limit in order to obtain accurate energetics, whereas, at the other extreme, Gn and CBS-x schemes employ B3LYP geometries, whose accuracy is often unsatisfactory. 23 In the last few years, a reliable and accurate computational protocol, referred to as cheap scheme (ChS) and devoid of any empirical parameter, has been developed and tested with remarkable success for structural and energetic data. [24][25][26] In conjunction with geometries and harmonic frequencies issuing from double hybrid functionals, ChS has given promising results also for the activation energies of some reactions of astrochemical interest. [27][28][29][30][31] More recently, an improved variant (referred to as jun-Cheap scheme, jChS) has been introduced, which, thanks to the use of the 'june' partially augmented basis set of the 'calendar' family, 32 provides very accurate results also for non-covalent interactions. 33,34 On these grounds, in this paper we provide a comprehensive benchmark of the jChS model chemistry for several classes of reactions for which accurate reference results are available or have been purposely computed. Together with electronic energies, we analyze also zero point energies, thermal contributions to enthalpies and entropies and overall reaction rates computed for elementary reactions in the framework of the master equation (ME) approach based on ab initio transition state theory (AITSTME). [35][36][37] The paper is organized as follows. In the first part we validate the jChS model chemistry with reference to some well-known databases, namely i) the 24 Next, the reliability of the jChS model chemistry for zero point energies and thermal contributions to enthalpies and entropies is assessed with respect to the new databases THCS21   and THOS10 containing accurate reference values for closed-and open-shell systems, respec-tively.
Finally, the role of different contributions in determining the overall accuracy of computed reaction rates is analyzed by means of some simple elementary reactions and two more complex reaction networks relevant for astrochemistry and atmospheric chemistry. Conclusions and perspectives are given in the last section.

Computational details
All the composite schemes employed in the present work extrapolate single point energies computed at suitable geometries (see next sections) using the cc-pV(n+d)Z (hereafter nZ) 43 or jun-cc-pV(n+d)Z (hereafter jnZ) 32 families of basis sets. The coupled cluster (CC) ansatz including single, double and (perturbative) triple excitations (CCSD(T)) 44 within the frozen-core approximation and in conjunction with 3Z or j3Z basis sets is always employed in the first step. Next, CBS extrapolation and core-valence correlation (CV) are added using either MP2 45 (leading, in conjunction with jnZ basis sets, to our standard jChS model) or CCSD(T). In the latter case, inclusion of higher-level terms (diagonal Born-Oppenheimer, 46-49 scalar relativistic, 50,51 full triple and perturbative quadruple excitations [52][53][54] ) and systematic use of nZ basis sets leads to the CBS-CVH scheme.
The effect of spin-orbit coupling is added to the energies of the O, OH, SH and Cl radicals, where the CBS term is and the core valence correction ∆E CV is the MP2 energy difference between all electron (ae) and frozen core (fc) calculations employing the cc-pwCVTZ basis set. 64 At this level, the extrapolation of Hartree-Fock (HF) and correlation contributions is performed with the same equation and basis sets since several tests have shown that this simplified recipe has a negligible impact on the overall accuracy of the results. Furthermore, scalar relativistic effects are neglected, which is not a serious approximation since the heaviest element involved in this study is Cl.

The CBS-CVH composite scheme
The CBS-CVH total electronic energies are obtained from single-point computations at geometries optimized by the jChS composite method described above for energies: In this case, HF and correlation energies are extrapolated separately. In particular, the HF CBS limit is estimated by using Feller's exponential formula 65 whereas the CBS limit of the correlation energy is obtained by the n −3 formula proposed by Helgaker and coworkers 66 The three-point extrapolation of HF energies employs 3Z,4Z and 5Z basis sets, whereas the two smaller basis sets are used in the two-point extrapolation of correlation energies.
The core valence correction ∆E CV is computed as the CCSD(T) energy difference between all electron and frozen core calculations employing the cc-pCVTZ basis set. 64 The diagonal Born-Oppenheimer correction ∆E DBOC 46-49 and the scalar relativistic contribution to the energy ∆E rel 50,51 are computed at the HF-SCF/aug-cc-pVDZ and CCSD(T)/augcc-pCVDZ level, after having checked their convergence with respect to contributions calculated with triple-zeta basis sets for a few stationary points.
Finally, the corrections due to full treatment of triple (∆E fT ) and perturbative treatment of quadruple (∆E pQ ) excitations are computed, within the fc approximation, as energy differences between CCSDT and CCSD(T) and between CCSDT(Q) and CCSDT calculations employing the cc-pVTZ and cc-pVDZ basis set, respectively.

Kinetic models
Global and channel-specific rate constants were computed solving the multi-well one-dimensional master equation using the chemically significant eigenvalues (CSEs) method within the Rice-Ramsperger-Kassel-Marcus (RRKM) approximation. 67 The collisional energy transfer probability is described using the exponential down model 68 with a temperature dependent ∆E down of 260 × (T /298) 0.875 cm −1 in an argon bath gas.
For channels ruled by a distinct saddle point, rate coefficients are determined by conventional transition state theory (TST) within the rigid-rotor harmonic-oscillator (RRHO) approximation 69 and including tunneling as well as non classical reflection effects by using the Eckart model. 70 Instead, rate constants for barrierless elementary reactions are computed employing phase space theory (PST), 71,72 again within the RRHO approximation.
The isotropic attractive potential V eff entering the PST is described by a C R 6 power law, whose C coefficient is obtained by fitting rev-DSD energies computed at various long-range distances of fragments. We obtained the following C coefficients for the PST calculations of barrierless channels: 230 a 0 6 E h for the H 2 S + Cl entrance channel, 64.2 a 0 6 E h for the CH 3 NH 2 + CN entrance channel on the methyl side and 94.4 a 0 6 E h for the CH 3 NH 2 + CN entrance channel on the nitrogen side.
The rate constants of the overall reactions evaluated in different temperature ranges are fitted by the three-parameter modified Arrhenius equation proposed by Kooij 73,74 : where A, n, and E are the fitting parameters, and R is the universal gas constant.

Results and discussion
In the original jChS model, geometries and force fields were computed with the B2PLYP double hybrid functional 75 augmented by empirical dispersion contributions (namely the D3(BJ) model) 76,77 in conjunction with partially augmented triple-zeta basis sets. 33 However, the recently developed rev-DSD model 56 delivers improved descriptions of non-covalent interactions and activation energies. 78,79 Therefore, we benchmarked the performances of this functional (still in conjunction with partially augmented triple-zeta basis sets) for geometrical parameters and vibrational frequencies obtaining results close to those delivered by the CCSD(T) ansatz in conjunction with comparable basis sets, but at much reduced computational cost. 80 As a consequence, the jChS model chemistry now uses by default rev-DSD geometries and force-fields.
If the spin contamination from higher spin states is large, the potential energy surfaces so that calculations at the CCSD(T) level are usually relatively insensitive to the choice of (restricted or unrestricted) orbitals. 84 However, in cases where higher spin contaminants become important, CCSD(T) can also fail. 83 On these grounds, all the jChS and CBS-CVH energies have been computed by the restricted open-shell approach.
Concerning DFT methods, it is well known that the extent of spin contamination in unrestricted versions of hybrid density functionals increases with the amount of HF exchange. 85 However, Menon and Radom 86 showed that in unrestricted double-hybrid procedures, the opposing behavior of UHF and UMP2 with respect to spin contamination leads to smaller differences between the energies predicted by unrestricted and restricted open-shell variants.
Although rev-DSD energies are not used in the present context, spin contamination can have an effect also on gradients and Hessians. We have, therefore, checked systematically the spin contamination and found that its effect is always negligible (within the target accuracy of the jChS model) except for the CN radical and the transition state ruling the reaction H • + F 2 −−→ HF + F • , which will be analyzed in detail in a following section.

Reaction barriers
The most well-known database of accurate reaction barriers is the DBH24 compilation 38,87 containing results mostly obtained at the CCSDTQ5/CBS level via W4 theory 88 for a statistically representative set including 3 prototypes for each of the following classes of reactions: heavy atom transfer, nucleophilic substitution, unimolecular and association reactions, and hydrogen-transfer reactions. Table 1 Table 1, this improvement is close to that obtained when going from CCSD(T)/j3Z (0.71 and 2.49 kcal mol −1 ) to jChS (0.36 and 0.80 kcal mol −1 ).
These trends suggest that either inclusion of explicit correlation or two-point extrapolation at the MP2 level are effective routes for improving significantly the accuracy of computed energy barriers, without introducing additional computational bottlenecks with respect to the underlying CCSD(T)/j3Z reference. As a matter of fact, already for reactions involving two heavy atoms (e.g., A7, A8, A9, A10 in Table 1  Two larger databases of prototypical reactions are also available for barriers related to transfers of hydrogen and non-hydrogen atoms (HTBH38 39 and NHTBH38, 40 respectively).
However, the reaction barriers not already included in the DBH24 set have been obtained at lower computational level (W1). We have thus decided to compute at the jChS level all the reactions of the above two sets not contained in the original DBH24 compilation using both rev-DSD and the original QCISD/MG3 geometries. Whenever significant discrepancies were found, the reactions were recomputed also at the CBS-CVH level.
The reactions from the NHTBH38 set not included in the DBH24 selection are collected in Table 2. It is noteworthy that rev-DSD energy barriers, although not directly used in the jChS model chemistry, show MUEs smaller than 2.0 kcal mol −1 , thus suggesting that the corresponding geometries should be sufficiently accurate for single-point energy evaluations at higher computational levels. This is confirmed by the finding that only for reaction NHT3, QCISD and rev-DSD geometries lead to significantly different results (cfr. columns 2 and 3 of Table 2). Geometry optimization at the jChS level provides results far  The reactions from the HTBH38 set not included in the DBH24 selection are collected in Table 3. Once again it is noteworthy that the rev-DSD energy barriers, although not directly used in the jChS model chemistry, do not show any unrealistic outlier. Only for reactions HT1 and HT5 QCISD and rev-DSD geometries lead to significantly different results (cfr. columns 2 and 3 of Table 3). Geometry optimization at the jChS level provides results close We have then selected six 'challenging reactions' among those mentioned above for further investigation. To this end, we report in Table 4    In order to extend the benchmark to larger and more complex systems we resorted to the BH28 set of ref. 41, which includes accurate (W3lite-F12) energy barriers for several peri-   Figure 3 and the corresponding forward and reverse reaction barriers (BH14 set) are collected in Table 5. The average and maximum errors are larger than those of the DBH24 set, but, closer inspection of the results shows that, as already pointed out in ref. 41, the role of full triple turbation theory in conjunction with a separate treatment of large amplitude motions. 58,98 In fact, a resonance-free expression for ZPEs of energy minima and transition states, 99,100 an unsupervised smoothing procedure (HDCPT2) for fundamental frequencies 101 and a fully automatic detection and treatment of torsional motions (hindered rotor, HR, approximation) 102 have been implemented in the Gaussian code 60 and validated. 59 As a consequence a fully black-box procedure is available for taking into account all these contributions.
Next, the so-called simple perturbation theory (SPT) 103 can be applied for computing partition functions without the need of performing explicit (or stochastic) summations of individual energy levels. In fact, the SPT retains the formal expression of the harmonic partition function, but employing the anharmonic ZPE and fundamental levels (∆ i ) issuing from HDCPT2 and HR computations.
This approximation provides results in remarkable agreement with accurate reference values and leads to analytical expressions for the different thermodynamic functions. 103 On these grounds, we will now analyze the performances of the jChS model chemistry in dealing with these terms starting from a benchmark of the RRHO approximation with reference to accurate quantum chemical results and then proceeding to take into account anharmonic contributions. For illustration purposes, we will focus our attention on ZPEs and absolute entropies (S), which are especially sensitive to high and low-frequencies, respectively.
To this end a new database has been built (ThCS21), which contains accurate experimental values for the ZPEs and absolute entropies of 21 semi-rigid closed-shell molecules, whose estimated errors are below 0.1 kcal mol −1 and 0.05 cal/(mol K), respectively. The results collected in Table 6 show that already at the harmonic level, the errors are well within the level of accuracy expected from the jChS model chemistry and the anharmonic results can be confidently used in conjunction with the most sophisticated models (e.g. CBS-CVH). Actually, the harmonic frequencies obtained at this level do not require any empirical correction to compensate for method and/or basis set deficiency, but only for genuine anharmonic effects, which, in turn, give significant contributions to ZPEs only for some XH bonds (X=C,N,O).
As a consequence, an empirical correction of 0.12 kcal mol −1 for each bond of this kind provides results very close to the anharmonic counterparts (see results in parenthesis in the first column of Table 6).  Accurate entropy values are also available for the same set of molecules and harmonic computations perform a remarkable job in reproducing the experimental values. However, entropy is exquisitely sensitive to low-frequency vibrations, so that a set of flexible molecules is collected in Table 7. It is apparent that the HRHO model (which does not add any computational burden with respect to the underlying RRHO model) performs a remarkable job for systems containing a single torsion. The situation is more involved for larger flexible systems due to the presence of several low-energy minima contributing to the overall thermodynamic functions. Although this aspect goes beyond the main topic of the present contribution, we point out that several strategies are being proposed, following systematic search, 107 stochastic 108 and, more recently, machine learning 109 approaches. Other kinds of large amplitude motions can be taken into account by means of one-dimensional variational or quasi-variational approaches 110 followed by SPT or direct count of energy levels. 98   Table 8 for a few representative systems (ThOS10 database) suggest that (in the absence of strong multireference effects) the expected accuracy is close to that reached for closed-shell systems.

Reaction rates
In this section we analyze the impact on reaction rates of the different ingredients discussed in the previous section, comparing the results issuing from different model chemistries including CBS-QB3, jChS and CBS-CVH. Starting from simple elementary mechanisms we proceed to more complex potential energy surfaces including several intermediates and transition states, possibly leading to different products.
The first test case is the high pressure limit of the reaction H • + CO, which has been recently investigated by Vichietti et al. 116 This reaction belongs to the HTBH38 set, whose jChS results have been discussed in the Section devoted to energy barriers. though the presence of a van der Waals pre-reactive complex has been suggested, its stability (if any) is so small that its impact on the computed reaction rates is negligible.
The reaction rates computed in the 50-4000 K temperature interval are shown in Figure 4 and the parameters of the corresponding Arrhenius-Kooij fits obtained by different electronic structure methods are collected in Table 9. The non-Arrhenius behaviour of the reaction is approach at least at low temperatures.   We next consider the BHPERI1 and CRBH4 reactions discussed in the section on the energy barriers (see Figure 3a and 3c). The rates computed in the 300-1000 K temperature interval by different electronic structure methods are shown in Figure 5a and 5b, whereas the parameters of the corresponding Arrhenius-Kooij fits are collected in Table 10. Both reactions are characterized by quite high energy barriers and their rates show a clear Arrhenius behaviour. In these circumstances the different electronic structure methods deliver comparable results over the whole temperature range.  The next example is the reactive potential energy surface for H 2 S + Cl (see Figure 6), which involves a van der Waals pre-reactive complex (RW) followed by the transition state TS leading to a product-like van der Waals complex (PW) and then to the products, i.e., HS + HCl. Since this reaction has been recently investigated at the CBS-CVH level, 28  The reaction rates issuing from jChS computations are compared in Figure 7 to the CBS-QB3 and CBS-CVH counterparts, whereas the parameters of the corresponding Arrhenius-Kooij fittings (see Equation 6) are collected in Table 11. The root mean square deviations reported in Table 11 demonstrate that the data are indeed well fitted by the Arrhenius-Kooij expression with a negative activation energy (E) at 0 K. The results issuing from jChS and CBS-CVH computations are virtually indistinguishable, whereas significantly larger rates are obtained at the CBS-QB3 level.   The last example is the quite complex reactive potential energy surface ruling the addition of CN to CH 3 NH 2 shown in Figure 8 together with the jChS energies of all the stationary points. The experimental reaction rate at different temperatures 117 has been recently well reproduced employing CBS-CVH electronic structure computations within a master equation treatment similar to that employed in the present paper. 118 This system represents, therefore, a challenging test for the jChS model.
The attack of CN on the nitrogen side of methylamine proceeds via a potential well associated with a pre-reactive complex, NC···NH 2 CH 3 (IC), which evolves in an inner (submerged) transition state (TS3) that, passing through an NCH···NHCH 3 intermediate (FC02), forms the HCN + NHCH 3 products (P3). Alternative channels, and in particular that leading to NH 2 CN + CH 3 , are ruled by non-submerged transition states and are, therefore, closed under the ISM conditions. The attack on the methyl side forms directly the FC01 complex, which, in turn, leads to HCN + NH 2 CH 2 (P1) without any potential energy barrier. In this case, the alternative two-step mechanism (TS0-RI-TS2-P2) leading to aminoacetonitrile + H is open since it involves only submerged transition states, but it is less favorable. The errors of the CBS-QB3 22 model are again larger than 1 kcal mol −1 , in agreement with the estimates of previous studies. 119 The reaction rates issuing from jChS computations are compared in Figure 9 with the CBS-QB3 and CBS-CVH counterparts, whereas the parameters of the corresponding Arrhenius-Kooij fits (see Equation 6) are collected in Table   12. Noted is that pressure does not influence the reaction rate, as the reactants always proceed to form the products without experiencing significant collisional stabilization in the investigated pressure range (0.001-1 bar).  A curved Arrhenius plot is obtained when the activation energy depends on the temper-ature and this behavior is captured by the Arrhenius-Kooij equation when this dependence is linear. The root mean square deviations reported in Table 12 demonstrate that the data for the different products are indeed well fitted by the Arrhenius-Kooij expression. Within this model, E represents the activation energy at 0 K and the activation energy at a generic temperature T is given by E + n RT 300 . In the present case, the activation energy is positive for P1 and P3, as a result of both the capture rate constant and the subsequent energy barriers for the unimolecular steps. The value is, instead negative for P2, but in this case the Arrhenius plot is essentially linear. The n parameter (the first derivative of the activation energy with respect to temperature) is positive for the P1 and P2 products, thus reflecting an increase of the activation energy with temperature, while the opposite behavior (n negative) is obtained for P3. Finally, the value of the pre-exponential factor A is typical for this kind of reactions and rules the branching ratio between the different channels.

Supporting Information Available
The Supporting Information is available free of charge at xxxxxxxx. Extended Table of re-sults for the HTBH38/08 dataset, cartesian coordinates of all the stationary points optimized in this work and MESS input files.