香京julia种子在线播放

The prerequisite of any accurate kinetic study is the characterization of all the stationary points (minima and transition states) of the reactive potential energy surface (PES) ruling the process under investigation. This step is conveniently performed using hybrid or double-hybrid density functionals that incorporate empirical dispersion corrections by means of the so-called D3 model combined with Becke-Johnson damping [D3BJ; (Grimme et al., 2010; Grimme, 2011)]. A preliminary investigation employing suitable hybrid density functionals [in our approach, B3LYP (Lee et al., 1988; Becke, 1993) or PW6B95 (Zhao and Truhlar, 2005)] in conjunction with (partially) augmented double-zeta basis sets [preferably jul-cc-pVDZ (Papajak et al., 2011)] is usually sufficient. Next, both geometries and energies are refined at a higher level. Concerning geometries, double-hybrid functionals in conjunction with (partially) augmented triple-zeta basis sets usually represent a remarkable compromise between accuracy and computational cost. In this connection, the B2PLYP-D3/maug-cc-pVTZ model (Grimme, 2006; Papajak et al., 2009) employed in our earlier studies has recently been replaced by the rev-DSD-PBEP86 functional (Santra et al., 2019) in conjunction with the jun-cc-pVTZ (Papajak et al., 2011) basis set (such a functional-basis combination is shortly denoted as rev-DSD in the following), which usually delivers improved results (with respect to the former) for transition states and non-covalent complexes without reducing the accuracy of other parameters. When more accurate geometries are sought (especially in connection with spectroscopic properties or inter-molecular distances of pre-reactive complexes), more accurate quantum-chemical approaches (vide infra) or a “LEGO bricks” approach can be employed (Ceselin et al., 2021; Melli et al., 2021).

The methodology presented in this work focuses on the energetics and implies that geometries are optimized using a double-hybrid density functional as explained above. The first model chemistry we address is the so-called ‘jun-Cheap’ scheme (Alessandrini et al., 2020; Lupi et al., 2021), which employs rev-DSD optimized geometries. While interested readers are referred to the specific literature for a complete account, here we briefly mention that it is based on the coupled-cluster (CC) singles, doubles and perturbative triples method [CCSD(T); (Raghavachari et al., 1989)] in conjunction with the jun-cc-pVTZ basis set [this implies that jun-cc-pV(T+d)Z is used for atoms of third period; (Papajak et al., 2011)]. To improve the fc-CCSD(T)/jun-cc-pVTZ level of theory, with fc standing for the frozen-core approximation, the extrapolation to the complete basis set limit (CBS) is performed and the core-valence correlation (CV) correction added. To reduce the computational cost of these two latter contributions, they are evaluated using Møller-Plesset perturbation theory to the second order [MP2; (Møller and Plesset, 1934)]. Within any ‘cheap’ scheme, as well tested in the literature (Lupi et al., 2020b; Alessandrini et al., 2021; Lupi et al., 2021), the MP2 energy is extrapolated using the n ⁻³ formula proposed by Helgaker and coworkers (Helgaker et al., 1997). Overall, the jun-Cheap energy is given by the following expression: E jun - Cheap = E ( C C S D ( T ) / T Z ) + Δ E M P 2 / C B S + Δ E M P 2 / C V , (1) where the first term on the right-hand side is the energy at the fc-CCSD(T)/jun-cc-pVTZ level. The second term is the contribution due to the MP2 extrapolation to the CBS limit, obtained –within the fc approximation– as follows: Δ E M P 2 / C B S = n 3 E MP 2 / n Z − n − 1 3 E MP 2 / n − 1 Z n 3 − n − 1 3 − E MP 2 / n − 1 Z , (2) where n = 4 (i.e., jun-cc-pVQZ and n-1 stands for jun-cc-pVTZ). Finally, the ΔE(MP2/CV) term is computed as energy difference between MP2 calculations correlating all electrons and within the fc approximation, both in the cc-p(w)CVTZ basis set (Woon and Dunning Jr., 1995; Peterson and Dunning Jr., 2002).

An alternative model (hereafter referred to as the jun-F12 scheme) replaces the conventional CCSD(T) method with its explicitly-correlated CCSD(T)-F12 counterpart (Knizia et al., 2009), still in conjunction with the jun-cc-pVTZ basis set. While the extrapolation to the CBS limit is avoided (because of the fast convergence of F12 methods), the CV correction is incorporated at the MP2-F12 level (Werner et al., 2007) employing the cc-p(w)CVTZ basis set (Woon and Dunning Jr., 1995; Peterson and Dunning Jr., 2002). In all CCSD(T)-F12 computations, the F12b approximation (Feller and Peterson, 2013) is employed, while the default approximation (3C with EBC) is used for MP2-F12. It has to be noted that a full F12 treatment of triple excitations is not yet available in commercial codes; thus, their contribution is essentially the same as in conventional CCSD(T).

To benchmark the performance of the two model chemistries introduced above, reference electronic energies are obtained by means of the so-called HEAT-like (HL) scheme, closely resembling the original HEAT methodology proposed by Stanton and co-workers (Tajti et al., 2004; Bomble et al., 2006; Harding et al., 2008; Lupi et al., 2020b), and entirely based on the CC ansatz. In this benchmark, the HL composite scheme is built on top of geometries optimized at the jun-Cheap level: E HL = E HF CBS + Δ E CCSD ( T ) CBS + Δ E CV + Δ E fT + Δ E pQ + Δ E rel + Δ E DBOC (3)

As evident from the formula above, HF and CCSD(T) correlation energies are extrapolated separately (differently from the jun-Cheap scheme). In particular, the HF CBS limit is estimated by using Feller exponential formula (Feller, 1992): E HF n = E HF CBS + B ⁡ exp − C n (4)

This is a three-point extrapolation expression used in conjunction with triple- quadruple-, and quintuple-zeta basis sets. Instead, the CBS limit of the CCSD(T) correlation energy is obtained using the two-point n ⁻³ formula (see above, the discussion made for the jun-Cheap scheme) in combination with triple- and quadruple-zeta basis sets. The CV term is computed at the CCSD(T) level as energy difference, as described above. Higher-order terms in the CC expansion are also incorporated. These are the corrections due to the full treatment of triples, ΔE _fT, and the perturbative treatment of quadruples, ΔE _pQ, which are computed, within the fc approximation, as energy differences between CCSDT (Noga and Bartlett, 1987; Scuseria and Schaefer, 1988; Watts and Bartlett, 1990) and CCSD(T) and between CCSDT(Q) (Bomble et al., 2005; Kállay and Gauss, 2008) and CCSDT calculations employing the cc-pVTZ and cc-pVDZ basis sets, respectively. The last two terms are the diagonal Born-Oppenheimer correction, ΔE _DBOC (Sellers and Pulay, 1984; Handy et al., 1986; Handy and Lee, 1996; Kutzelnigg, 1997), and the scalar relativistic contribution, ΔE _rel (Cowan and Griffin, 1976; Martin, 1983), which are computed at the HF-SCF/aug-cc-pVDZ and CCSD(T)/aug-cc-pCVDZ levels, respectively. If required, the effect of spin-orbit coupling is added to the electronic energies, mainly relying on the available experimental data (Moore, 1949).

A term not included in Eq. 3 is the zero-point energy (ZPE) correction. Within the harmonic approximation, this requires the evaluation of harmonic vibrational frequencies, which are computed at the same level of theory employed in the structural determinations. Indeed, this type of calculations is needed to check the nature of the stationary point under consideration, and in particular, for discriminating between minima and transition states. In this work, harmonic force-filed computations have been carried out at the rev-DSD level exploiting analytical second derivatives (Biczysko et al., 2010). The ZPE term can be improved by incorporating anharmonic corrections, which can be computed using hybrid or double-hybrid functionals within the generalized second-order vibrational perturbation theory [GVPT2; (Barone, 2005; Barone et al., 2014; Frisch et al., 2016)].

As far as technical details are concerned, all DFT, MP2 and single-point energy CCSD(T) computations are performed using the Gaussian program package (Frisch et al., 2016), while the CFOUR program (Stanton et al., 2016) is employed whenever geometry optimizations at the CCSD(T) level are needed as well as for relativistic and diagonal Born-Oppenheimer corrections. For explicitly-correlated calculations, the MOLPRO program (Werner et al., 2020) is used. Finally, CCSDT and CCSDT(Q) energy evaluations require the MRCC program (Kállay et al., 2018).

Global and channel-specific rate constants are computed solving the multi-well one-dimensional master equation using the chemically significant eigenvalues (CSEs) method within the Rice-Ramsperger-Kassel-Marcus (RRKM) approximation (Georgievskii et al., 2013). The collisional energy transfer probability is described using the exponential down model (Tardy and Rabinovitch, 1966) with a temperature dependent ΔE _down of 260 × (T/298)^0.875 cm⁻¹ in an argon gas bath. The ratio between the specific rate constant for each reaction product and the overall rate constant accounting for all products leads to the reaction branching ratios.

For channels ruled by a distinct saddle point, rate coefficients are determined using conventional transition state theory (TST) within the rigid-rotor harmonic-oscillator (RRHO) approximation (Fernández-Ramos et al., 2006), also incorporating tunneling and non-classical reflection effects by means of the Eckart model (Eckart, 1930). Instead, rate constants for barrierless elementary reactions (usually, the entrance channels) are evaluated employing the phase space theory [PST; (Hunter et al., 1993; Skouteris et al., 2018)] within the RRHO approximation. The isotropic attractive potential V _eff entering the PST is described by a C R 6 power law, with the C coefficient obtained by fitting rev-DSD energies computed at various long-range distances of fragments.

To model the temperature dependence, the rate constants are evaluated at different temperatures and fitted using the three-parameter modified Arrhenius equation proposed by Kooij (Kooij, 1893; Laidler, 1996): k T = A T 300 n ⁡ exp − E a R T (5) where A, n, and E _a are the fitting parameters, and R is the universal gas constant.

All kinetic computations presented in this manuscript have been performed with the MESS code (Georgievskii et al., 2013).

As a first benchmark, we focus our attention on the electronic energy contributions to reaction and activation energies employing as references the very accurate results of the well-known DBH24 database (Karton et al., 2008; Zheng et al., 2009), which contains results from state-of-the-art quantum-chemical approaches (estimated accuracy: 0.1–0.2 kJ/mol) for a statistically representative set of heavy atom transfer (NH), nucleophilic substitution (NS), unimolecular and association (UA) and hydrogen transfer (HT) reactions. Table 1 collects the reaction energies obtained at different levels. It is quite apparent that both the jun-F12 and jun-Cheap schemes fulfill the target of sub-chemical accuracy without the need of any empirical parameter. Furthermore, owing to the inclusion of the CV correlation contribution, the accuracy of the results for molecules containing third-row atoms is comparable to that of the first- and second-row-atom counterparts.

From the inspection of Table 1 and, in particular, from the comparison of the first (fc-CCSD(T)/jun-cc-pVTZ) and the third (jun-Cheap) column, it is noted that the rather inexpensive MP2 estimate of the CBS and CV contributions reduces to about one third the deviation from the reference values. Moving from the jun-Cheap scheme to the jun-F12 model leads to a lowering of the mean unsigned error (MUE), while the maximum unsigned error (MAX) remains nearly unchanged. These results indicate that either the inclusion of an explicit correlation treatment or the extrapolation to the CBS limit at the MP2 level for conventional methods offer reliable routes for obtaining accurate reaction energies without any significant computational increase with respect to the underlying CCSD(T) step. In particular, already for reactions involving two non-hydrogen atoms, jun-Cheap computations require no more than twice the computer time of the CCSD(T) step. The effectiveness of the jun-Cheap model steeply increases with the dimension of the target system because of the favorable scaling of MP2, which can be further enhanced by the resolution of identity and other acceleration techniques. In the same vein, inclusion of explicit correlation through the F12 ansatz increases by 20–30% the computer time with respect to the standard CCSD(T) model employing the same basis set.

In Table 2, the energy barriers for the same DBH24 set as above, obtained at different levels of theory, are reported. It is remarkable that the errors associated to both the jun-F12 and jun-Cheap approaches only slightly increase with respect to those of the corresponding reaction energies, with MUEs of 1.2 and 1.5 kJ/mol and MAXs of 3.8 and 3.3 kJ/mol, respectively. This outcome paves the route toward the systematic use of the jun-F12 and jun-Cheap schemes in place of the widely used CCSD(T)/cc-pVTZ model with a remarkable increase of the accuracy, but without any significant increase of the computational cost. Furthermore, it has to be noted that the results are much improved with respect to those delivered by well-known (and largely employed) parametrized model chemistries such as CBS-QB3 (Montgomery et al., 1999) and G4 (Curtiss et al., 2007), whose MUEs are around 2.5 kJ/mol (Zheng et al., 2009). To put things in a wider perspective, we mention that the MUEs and MAXs of jun-F12 and jun-Cheap are close to those obtained employing larger basis sets, which in turn imply computational costs that might not be affordable for large systems. In detail, the MUE and MAX for CCSD(T)/aug-cc-pVQZ and CCSD(T)-F12/aug-cc-pVQZ are 1.2 and 4.1 kJ/mol and 1.1 and 2.5 kJ/mol, respectively (Zhang and Valeev, 2012).

Together with the accuracy of the electronic energies, also the quality of the geometrical structures employed in their evaluation can play a role in the final results (Puzzarini and Barone, 2018; Lupi et al., 2020b). The numerical evidence of this aspect has been investigated for the specific case of a radical attack to neutral molecules considering three extraction (Ex) and three addition (Add) reactions. The results are collected in Table 3. A re-optimization of the rev-DSD geometries at the jun-Cheap level has been performed. This slightly affects the energetic results, with differences within 0.6 kJ/mol, as evident in the ΔGeom column of Table 3. Therefore, rev-DSD structures can be safely employed in the evaluation of the energetics ruling astrochemical reactions.

Another important issue to address is the effect of improved quantum-chemical treatments on the energetics and, for our discussion, we again resort to the “Ex” and “Add” reactions of Table 3. We note that the replacement of MP2 by CCSD(T) in both the extrapolation to the CBS limit and the evaluation of core-valence correlation (ΔCBS-CV) leads to very small discrepancies. The column ΔHL of Table 3 collects the differences arising from the incorporation of higher-order contributions in the coupled-cluster expansion (full treatment of triples and perturbative inclusion of quadruples) as well as the diagonal Born-Oppenheimer and relativistic corrections. As often noted in the literature (Puzzarini and Barone, 2020; Puzzarini et al., 2020), the sum of these contributions leads to a small corrective term, which is generally well within 1 kJ/mol. Overall, these results indicate that the jun-Cheap scheme is robust and very accurate, and that its improvement requires expensive calculations (and, often, it is not worth it). Furthermore, all the other contributions discussed above should be considered together because none of them plays a dominant role and all of them are rather small.

As already mentioned, electronic energies need to be corrected for the ZPE contributions in order to obtain the data to be used in the subsequent kinetic study. In Table 3 harmonic and anharmonic ZPE contributions to prototypical extraction and addition reactions are compared. As expected, these corrections are relevant, indeed ranging from a few kJ/mol to contributions as large as 19 kJ/mol. However it is noted that one can resort to harmonic ZPEs because the corresponding anharmonic corrections are well within 1 kJ/mol and can thus be safely neglected whenever –as in the present case– kJ/mol accuracy is sought.

However, further details on the vibrational frequencies determining ZPE corrections and also thermal effects (the latter via the vibrational partition functions) are deserved. At the harmonic level, several benchmark studies have shown that hybrid and, especially, double-hybrid functionals deliver quite reliable results, provided that empirical scaling factors (that should be different for ZPE and vibrational partition function) are employed (Kesharwani et al., 2015). Empirical parameters can be avoided by resorting to GVPT2 for small amplitude motions (Barone, 2005) together with a separate treatment of individual large amplitude motions, if any (Ayala and Schlegel, 1998). As an example of the accuracy reachable, a recent benchmark study for a database of 21 closed-shell molecules for which accurate reference values are available led to average errors of 0.2 kJ/mol and 0.2 cal/(mol K) for ZPE corrections and entropies, respectively (Barone et al., 2021). Experimental data are available only for a limited number of open-shell species, but the 10 small radicals investigated in Barone et al. (2021) suggest that a comparable accuracy can be obtained.

As a last comment, we address the impact of the level of theory on the computed reaction rates. To illustrate this, we consider the rate of the H^• + CO reaction in the 50–4000 K temperature range, evaluated employing the CBS-QB3, jun-Cheap, jun-F12 and HEAT-like models. The results are graphically shown in Figure 2. From the inspection of this figure, the non-Arrhenius behavior of the reaction is apparent and the limits of the CBS-QB3 model at low temperatures (i.e., those of interest for astrochemical studies) are evident. Instead, the jun-Cheap, jun-F12 and HEAT-like results are very similar to each other in the whole temperature range and are similar to such an extent that only the jun-Cheap data are depicted in Figure 2, the jun-F12 and HEAT-like curves being indistinguishable from the jun-Cheap one at the scale of the figure. Finally, it is interesting to note that they are also in quantitative agreement with the available experimental data (Vichietti et al., 2020).

The results discussed above give full support to the reliability of the jun-Cheap and jun-F12 computational models for studying reactions of astrochemical interest. Since the explicit correlation route (jun-F12) surely deserves further investigations, in the following we retain the jun-Cheap model. One of its advantages is that it can be accessed using a large panel of available computer codes.

Among the iCOMs, the species containing the -C ≡N moiety have a strong prebiotic character because they represent key intermediates toward the formation of biomolecule building blocks such as aminoacids and nucleobases. A large number of cyanides have been detected in the ISM (McGuire, 2018; Müller, 2021), and it has been demonstrated that they can be formed not only by reactions on dust-grain surfaces (Lovas et al., 2006; Krim et al., 2019), but also directly in the gas phase (Vazart et al., 2015; Tonolo et al., 2020; Lupi et al., 2020a). In the following, we focus on two specific cyano-compounds: acrylonitrile and cyanomethanimine. As it will be shown, they can be obtained by a general addition/elimination mechanism following the attack of the CN radical to the isoelectronic ethylene and methanimine molecules, respectively.

Figure 3 shows the stationary points ruling the addition of CN to ethylene. As far as the reaction mechanism is concerned, hydrogen abstraction could be in principle competitive with addition/elimination (Bowman et al., 2020). Both mechanisms start with the formation of the very stable intermediate 1 (Figure 3) and then proceed toward products through a number of steps. However, the transition state ruling hydrogen abstraction (leading to vinyl radical and HCN) lies above the reactants; therefore, this channel is closed in the ISM conditions and not shown in Figure 3. The only open channel leads to acrylonitrile (AN) and atomic hydrogen by either a one- or a two-step mechanism, with the former route being favored at low temperatures.

The temperature dependence of the resulting rate constant is shown in Figure 4. It is apparent that, in the 30–300 K range, the rate coefficient varies a little (from 4×10⁻¹⁰ to 5.9 ×10⁻¹⁰ cm³molec⁻¹ s⁻¹) and follows the expected Arrhenius behavior. Furthermore, different electronic structure methods provide comparable results, as seen in the figure for the jun-Cheap and CBS-QB3 models.

This reaction has been investigated with a variety of experimental techniques (Sims et al., 1993; Balucani et al., 2000; Leonori et al., 2012), which demonstrated that the reaction is very fast, approaching the gas-kinetics limit, also at very low temperatures and that the main reaction channel is the one leading to AN via an addition/elimination mechanism (Balucani et al., 2000). The computed reaction rate at low temperatures is in quantitative agreement with the experimental results by Sims et al. (1993). However, the experimental temperature dependence shows a rather surprising non-Arrhenius behavior, which leads, at room temperature, to a disagreement by a factor of 2 or more. While further investigations of this aspect is surely warranted [for instance, additional reaction channels should become open when increasing the temperature, as suggested by Leonori et al. (2012)], only the low-temperature region is of interest for astrochemical studies and, furthermore, a factor of 2 error on computed reaction rates is usually more than acceptable.

The situation is more involved for methanimine because of the presence of two possible attack sites: the C-atom and N-atom of CH₂NH. Analogously to the C₂H₄ + CN reaction, the contribution of the abstration reaction is negligible, with test computations indicating that the corresponding rate constant is at least one order of magnitude smaller than that of the addition-elimination route. The addition/elimination mechanism is depicted in Figure 5, with the relative energies being computed at the jun-Cheap level. The comparison of these results with those recently obtained using an HEAT-like approach (Puzzarini and Barone, 2020) confirms the accuracy of the jun-Cheap model. The intermediate obtained upon attack to the N-side, 1N, is slightly more stable than those from the C-side attack, i.e., 1Z and 1E. However, the reaction channels originating from 1N are ruled by transition states significantly less stable (albeit submerged) than those issuing from 1Z or 1E.

Concerning the attack at the C-side, starting from the very stable 1Z (or 1E) pre-reactive complex, direct loss of a hydrogen atom, leading directly to the Z (or E) isomer of C-cyanomethanimine (C-CMIM), involves an exit barrier of 163.5 (or 164.1) kJ/mol. On the other hand, the stabilizing CN moiety on the carbon atom can promote hydrogen migration, which can lead to localize the unpaired electron on this atom. This migration is ruled by the transition state TS-1Z2 (TS-1E2 for the E form of C-CMIM) lying 140.1 kJ/mol above 1Z (144.0 kJ/mol above 1E for the E-route) and leads to 2, the most stable intermediate of the whole PES. Next, loss of hydrogen again leads to the Z (or E) form of C-CMIM and atomic hydrogen (denoted as PZ and PE in Figure 5) through the submerged transition state TS-2PZ (TS-2PE) lying about 49.6 kJ/mol (47.4 kJ/mol for the E route) below the reactants. As mentioned above, attack of the CN radical at the N-side of methanimine produces the very stable 1N intermediate from which two different paths can be followed in analogy with what is found for the C-side attack. While 1N is more stable than both 1Z and 1E, all the other stationary points of the N-side route are less stable than the C-route counterparts. The final product, i.e., N-cyanomethanimine (N-CMIM) + H (denoted as PN in Figure 5), is less stable than the products of the C-side attacks (PE or PZ) by about 30 kJ/mol.

Moving from thermochemistry to kinetics, the comparison between the jun-Cheap and CBS-QB3 rate constants and their temperature dependence is depicted in Figure 6 for the N-attack and in Figure 7 for the C-attack. From the inspection of these two figures, it is evident that the rate constant for the channel leading to N-CMIM is one order of magnitude smaller than those for the formation of the Z/E forms of C-CMIM. We also note that the discrepancies between jun-Cheap and CBS-QB3 results are significant (about a factor of 3–4). Focusing on the C-side attack, the prevalence of the PZ product over PE is due to the slightly lower energy barriers ruling the former path, with back-dissociation to reactants being always negligible. The overall rate constant for the formation of cyanomethanimine shows an increase with the temperature, which is accompanied by an increasing deviation from the Arrhenius behavior. At temperatures below 300 K, the value of the overall rate constant is determined by the faster one-step channel ruled by the TS-1ZPZ (TS-1EPE) transition state, while at higher temperatures the two-step mechanism starts to give non-negligible contributions.

The computed branching ratio (0.42/0.58) for the E and Z isomers at 150 K [the estimated temperature of G+0.693 where both isomers have been recently detected (Rivilla et al., 2018)] is much lower than the value of the abundance ratio derived from astronomical observations (which is about 6). While the observed ratio is close to that computed from a thermodynamic estimate based on the relative stability of the E and Z isomers (5.5), a fully satisfactory explanation of this puzzle requires the investigation of all possible destruction pathways for CMIM isomers and/or the demonstration that the thermodynamic equilibrium can be effectively established between the two isomers, even at low temperatures. Concerning the first issue, previous studies showed that all dissociation processes (also considering the cyclization to azirin-2-imine) require too much energy to occur in the ISM (Vazart et al., 2015). As regards the second issue, a recent study (García de la Concepción et al., 2021a) demonstrated that a multidimensional small curvature treatment of quantum-tunneling effects provides quite large (and nearly constant) inter-conversion rates, these being in the order of 10^–11 sec⁻¹ at temperatures between 10 and 150 K. Above the latter temperature, the contribution of tunneling becomes negligible with respect to thermal activation and the rate follows an Arrhenius behavior. The abundance ratio obtained when accounting for tunneling (6.0) is actually close to the value suggested by observations and to that computed from the thermodynamic equilibrium [provided that accurate energetics is used; see Puzzarini and Barone (2020)]. The situation being different for deuterated species for which the contribution of tunneling is always negligible, their astronomical observations would be of great interest and great help.