State-of-the-Art Quantum Chemistry Meets Variable Reaction Coordinate Transition State Theory to Solve the Puzzling Case of the H2S + Cl System.

The atmospheric reaction of H2S with Cl has been reinvestigated to check if, as previously suggested, only explicit dynamical computations can lead to an accurate evaluation of the reaction rate because of strong recrossing effects and the breakdown of the variational extension of transition state theory. For this reason, the corresponding potential energy surface has been thoroughly investigated, thus leading to an accurate characterization of all stationary points, whose energetics has been computed at the state of the art. To this end, coupled-cluster theory including up to quadruple excitations has been employed, together with the extrapolation to the complete basis set limit and also incorporating core-valence correlation, spin-orbit, and scalar relativistic effects as well as diagonal Born-Oppenheimer corrections. This highly accurate composite scheme has also been paralleled by less expensive yet promising computational approaches. Moving to kinetics, variational transition state theory and its variable reaction coordinate extension for barrierless steps have been exploited, thus obtaining a reaction rate constant (8.16 × 10-11 cm3 molecule-1 s-1 at 300 K and 1 atm) in remarkable agreement with the experimental counterpart. Therefore, contrary to previous claims, there is no need to invoke any failure of the transition state theory, provided that sufficiently accurate quantum-chemical computations are performed. The investigation of the puzzling case of the H2S + Cl system allowed us to present a robust approach for disclosing the thermochemistry and kinetics of reactions of atmospheric and astrophysical interest.

Introduction

Transition state theory (TST) and its variational extension are accepted to be the “workhorses” in computational kinetics for several fields ranging from combustion to atmospheric chemistry and astrochemistry (see, e.g., refs (1−3)). It is thus fundamental to understand the features of a reactive potential energy surface (PES) that can introduce a breakdown of such theory. In other words, it is important to rationalize its possible limitations. In this respect, some cases of breakdown of TST have been reported in the literature. For instance, Hase and co-workers pointed out that many recrossings can take place in the Cl– + CH3Cl gas-phase reaction,[4] and the same behavior was observed for the unimolecular isomerizations of NCCN and CH3CN.[5] Another interesting example is offered by the roaming mechanism in H2CO photodissociation.[6,7] The unsatisfactory results obtained for the H2S + Cl reaction by applying canonical variational transition state theory (CVTST) have been interpreted as another failure of TST.[8] However, the quantum-chemical approach employed there calls for a deeper reinvestigation of the reaction.

To provide a definitive elucidation of the mechanism of the H2S + Cl reaction, in the present work, we have investigated its reactive PES by means of state-of-the-art quantum-chemical computations and coupled them with sophisticated kinetic models still rooted in the TST. In the framework of an ab initio-transition-state-theory-based master equation (AITSTME) treatment, different approaches have been employed to deal with barrierless reactions. The investigation of this specific reaction will lead us to the definition of an accurate protocol for the investigation of the thermochemistry and kinetics of challenging reactions.

In addition to offer a puzzling case study, the H2S + Cl reaction plays an important role in atmospheric chemistry and might be of relevance also for the investigation of other planetary atmospheres. Indeed, the atmospheric sulfur cycle has been the subject of intensive investigations for a long time, mostly because of the need of continuously assessing and monitoring the contribution of anthropogenic sulfur compounds to problems such as acid rain, visibility reduction, and climate modification.[9] In particular, reduced sulfur-containing species are important in the chemistry of the atmosphere; among them, hydrogen sulfide (H2S) is one of the simplest, but yet it plays an important role in earth’s and planetary atmospheres. Concerning the terrestrial environment, its concentration in the atmosphere is significantly due to the decomposition of organic matter and volcanic eruptions, which can inject H2S directly into the stratosphere.[10,11] For example, measurements of H2S concentration by UV spectroscopy at volcanic sites in Italy have shown that this quantity can be on the order of 100 ppm (i.e., much larger than its average atmospheric concentration), H2S thus being 2–3 times more abundant than SO2.[12,13] However, the anthropogenic emission of H2S should not be overlooked.[14,15]

In earth’s atmosphere, H2S is mainly removed by the hydroxyl radical (OH) by means of the gas-phase reaction[16]However, in some marine remote boundary layers and coastal urban areas, the concentration of the chlorine radical (Cl) is larger than that of OH.[17] Therefore, the reaction of H2S with Cl, namelyis also important. Concerning planetary systems, reaction can play a role in Venus’ atmosphere, the latter being rich in H2S.[18−20] In addition to its importance in atmospheric processes, reaction is of great interest as a prototype for heavy–light–heavy atom reactive systems,[21] and, furthermore, it leads to the production of vibrationally excited HCl molecules, which can be used in infrared chemiluminescence and laser-induced fluorescence studies.[22]

Reaction has been studied since 1980, both experimentally and theoretically. Nevertheless, there is still some uncertainty concerning its detailed mechanism. At room temperature, the experimental rate constant spans from 3.7 × 10–11 to 10.5 × 10–11 cm3 molecule−1 s−1.[23−32] According to the review study by Atkinson et al.,[33] the most reliable results are those by Nicovich et al.,[30] who reported an extensive investigation over a wide range of experimental conditions. In that work, the value of the rate constant at room temperature was found to be pressure-independent over a wide range, the latter being 33–800 mbar (i.e., 25–600 Torr).

One proposed mechanism is the direct hydrogen abstraction, the reaction thus proceeding through a transition state. Another possibility is offered by an addition/elimination mechanism, which has been discussed by various groups and nowadays seems to be widely accepted. Despite this, to the best of our knowledge, no theoretical work was able to correctly reproduce and interpret the experimental data. Indeed, the computational work by Wilson and Hirst pointed out the existence of a H2S···Cl adduct, but without being able to connect it with the products, i.e., HS and HCl.[34] In ref (34), the rate constant for the direct abstraction was computed using conventional transition state theory (cTST), which led to a value, at room temperature, of 2.8 × 10–12 cm3 molecule–1 s–1. However, such rate constant is 1 order of magnitude smaller than the experimental datum. Resende et al. instead investigated the addition/elimination mechanism.[8] They obtained a rate constant of 1.2 × 10–9 cm3 molecule–1 s–1, 1 order of magnitude larger than the experimental value, and, as already mentioned above, they ascribed this discrepancy to a “breakdown of transition state theory”.

To solve this puzzle, we have undertaken a comprehensive analysis of the whole reaction mechanism using state-of-the-art electronic structure and kinetic models. The paper is organized as follows. In the next section, the computational methodology is described in some detail, thus introducing the different approaches employed for the electronic structure and kinetic calculations. Then, the results will be reported and discussed: first, the characterization of the reactive PES will be provided, followed by the accurate evaluation of its thermochemistry; then, the reaction rate constants will be addressed. Finally, the major outcomes of this work will be summarized in the concluding remarks.

Computational Methodology

In this section, the methodology employed for the characterization of the H2S + Cl PES and its energetics will be first of all introduced. Then, we will move to the definition of the models used for the accurate interpretation of the kinetic aspects of the title reaction.

Electronic Structure Calculations

Several works have shown that double-hybrid functionals in conjunction with basis sets of at least triple-ζ quality represent a remarkable compromise between accuracy and computational cost.[35−38] For the calculation of equilibrium geometries and vibrational frequencies, the B2PLYP functional often approaches, and in some cases even overcomes, the accuracy of the much more computationally expensive CCSD(T) method, when used in conjunction with comparable basis sets (see, e.g., refs (39) and (40)). CCSD(T), often denoted as the “gold standard” for accurate calculations, stands for the coupled-cluster (CC) method including a full account of single and double excitations, CCSD,[41] and a perturbative estimate of triple excitations (CCSD(T)).[42] In this respect, the recent work by Martin’s group has led to the development of the revDSD-PBEP86 functional.[43] This represents a significant improvement with respect to B2PLYP, especially for transition states and noncovalent interactions, also showing very good performances (as B2PLYP) for equilibrium geometries and, especially, vibrational frequencies. Although the very recent D4 model for dispersion contributions[44] provides some improvement on energy evaluations, the D3(BJ) model[45,46] is already remarkably accurate. Therefore, since a full analytical implementation of second derivatives of energy is available for the latter,[47] we have decided to rely on the D3(BJ) scheme for incorporating dispersion effects. For all of these reasons, in this paper, we have characterized all stationary points of the reactive PES under consideration with the revDSD-PBEP86-D3(BJ) functional in conjunction with the jun-cc-pV(T+d)Z basis set.[48−50]

Subsequently, the energetics of all stationary points was accurately determined by exploiting the composite scheme denoted “HEAT-like”, which is a state-of-the-art approach and will be described in detail in the following. Using this model as reference, the performance of different variants of the so-called “cheap” composite scheme[51,52] (described in the following as well), denoted as ChS, will be investigated. For comparison purposes, the CBS-QB3 model[53,54] will also be considered because it is extensively used in the evaluation of the thermochemistry of reactive systems.

For all levels of theory, the spin–orbit (SO) corrections for the Cl radical as well as for all open-shell species have been computed, within the state-interacting approach implemented in the MOLPRO program,[55−57] at the complete active space self-consistent field (CASSCF)[58,59] level in conjunction with the aug-cc-pVTZ basis set[60,61] and the full valence as active space. For Cl, calculations of the SO corrections have also been carried out using the multireference configuration interaction (MRCI) method[62−64] in conjunction with the aug-cc-pVnZ sets, with n = T, Q, and 5. The computations for Cl have been performed to calibrate the level of theory to be used for the other radicals. Since it has been noted that the CASSCF values are almost independent of the basis set used (from 274.5 cm–1 with aug-cc-pVTZ to 275.7 cm–1 with aug-cc-pV5Z) and very close to the results at the MRCI level (276.1 cm–1 with aug-cc-pVTZ and 279.5 cm–1 with aug-cc-pVQZ), the cheapest level of theory has been chosen for all computations. In passing, we note that only at the MRCI/aug-cc-pV5Z level a value of 295.2 cm–1, in very good agreement with the experimental result of 293.663 cm–1,[65] was obtained. In conclusion, the CASSCF/aug-cc-pVTZ level is expected to provide SO corrections affected by uncertainties not exceeding 0.2 kJ mol–1.

Finally, electronic energies need to be corrected for the zero-point vibrational energy (ZPE) contribution. These corrections have been obtained both within the harmonic approximation and at the anharmonic level. In both cases, they have been computed using the double-hybrid revDSD-PBEP86-D3(BJ) functional in conjunction with the jun-cc-pV(T+d)Z basis set. Second-order vibrational perturbation theory (VPT2)[66] has been exploited for the evaluation of anharmonic ZPEs. Density functional theory (DFT) geometry optimizations and force field computations have been performed with Gaussian 16 quantum-chemical software.[67]

Reference Structures

The first issue that needs to be addressed is the effect of the reference geometries on the energetics. Indeed, a reactive PES can be very complicated and characterized by several alternative mechanisms. Therefore, the search of the stationary points of a reactive PES and their geometry optimizations might become the computational rate-determining step of the investigation. In this view, it is important to rely on a suitable and reliable level of theory for this purpose. The revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level seems indeed to meet the requirements. However, further calculations to check the accuracy obtainable in structural determinations and their suitability for energetics have been performed. Concerning the latter point, the ChS approach applied to geometry optimizations has been employed, with the computational details being provided in the specific section.

Focusing on the reactant-well (RW) adduct (see Figure ), as a test case (a radical species with a noncovalent bond), the equilibrium structure has been accurately evaluated by means of a combination of gradient and geometry approaches entirely based on coupled-cluster (CC) theory.[41] First of all, the so-called CCSD(T)/CBS+CV equilibrium structure has been obtained by minimizing the following gradientwhere the first two terms on the right-hand side are the energy gradients for the extrapolation to the complete basis set (CBS) limit and the last term incorporates the effect of core–valence (CV) correlation. The exponential formula introduced by Feller[68] and the two-point n–3 expression by Helgaker et al.[69] are used for the extrapolation to the CBS limit of the Hartree–Fock self-consistent field (HF-SCF) energy gradient and the CCSD(T) correlation contribution, respectively. The cc-pVnZ basis sets[49,60,70,71] have been employed, with n = Q, 5, and 6 being chosen for the HF-SCF extrapolation and n = Q and 5 for CCSD(T). Since the extrapolation to the CBS limit is performed within the frozen-core (fc) approximation, CV correlation effects have been considered by adding the corresponding correction, dΔECV/dx, where the all-electron–frozen-core energy difference is evaluated employing the cc-pCVQZ basis set.[72,73] It is to be noted that the CCSD(T)/CBS+CV equilibrium structure employing the aug-cc-pVnZ sets (n = T, Q, 5 for HF-SCF and n = T, Q for CCSD(T)) for the extrapolation to the CBS limit and the cc-pCVTZ basis set for the CV contribution has been obtained for all minima: H2S, HS, HCl, and the reactant-well (RW) and product-well (PW) adducts.

Figure 1

Reaction mechanism for the H2S + Cl reaction. SO- and ZPE-corrected HEAT-like energies are reported.

The contributions due to the full treatment of triple (Δr(fT)) and quadruple (Δr(fQ)) excitations have been obtained at the “geometry” level, by adding the following differences to the CCSD(T)/CBS+CV geometrical parameterswhere r denotes a generic structural parameter. The cc-pVTZ basis set has been used for the fT correction and the cc-pVDZ set for the fQ contribution. This implies that geometry optimizations at the fc-CCSDT[74−76]/cc-pVTZ, fc-CCSD(T)/cc-pVTZ, fc-CCSDTQ[77]/cc-pVDZ, and fc-CCSDT/cc-pVDZ levels have been performed.

HEAT-like Approach

The CC-based approach that has been denoted as HEAT-like takes the HEAT protocol[78−80] as reference and starts from the CCSD(T) method. In detail, the scheme (that will be denoted as “CBS+CV+DBOC+rel+fT+fQ” in the following) can be summarized as followsIn the expression above, “CBS” means that CCSD(T) energies, obtained within the fc approximation, have been extrapolated to the CBS limit. Analogously to what has been done for the CCSD(T)/CBS+CV gradient scheme, the extrapolation to the CBS limit has been performed in two steps, i.e., HF-SCF and the CCSD(T) correlation energies have been extrapolated separately. The HF-SCF CBS limit has been evaluated by exploiting the exponential expression introduced by Feller[68]For the CCSD(T) correlation contribution, the extrapolation to the CBS limit has been carried out using the n–3 formula by Helgaker and co-workers[69]The cc-pVQZ and cc-pV5Z basis sets have been employed for the CCSD(T) correlation energy, whereas the cc-pVnZ sets, with n = Q, 5, and 6, have been used for HF-SCF.

By making use of the additivity approximation, the CV effects have been taken into account by means of the following expressionthus incorporating the CCSD(T) energy difference obtained from all electrons (ae) and fc calculations, both in the cc-pCVQZ basis set.[72,73]

The diagonal Born–Oppenheimer correction, ΔEDBOC ,[81−84][81−84] and the scalar relativistic contribution to the energy, ΔErel , have been computed at the HF-SCF/aug-cc-pVTZ and MP2/unc-cc-pCVQZ (correlating all electrons) levels, respectively, where MP2 stands for Møller–Plesset theory to second order[85] and “unc” denotes the use of the uncontracted basis set. The contributions due to relativistic effects have been evaluated using the lowest-order direct perturbation theory (second order in 1/c, DPT2).[86]

In analogy to eq , corrections due to a full treatment of triples, ΔEfT , and of quadruples, ΔEfQ , have computed, within the fc approximation, as energy differences between CCSDT and CCSD(T) and between CCSDTQ and CCSDT calculations employing the cc-pVTZ and cc-pVDZ basis sets, respectively.

In addition to the scheme described above, a variant including the less expensive CCSDT(Q) method[87−89] has also been consideredwhere the only difference with respect to the former approach is the use of CCSDT(Q), which incorporates the quadruple excitations in a perturbative manner, instead of CCSDTQ. This approach will be denoted as “CBS+CV+DBOC+rel+fT+pQ”.

It should be noted that the HEAT-like schemes also provide all possible intermediate approaches, such as CCSD(T)/CBS, CCSD(T)/CBS+CV, and CCSD(T)/CBS+CV+DBOC+rel, together with the results for the single-basis calculations like fc-CCSD(T)/cc-pVQZ and fc-CCSD(T)/cc-pV5Z.

All computations within the HEAT-like schemes have been carried out using the quantum-chemical CFOUR program package,[90] with the MRCC code[91] being interfaced to CFOUR to perform the calculations including quadruple excitations.

For all schemes described above, the geometries of the stationary points have been optimized using the double-hybrid revDSD-PBEP86-D3(BJ) in conjunction with the jun-cc-pV(T+d)Z basis set.

ChS Variants

The so-called cheap geometry approach (as already mentioned, denoted as ChS) was initially developed for accurate molecular structure determinations[51,92] and then extended to energetic evaluations[52] for medium-sized systems. Recently, this scheme has been improved to accurately describe intermolecular, noncovalent interactions (thus leading to the definition of the jun-ChS approach).[93]

The starting point of this composite scheme is the CCSD(T) method in conjunction with a basis set of triple-ζ quality within the fc approximation. To improve the accuracy of this level of theory, the ChS model requires the incorporation of the extrapolation to complete basis set (CBS) and the core–valence (CV) correlation contribution, with all of them taken into account using MP2. Therefore, the general ChS model can be described by the following expressionIn all cases, the extrapolation to the CBS limit has been performed:

In two steps, as done in the framework of the HEAT-like approach. For HF-SCF, n = T, Q, and 5 have been used, while n = T and Q have been considered for MP2.

In one step, according tothus employing triple- and quadruple-ζ basis sets.

In analogy to the HEAT-like approach, the CV contribution, ΔECV, has been evaluated as energy (at the MP2 level) difference between ae and fc calculations. A basis set of triple-ζ quality has been employed for all ChS models.

The ChS variants employed in the present study can be classified as follows:For both standard ChS approaches, d-augmented basis sets have been employed for CCSD(T) calculations as well as for the extrapolation to the CBS limit: the cc-pV(n+d)Z[49,50,60] basis sets for ChS and jun-cc-pV(n+d)Z[48−50] for jun-ChS. The CV correlation correction has been evaluated employing the weighted-core–valence cc-pwCVTZ basis set.[73]

Standard, based on standard CCSD(T) and MP2 methods: ChS[51,52] and jun-ChS.[93]

F12, based on explicitly correlated F12 methods: ChS-F12.

While the standard ChS approaches have already been largely applied (the reader is referred to the cited references for more details), the ChS-F12 model has been introduced for the first time. In detail, the general expression of eq becomes the followingFor evaluating the CC term (using CCSD(T)-F12[94−96] within the fc approximation), the cc-pVTZ-F12 basis set[97] has been used. As in the case of ChS and jun-ChS, the extrapolation to the CBS limit has been performed both in one and two steps. In the case of the two-step procedure, the HF/CBS contribution has been taken from ChS, while for the extrapolation of the correlation contribution with the MP2-F12 method,[98] two sets of basis sets have been considered:[97] cc-pVDZ-F12/cc-pVTZ-F12 (extrapolation D,T) and cc-pVTZ-F12/cc-pVQZ-F12 (extrapolation T,Q). For the one-step extrapolation to the CBS limit, the two sets of bases have been used as well. In the case of the D,T extrapolation, formula has been modified accordingly. The CV contribution has been computed using the cc-pCVTZ-F12 basis set.[99] Furthermore, both the F12a and F12b variants[94] have been employed, and, as described in ref (100), the geminal exponent (β) was set to 1.0 for all the F12 basis sets considered.

The final remark concerns the reference geometries used in the energy evaluations. In the case of the standard ChS approaches, the molecular structures employed have been obtained at the same level of theory, the only exception being the tests made on the effect of the reference geometries. This means that an expression analogous to eq has been exploited to obtain the geometries of the stationary points. The only difference lies in the fact that, for structural determinations, the extrapolation to the CBS limit has been performed only in one step, according to eq . Conversely, for ChS-F12, the reference geometries have been determined at the fc-CCSD(T)-F12/cc-pVDZ-F12 level. As mentioned above, to investigate the structural effects on energetics, the jun-ChS model has also been applied to revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z geometries.

For the ChS approaches, the geometry optimizations and energy evaluations have been performed with the MOLPRO quantum-chemistry package.[55,56]

Kinetic Models

The reactive PES for the H2S + Cl reaction is summarized in Figure . As evident from this figure, the Cl atom reacts with H2S, thus leading to the H2S···Cl intermediate well, denoted as RW. From here onward, the only available channel is, at least under atmospheric conditions, the addition/elimination reaction: RW isomerizes via the transition state, TS, to the H-bonded HS···HCl product well, denoted as PW, which then evolves to the products, i.e., HS + HCl. This process can be described in the framework of the AITSTME approach through a three-channel, two-well master equation. Depending on the examined temperature and pressure conditions, the global reaction rate is controlled either by the rate of conversion of RW into PW or by both its rate and that of formation of RW. The rate of decomposition of PW is generally fast with respect to the two other reaction channels and thus does not impact the global reaction rate.

The rate of formation of RW, which proceeds over a barrierless PES, has been determined at three different levels of theory. At the lowest level, the rate constant has been determined using phase space theory (PST),[2] assuming a attractive potential, with the coefficient C obtained by fitting the energies computed at various long-range distances (R) of the fragments using the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level of theory. Alternatively, the rate constant has been determined using variational transition state theory (VTST) and variable reaction coordinate transition state theory (VRC-TST),[2] which, among the considered theoretical approaches, is the most accurate to describe properly the large amplitude motions that characterize this reaction channel.

VTST calculations have been performed within the rigid-rotor harmonic-oscillator (RRHO) approximation over a PES scanned as a function of the distance between Cl and H2S in the 3.0–4.6 Å interval using a 0.2 Å step. Structures and vibrational frequencies of each point along the PES have been determined at the CASPT2/cc-pVDZ[101−103] level using a (9e,7o) active space composed of the four H–S σ bonding and antibonding orbitals (4e,4o) and of the three p valence orbitals of chlorine (5e,3o). Higher-level energies have been determined at the CASPT2/aug-cc-pVTZ level over a (21e,13o) active space, consisting of the same active space used for the geometry optimization with the addition of the remaining valence electrons of Cl and S and of the three 2p orbitals and electrons of S. All CASPT2 calculations have been performed by state averaging over three states using a 0.2 level shift, the only exception being made for high-level calculations, which used a 0.25 IPEA shift. The IPEA shift was chosen as it allows for reproducing with relatively good accuracy the HEAT-like electronic energy of the entrance van der Waals well RW: −43.7 vs −41.89 kJ mol–1.

VRC-TST calculations have been performed sampling the PES over multifaceted dividing surfaces constructed using three pivot points, positioned as schematized in Figure . Two pivot points (P11 and P12) were placed in proximity of the S atom, symmetrically displaced along the axis perpendicular to the H2S plane and passing from S by a distance d, while the third pivot point was centered on Cl. The multifaceted dividing surface is constructed varying the distance r between Cl and the pivot points, as described by Georgievskii and Klippenstein in ref (104), between 5 and 11 a0 with 1.0 a0 step. Calculations were repeated changing the position of the P1 pivot points varying d between 0.01 and 0.6 a0, using 0.1 a0 steps. Reactive fluxes were computed through Monte Carlo sampling using a 5% convergence threshold.

Figure 2

Scheme describing the position of the pivot points used to construct the multifaceted dividing surface for VRC-TST calculations. The H2S plane is perpendicular to the plane of the picture, and the second H atom is hidden behind the S atom.

A minimum of 200 sampling points have been taken for each dividing surface. The calculated reactive fluxes have then been multiplied by a flat 0.9 factor to correct for recrossing dynamical effects. The 0.9 correction factor comes from the comparison of VRC-TST with trajectory calculations, which showed that VRC-TST total rate coefficients are generally about 10% greater than those obtained from related trajectory simulations,[105] when VRC-TST is applied at the level of theory used in the present study.[106] Energies have been computed at the CASPT2/cc-pVDZ level using a (5e,3o) active space consisting of the p orbitals of Cl, keeping the structures of H2S frozen in its minimum-energy configuration and state averaging over three states. The sampled energies have been corrected, using a one-dimensional potential function of the distance between the S and Cl atoms, for geometry relaxation, energy accuracy, and SO effects. The correction energy for geometry relaxation (ΔEGEOM) has been determined at the level of theory used to optimize geometries for VTST calculations, while that for the energy accuracy (ΔEHL) has been estimated at the higher level used for VTST calculations (CASPT2/aug-cc-pVTZ level over a (21e,13o) active space). Spin–orbit corrections along the PES (ΔESO) have been computed using the state interacting method with a Breit–Pauli Hamiltonian and a CASSCF wave function determined with a (9e,7o) active space state and six electronic states. The VRC-TST energy is thus computed asThe rate of conversion of RW into PW has been calculated using both conventional and variational TST, as the reaction proceeds through a well-defined saddle point. Energies and structures of the stationary points (wells and saddle point) have been evaluated as described in Section . VTST calculations have been performed computing frequencies and structures along a minimum-energy path (MEP) determined through intrinsic reaction coordinate (IRC) calculations, at the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level, taking 10 steps of 0.03 a0 in both directions. Frequencies along the MEP have been computed both in Cartesian and in internal curvilinear coordinates.[107] Tunneling corrections have been taken into account by means of both small curvature theory (SCT)[108] and the Eckart model,[109] using an asymmetric potential.

VRC-TST calculations have been performed using the VaReCoF software,[106,110] while master equation simulations have been carried out using MESS.[111] VTST calculations (in Cartesian and internal curvilinear coordinates) as well as the determination of stationary points along the MEP to determine the VRC-TST potential and the construction of the input files for VaReCoF computations have been performed using a modified version of EStokTP.[112] All CASPT2, CASSCF, and SO calculations have been carried out using the MOLPRO program.[56] For all TST variants, where required, the energies of the stationary points of the PES have been considered at the CBS+CV+DBOC+rel+fT+fQ level, also including SO corrections and anharmonic ZPE contributions.

Results and Discussion

The reaction mechanism for the H2S + Cl reaction is shown in Figure and, as mentioned in the Introduction, none of the previous computational works provided a complete picture for it. In particular, the presence of both H2S···Cl and HS···HCl has been thoroughly investigated only in the present study. In the following, the molecular structures of the stationary points are first presented. Then, the associated thermochemistry and kinetics are reported and discussed.

Molecular Structures

The structural parameters of the stationary points located on the H2S + Cl reactive PES, optimized at different levels of theory, are collected in Table . From the inspection of this table, the first conclusion that can be drawn is that for covalently bonded compounds, i.e., H2S, HS, and HCl, there is a perfect agreement among ChS, jun-ChS, CCSD(T)/CBS+CV, and revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z results. For the intermediate adducts, i.e., RW and PW, we note again that the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z covalent bond lengths agree very well with those determined by means of composite schemes. The discrepancies are evident only for the noncovalent distances, i.e., S···Cl in RW and S···H in PW. However, the differences are well within 0.01–0.02 Å, meaning that their impact on energetics is expected to be negligible (as will be demonstrated in the next section).

Table 1

Structural Parameters of the Stationary Points of the H2S + Cl Reaction at Different Levels of Theoryg

		ChS	jun-ChS	CC-F12a	DSD-D3b	CBS+CVc	CBS+CV+fT+fQd	QCISDe	experimentf
H2S	r(H–S)	1.336	1.336	1.337	1.338	1.335	1.335		1.3356
				(1.338)		(1.335)
	θ(H–S–H)	92.2	92.2	92.2	92.5	92.3	92.3		92.11
				(92.3)		(92.3)
H2S···Cl	r(H–S)	1.337	1.336	1.338	1.339	1.336	1.336	1.337
				(1.339)		(1.336)
	r(S–Cl)	2.582	2.586	2.585	2.595	2.567	2.568	2.670
				(2.584)		(2.585)
	θ(H–S–H)	93.0	92.8	92.6	92.9	92.9	92.9
				(92.8)		(92.8)
	θ(H–S–Cl)	87.4	87.6	87.5	88.0	87.4	87.4	88.2
				(87.5)		(87.5)
TS	r(H1–Cl)	1.645	1.644	1.630	1.633	1.642		1.614
				(1.642)		(1.635)
	r(H1–S)	1.470	1.469	1.476	1.472	1.468		1.478
				(1.471)		(1.469)
	r(H2–S)	1.339	1.339	1.340	1.341	1.340		1.339
				(1.341)		(1.338)
	θ(Cl–H1–S)	127.0	128.2	129.8	129.6	127.2		137.0
				(128.0)		(127.4)
	θ(H1–S–H2)	90.7	90.7	90.6	91.3	91.0
				(90.8)		(90.8)
	φ(Cl–H1–S–H2)	281.0	281.0	281.4	280.4	280.9
				(281.0)		(281.1)
HS···HCl	r(H2–S)	2.508	2.506	2.484	2.481	2.492	2.492
				(2.506)		(2.492)
	r(H1–S)	1.341	1.341	1.343	1.343	1.341	1.341
				(1.343)		(1.341)
	r(H2–Cl)	1.284	1.284	1.285	1.288	1.284	1.284
				(1.286)		(1.284)
	θ(H1–S–H2)	92.4	92.3	91.8	92.6	92.0	92.0
				(91.8)		(91.9)
	θ(Cl–H2–S)	176.5	176.1	175.8	176.6	175.8	175.8
				(176.0)		(175.8)
HS	r(H–S)	1.340	1.340	1.341	1.342	1.340	1.340		1.3406194(3)
				(1.342)		(1.340)
HCl	r(H–Cl)	1.274	1.274	1.274	1.276	1.274	1.274		1.274565598(53)
				(1.276)		(1.274)

CC-F12 stands for fc-CCSD(T)-F12 in conjunction with the cc-pVDZ-F12 basis set. Values within parentheses have been obtained in conjunction with the cc-pVTZ-F12 basis set.

DSD stands for revDSD-PBEP86-D3(BJ) in conjunction with the jun-cc-pV(T+d)Z basis set.

CBS+CV stands for CCSD(T)/CBS+CV, with the aug-cc-pVnZ sets (n = T, Q) used for the extrapolation to the CBS limit and cc-pCVTZ for the evaluation of the CV contribution. Within parentheses are given the results for cc-pVQZ and cc-pV5Z being used for CBS and cc-pCVQZ for CV. See text.

CBS+CV+fT+fQ stands for CCSD(T)/CBS+CV augmented by fT and fQ contributions. See text.

In conjunction with the cc-pV(T+d)Z basis set. Values are taken from ref (8).

H2S: ref (113); HCl: ref (114); HS: ref (115).

Distances are in angstroms, and angles are in degrees.

To further test the performance of different ChS approaches and of the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level, as mentioned in the Section , the CCSD(T)/CBS+CV+fT+fQ scheme has been exploited for the RW adduct. This allowed us to confirm beyond all doubts the accuracy of the ChS models and the suitability of the revDSD-PBEP86-D3(BJ) functional for the characterization of noncovalently bonded systems. Moreover, the results of Table point out the very good performance of the CCSD(T)-F12/cc-pVDZ-F12 level of theory, which is only marginally improved by the use of the cc-pVTZ-F12 basis set. Indeed, also in the case of fc-CCSD(T)-F12 calculations, a very good agreement with the structural parameters obtained from composite approaches can be noted.

Finally, it is apparent that the results at the QCISD/cc-pV(T+d)Z level from ref (8) show remarkable deficiencies in the description of the geometrical parameters. In particular, the θ(Cl–H1–S) angle in the transition state is about 10° larger than the corresponding best-estimated value. Another important deviation is observed for the RW adduct, the S···Cl distance being too long by ∼0.1 Å and the θ(H–S–Cl) too large by ∼1°.

H2S + Cl Reaction: Previous Works

In the following, a critical analysis of the previous computational investigations on the H2S + Cl reaction is reported.

In the work by Wilson et al.,[34] the RW and TS stationary points were characterized by evaluating the electronic energies at the ae-MP4/6-311+G(2df,p) level of theory,[116,117] using ae-MP2/6-311G** reference structures. Moreover, the addition and abstraction pathways were treated as separate processes instead of combining them in an addition/elimination mechanism, as done in this work. As a consequence, the rate constant was calculated with standard TST considering only the abstraction step. However, a value 1 order of magnitude lower than the experimental one is mainly explained by an incorrect evaluation of the barrier height.

In a more recent work, in which the PW adduct was not considered as well, Resende et al.[8] optimized the geometries of the stationary points at the QCISD/cc-pV(T+d)Z level of theory.[118,119] Using these reference structures, the energies were subsequently evaluated at the PMP2 level of theory[120−123] in conjunction with cc-pV(n+d) basis sets (with n = D, T, Q), extrapolated to the CBS limit using the exponential expression introduced by Feller, and finally added to the corresponding ae-CCSD(T)/cc-pV(T+d)Z energies. Resende et al. thus employed for the energetics a composite scheme similar to our ChS. Using the aforementioned geometries and energies, the rate constant was subsequently calculated with VTST, thereby obtaining a value that was 1 order of magnitude higher than the experimental datum. To justify their overestimated result, the authors advocated a failure of transition state theory. The role of explicit dynamical effects (e.g., recrossing) was analyzed by performing some explicit trajectory computations. Although the results obtained were not conclusive, they induced the authors to advocate the inadequacy of TST in reproducing the significant role of vibrationally excited H2S in stabilizing the H2S···Cl adduct.

H2S + Cl Reaction: Thermochemistry

The relative electronic energies, obtained at various computational levels, are detailed in Table . In particular, fc-CCSD(T) results in conjunction with basis sets of increasing size up to the CBS limit are collected together with the CBS+CV, CBS+CV+DBOC, CBS+CV+DBOC+rel, and CBS+CV+DBOC+rel+fT+fQ values. These allow us to inspect the trend of the relative energies as a function of the basis set as well as the role of the various contributions.

Table 2

Relative Electronic Energiesa for the H2S + Cl Reactiong

	reactants	RW	TS	PW	products
	H2S + Cl	H2S···Cl		HS···HCl	HS + HCl
CCSD(T)/VTZ	0.00	–28.29	7.74	–56.10	–44.90
CCSD(T)/VQZ	0.00	–37.57	1.30	–58.23	–46.35
CCSD(T)/V5Z	0.00	–41.32	–0.84	–59.11	–46.95
CCSD(T)/CBS	0.00	–44.84	–3.10	–60.10	–47.51
CBS+CV	0.00	–45.09	–3.42	–60.30	–47.65
CBS+CV+DBOC	0.00	–45.09	–2.05	–60.15	–47.64
CBS+CV+DBOC+rel	0.00	–45.11	–2.31	–59.98	–47.45
CBS+CV+DBOC+rel+fT+pQ	0.00	–45.14	–3.41	–60.09	–47.45
CBS+CV+DBOC+rel+fT+fQ	0.00	–45.11	–3.33	–60.08	–47.44
ChSb	0.00	–42.81	–2.49	–60.29	–47.94
		[−42.94]	[−2.54]	[−60.41]	[−47.44]
jun-ChSb	0.00	–43.89	–3.46	–60.19	–47.58
		(−43.81)	(−3.45)	(−60.18)	(−47.59)
		[−43.89]	[−3.46]	[−60.19]	[−47.58]
ChS-F12a/F12b CBS(D,T)c	0.00	–44.64/–44.58	–1.97/–1.93	–59.80/–59.80	–47.11/–47.15
		[−43.97/–43.91]	[−2.74/–2.70]	[−60.10/–60.10]	[−47.06/–47.09]
ChS-F12a/F12b CBS(T,Q)c	0.00	–43.87/–43.81	–3.32/–3.29	–60.33/–60.33	–47.42/–47.46
		[−43.54/–43.48]	[−3.39/–3.35]	[−60.27/–60.27]	[−47.28/–47.31]
CBS-QB3d	0.00	–48.16	–8.12	–61.13	–53.43
revDSD-PBEP86-D3(BJ)e	0.00	–47.82	–3.85	–63.09	–46.11
		(−45.50)	(−1.80)	(−59.43)	(−46.00)
QCISD/cc-pV(T+d)Zf	0.00	–23.14	20.17		–43.76
PMP2/CBSf	0.00	–46.69	–9.08		–52.72
ae-CCSD(T)/CBSf	0.00	–46.48	–7.57		–48.16

Unless otherwise stated, reference geometries at the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level.

Reference geometries at the same level as the energy evaluation. Extrapolation to the CBS limit in one step. Values within parentheses: revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z geometries as reference. Values within square brackets: extrapolation to the CBS limit in two steps.

Reference geometries at the fc-CCSD(T)-F12/cc-pVDZ-F12 level. Extrapolation to the CBS limit in one step. Values within square brackets: Extrapolation to the CBS limit in two steps

Reference geometries at the B3LYP/6-31G(d) level.

Values obtained in conjunction with the jun-cc-pV(Q+d)Z basis set. Within parentheses: results for the jun-cc-pV(T+d)Z basis set. In both cases: reference structures at the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level.

Results from ref (8).

Values in kJ mol–1.

For all stationary points, it is noted that even for a basis set as large as cc-pV5Z the results are not yet converged, with values being quantitatively accurate only at the CBS limit. A particular remark is warranted for the fc-CCSD(T)/cc-pVTZ level because it is often used in the investigation of reactive PESs, in particular for astrophysical purposes. This level predicts the RW adduct to lie ∼28 kJ mol–1 below the reactants, which means less stable by about 17 kJ mol–1 with respect to what evaluated by more accurate calculations. Furthermore, at this level, the transition state is predicted to emerge above the reactants by more than 7 kJ mol–1. By enlarging the basis set, we note that TS is still emerged with the quadruple-ζ set and becomes barely submerged only when the quintuple-ζ basis is used.

Moving from the fc-CCSD(T)/CBS level on, it is observed that the CV corrections are small, these being, on average, of the order of 0.2 kJ mol–1. Roughly of the same order of magnitude are the combined DBOC and scalar relativistic contributions, with the only exception being the transition state, for which the correction amounts to about 1 kJ mol–1. Analogous is the situation for the fT and fQ contributions, with the corresponding correction for TS being, however, in the opposite direction. The overall conclusion is that the CBS+CV level provides results in very good agreement with the CBS+CV+DBOC+rel+fT+fQ model. Furthermore, it is noted that when basis sets more suitable for describing the third-row elements are used (i.e., cc-pV(n+d)Z), for the CBS+CV model, smaller basis sets can be employed. In another test, the extrapolation to the CBS limit has been applied to ae-CCSD(T)/cc-pCVnZ (n = Q, 5) energies. Since the differences with the CBS+CV level lie within 0.2 kJ mol–1, the validity of the additivity approximation for the CV contribution has been confirmed. From Table , it is also evident that the perturbative (instead of full) treatment of quadruple excitations marginally affects the relative energies.

Moving to less expensive approaches, the very good performance of all ChS variants deserves to be highlighted. Indeed, ChS, jun-ChS, and ChS-F12 provide very similar results, which deviate, in the cases, by about 2 kJ mol–1 from our best estimates (i.e., CBS+CV+DBOC+rel+fT+fQ). The jun-ChS values obtained using also the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z reference geometries are reported, within parentheses, in Table . A perfect agreement between the two sets of jun-ChS relative energies is observed, the differences being well below 0.1 kJ mol–1. Such a comparison confirms the accuracy of the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z structures for thermochemical studies. In addition, the revDSD-PBEP86-D3(BJ) functional provides accurate results also for the energetics, as clear from the results collected in Table . We note that this functional performs better when combined with the jun-cc-pV(T+d)Z basis set for all minima, while the jun-cc-pV(Q+d)Z set seems to be required for correctly evaluating the relative energy of the transition state.

In Table , for all ChS variants, relative energies based on extrapolations to the CBS limit performed in both one and two steps are reported. It is worth noting that the two approaches provide very similar results. This is an important outcome because the extrapolation in one step avoids the problems related to basis sets of quintuple-ζ quality for which convergence of the HF-SCF energy can be troublesome, especially for open-shell species. Furthermore, the computation of the corresponding integrals might become particularly expensive for large systems.

As mentioned in Section , both F12a and F12b variants have been employed in the ChS-F12 model. From the results of Table , it is evident that the two approximations provide nearly coincident results. For this scheme, the extrapolation to the CBS limit using double- and triple-ζ basis sets has also been tested. It is apparent that even in this case the results are very good, showing deviations from the best-estimated values well within 2 kJ mol–1. However, this ChS variant is the only one presenting some difference (about 1 kJ mol–1) between the one- and two-step procedures for the extrapolation to the CBS limit.

As already pointed out in the literature (see, e.g., refs (124) and (125)), despite its widespread use, the CBS-QB3 model provides disappointing results for energetics and, in particular, for barrier heights. In fact, the transition state lies too low in energy by about 5 kJ mol–1 with respect to the best estimate.

Finally, the comparison of our results with those from ref (8) is deserved. First, it is noted that the QCISD/cc-pV(T+d)Z level, whose unsuitability in the determination of reference structures has been previously pointed out, provides unreliable energetics, the transition state being predicted ∼20 kJ mol–1 above the reactants. The situation improves on moving to the PMP2 level and, in particular, to the level denoted in Table as ae-CCSD(T)/CBS, which corresponds, as previously described, to a composite scheme similar to the ChS model. Despite this improvement, the barrier height is underestimated by about 3 kJ mol–1.

To complete the thermochemistry of the H2S + Cl reaction, the relative energies need to be corrected for SO coupling (which was neglected in ref (8)). Doing so for our best level of theory leads to the values collected in the first row of Table . The major consequence is that the transition state is no longer submerged and lies above the reactants by 0.12 kJ mol–1. However, a further correction should be taken into consideration. Once ZPE is incorporated, the transition state is again submerged by about 5 kJ mol–1. Furthermore, it has to be noted that ZPE corrections evaluated within the harmonic approximation are in very good agreement with the anharmonic values.

Table 3

Best-Estimated Relative Electronic Energies (Including Spin–Orbit) together with ZPE and Thermochemical Correctionsd

	reactants	RW	TS	PW	products
	H2S + Cl	H2S···Cl		HS···HCl	HS + HCl
CBS+CV+DBOC+rel+fT+fQa	0.00	–41.89	–0.05	–56.79	–46.25
	(3.28)	(0.06)	(0.01)	(0.00)	(2.10)
anharm-ZPEb	0.00	5.54	–5.85	–1.43	–5.61
harm-ZPEb	0.00	5.86	–5.27	–1.55	–5.86
ΔH° – ΔH0° c	0.00	–2.78	–3.44	–0.58	0.94

Within parentheses, the SO corrections (at the CASSCF/aug-cc-pVTZ level) are given.

Relative ZPE corrections at the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level.

Standard state: 1 atm, 298 K; at the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level.

Values in kJ mol–1.

The computed standard formation enthalpies (298.15 K, 1 atm) of H2S and products (HS and HCl), obtained from our electronic structure calculations and the experimental formation enthalpies of the H, S, and Cl atomic species,[126] are in remarkable agreement with the most accurate experimental data taken from ref (126), except for that of HS,[127] namely, 141.86 vs 141.87, −91.82 vs −92.17, and −20.34 vs −20.50 kJ mol–1 for HS, HCl, and H2S, respectively. All details are provided in the Supporting Information (SI) (see Table S1). From these results, we obtain a sub-kJ mol–1 accuracy also for the reaction enthalpy of the title reaction (50.92 vs 51.24 kJ mol–1).

H2S + Cl Reaction: Rate Constants

The rate constant for the H abstraction from H2S by Cl has been computed at different levels of theory, as described in detail in Section , by solving the multiwell one-dimensional master equation using the chemically significant eigenvalue (CSE) method within the Rice–Ramsperger–Kassel–Marcus (RRKM) approximation, as detailed by Miller and Klippenstein in ref (128).

Although the use of a master equation model is generally not necessary to evaluate the rate constant of an abstraction reaction, there are several reasons that warrant its employment, as described by Cavallotti et al. in ref (112). In fact, it allows one to properly limit the contributions to the reactive flux from energy states below the asymptotic energies of the fragment, as well as to account for limitations to the reaction fluxes determined by an eventually slow rate of formation of the precursor complex. This is indeed the situation for the present system, for which the flux through the outer transition state, leading to the formation of the precursor complex, and that through the inner transition state, leading to H abstraction and the formation of the product complex, have comparable values. The consequence is that part of the flux passing through the outer transition state is reflected from the inner transition state.

As mentioned in Section , reactive fluxes through the outer transition state, which occur on a barrierless MEP, have been computed at three different levels of theory: VRC-TST, VTST, and PST. The Cl–H2S interaction potential used to determine the VRC-TST flux uncorrected for geometry relaxation, active-space size and basis-set size, and at different levels of corrections is reported in Figure . It is interesting to observe that the corrected CASPT2 potential used in VRC-TST calculations is similar, though not exactly overlapped with that determined at the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level of theory. This suggests that this functional may be used to perform VRC-TST calculations for open-shell systems for which it may be cumbersome, e.g., for difficulties in converging to a proper active space, to determine the interaction potential at the CASPT2 level.

Figure 3

Interaction potential between Cl and H2S calculated at different levels of theory: the uncorrected CASPT2/cc-pVDZ level; the CASPT2/cc-pVDZ level corrected for geometry relaxation, high-level energy contributions, and SO effects (see eq ); the CASPT2/cc-pVDZ level corrected for geometry relaxation and high-level energy contributions; and the revDSD-PBEP86-D3(BJ)/jun-cc-pV(T+d)Z level corrected for SO effects.

As already mentioned, master equation simulations have been performed using MESS software. The collisional energy transfer probability has been described by means of the single exponential down model[129] with a temperature dependence ⟨ΔE⟩down of 260(T/298)0.875 cm–1 in an argon bath gas. Different models have been employed to compute reaction fluxes through the inner and the outer transition states. The highest-level simulations have been obtained using VRC-TST for the outer TS and VTST in curvilinear internal coordinates for the inner TS (referred as VTSTin), the results being reported in Table and compared with experimental data in Figure .

Table 4

Rate Coefficients for the H2S + Cl Reaction at Various Temperatures (Pressure = 1 atm)a,b

T (K)	VRC-VTSTin	VTST-VTSTin	PST-VTSTin
200	1.01 × 10–10	1.02 × 10–10	1.31 × 10–10
225	9.49 × 10–11	9.53 × 10–11	1.22 × 10–10
250	8.94 × 10–11	9.03 × 10–11	1.14 × 10–10
275	8.51 × 10–11	8.64 × 10–11	1.08 × 10–10
300	8.16 × 10–11	8.33 × 10–11	1.03 × 10–10
325	7.89 × 10–11	8.07 × 10–11	9.84 × 10–11
350	7.63 × 10–11	7.87 × 10–11	9.50 × 10–11
375	7.44 × 10–11	7.69 × 10–11	9.21 × 10–11
400	7.26 × 10–11	7.56 × 10–11	8.98 × 10–11
425	7.13 × 10–11	7.45 × 10–11	8.78 × 10–11
450	6.99 × 10–11	7.36 × 10–11	8.62 × 10–11
475	6.89 × 10–11	7.29 × 10–11	8.49 × 10–11
500	6.80 × 10–11	7.24 × 10–11	8.38 × 10–11
600	6.56 × 10–11	7.17 × 10–11	8.16 × 10–11
700	6.45 × 10–11	7.24 × 10–11	8.16 × 10–11
800	6.41 × 10–11	7.41 × 10–11	8.31 × 10–11
900	6.44 × 10–11	7.65 × 10–11	8.57 × 10–11

The various prefixes stand for the theoretical methods used to handle the barrierless entrance channel, while the VTSTin suffix means that the inner TS is handled with VTST in curvilinear internal coordinates. The barrierless exit channel is always treated with PST.

Values in cm3 molecule–1 s–1.

Figure 4

H2S + Cl global rate constant: comparison between computed (VRC-TST theory for the outer TS, VTST with vibrational frequencies evaluated using curvilinear internal coordinates for the inner TS, small curvature theory for tunneling) and experimental data.

The comparison between calculated and experimental data shows an excellent agreement at 300 K, the temperature at which most of the measurements were made. Indeed, the calculated 7.76 × 10–11 cm3 molecule–1 s–1 value is in quite good agreement with the 7.4 × 10–11 cm3 molecule–1 s–1 datum recommended by Atkinson et al.[33] on the basis of an extensive review. Furthermore, the calculated rate is in excellent agreement with the rate constant measured in the 200–433 K temperature range by Nicovich et al.,[30] from which it differs by about 12% at most. The calculated global rate constant is almost pressure-independent in the considered conditions, an outcome that agrees with the experimental data that show that there is no measurable collisional stabilization of the entrance well.[32] However, the temperature trend is not perfectly reproduced, as the experimental data are slightly underestimated at low temperatures and slightly overestimated at high temperatures. The discrepancy is more evident for the recent measurements by Gao et al.,[32] which are overestimated by a factor of 1.5 at 900 K. Even if such a disagreement is relatively small and it has been observed only with respect to a single set of experimental data, it is anyway useful to try to understand its origin. It is first of all noted that the structure of the saddle point of the inner TS has an optical isomer, which thus effectively doubles the density of states (DOS) of the TS and, consequently, the rate constant. The two optical isomers are separated by a second-order saddle point with a barrier of about 10 kJ mol–1. Therefore, it is likely that, as the temperature increases, the two isomers interconvert among themselves. If this is the case, then the DOS of the TS is overestimated by up to a factor of 2. To investigate whether this can be the case, the one-dimensional (1D) PES for the conversion between the two isomers has been determined as a function of the Cl–H–S–H dihedral angle and the partition function of the corresponding vibrational internal motion has been replaced with a 1D hindered rotor model. The rate constants calculated at different temperatures are compared with the experimental temperature-dependent data in Figure . The computed results are now in quantitative agreement with experiments at high temperature and differ at most by a factor of 1.38 at 200 K. This outcome confirms, as it is well known in the literature, that one of the key aspects in the estimation of an accurate rate constant using TST is the proper description of anharmonic internal motions and that it is sometimes necessary to use different models depending on the investigated temperature and pressure conditions.

Figure 5

H2S + Cl global rate constant: comparison between computed (same level as in Figure , but modeling of the internal motion for the interconversion between the optical isomers of the TS with a 1D hindered rotor) and experimental data.

To compare the contribution of the entrance and inner channels to the global rate constant, it is interesting to report the rates of each channel computed solving the ME fictitiously enhancing the rate of the other channel (Figure ). As it can be observed, both rates exhibit a negative activation energy, in agreement with experimental observations, and their values are comparable, though the rate of the inner channel is smaller, and it thus impacts more significantly the global reaction flux.

Figure 6

Rate constants of the outer and inner channels computed using VRC-TST and VTST, respectively.

The impact of the level of theory chosen to compute the inner and outer TS fluxes on the global rate constant is analyzed in Figure . It can be observed that, for the outer channel, VRC-TST and VTST give similar results that differ at most by a factor of 1.1, while PST predictions deviate by up to a factor of 1.2. Despite this, it is interesting to notice that there is a slight but significant qualitative difference in the temperature dependence between the rate constants computed using VRC-TST, in better agreement with the experimental trend, and those determined at the other theoretical levels. It should also be recalled that the global rate constant is mostly controlled by the rate of the inner TS, so that the differences between the levels of theory used for the outer TSs are mitigated. The analysis of the impact of the chosen theoretical level for the inner TS shows that variational effects have a minor impact, though the rate constant computed using the internal coordinate model is in better agreement with experimental data. The most relevant effect on the rate constant, as commented above, is given by the use of the 1D hindered rotor model. Finally, though not shown, it has been found that using the Eckart model rather than small curvature theory to compute the tunneling contributions has a negligible impact on the rate constant evaluation. The reason is that tunneling corrections are small for this system because the RW adduct is not significantly collisionally stabilized in the examined temperature and pressure conditions and the energy barrier is submerged with respect to reactants.

Figure 7

Global rate constants computed using different theoretical approaches to determine the fluxes through the inner and outer transition states. The nomenclature is “outer TS–inner TS”: (1) the outer TS flux: VRC-TST (VRC), VTST (VTST), or PST (PST); (2) the inner TS flux: the VTST level with vibrational frequencies computed using internal curvilinear coordinates (VTSTin) or Cartesian coordinates (VTSTcar), conventional TST (cTST), VTST with one internal mode modeled as a 1D hindered rotor (vTST HR).

It can be concluded that the elementary processes that contribute to the reactive fluxes change depending on the ranges of temperatures and pressures that are investigated and that their rate must be determined at a suitable level of theory to obtain quantitative agreement with experimental data. For example, canonical VTST is not apt to study the entrance channel for this system, not much because it assumes a thermal distribution (thus allowing for computing a canonical k(T) rate constant instead of the microcanonical E, or E,J resolved rates, k(E) or k(E,J)), but rather because, as it is often implemented in the literature, it uses the harmonic approximation to evaluate reactive fluxes along the minimum-energy path, which is improper for loosely interacting fragments. A more “proper” evaluation of the reactive fluxes is that given by VRC-TST.

Conclusions

The kinetics of radical–molecule reactions is of remarkable interest in several fields including, inter alia, atmospheric- and astrochemistry. However, obtaining quantitative rate constants for such reactions by means of theoretical methods is challenging because of the difficulties that can be faced in the accurate description of some stationary points (intermediates and/or transition states). Indeed, they might show strong correlation effects, the situation being more involved when third-row atoms are present. Furthermore, SO coupling might be relevant for open-shell species. On these grounds, the first aim of this paper was to investigate a prototypical reaction of this kind, namely, the H2S + Cl addition/elimination reaction, beyond the usual “gold standard” of quantum-chemical calculations, represented by CCSD(T) possibly including the extrapolation to the CBS limit. To this aim, a HEAT-like approach, which includes the full treatment of triple and quadruple excitations together with diagonal Born–Oppenheimer corrections and relativistic effects, combined with a proper treatment of the SO coupling has been employed. This level of theory, in conjunction with anharmonic ZPE corrections evaluated using the double-hybrid revDSD-PBEP86-D3(BJ) functional in the framework of the VPT2 model, is expected to fulfill a sub-kJ mol–1 accuracy, thus allowing an unbiased analysis of the ability of different kinetic models in reproducing the experimental reaction rates. In the present work, different approaches of increasing accuracy have been employed to describe the barrierless entrance channel of the reaction, whose role (in evaluating the global reaction rate) depends on the examined temperature and pressure conditions.

In this connection, even the quite simple PST leads to results within a factor of 2 with respect to their experimental counterparts, whereas the more refined VTST and, especially, VRC-TST models lead to results in quantitative agreement with experiment. These outcomes show unambiguously that this reaction can be well described by models based on the transition state theory, provided that the underlying electronic structure computations are sufficiently accurate and that barrierless channels are properly described.

Furthermore, in view of extending the accuracy reached by our approach to reactive PESs involving larger systems, we have tested the performance of computationally less expensive composite schemes, which would become indeed unavoidable in such cases. Our conclusion is that different variants of the so-called cheap approach perform remarkably well. At the same time, last-generation double-hybrid functionals can be profitably used to optimize geometries and evaluate vibrational contributions. In summary, in our opinion, a promising route for computing reaction rates in semiquantitative agreement with experiment (i.e., well within a factor of 2) for quite large molecular systems can be based on ChS energy evaluations of the relevant stationary points coupled to TST for the activated steps and to PST for barrierless steps. All of these ingredients must be finally introduced in a master equation model of the overall reaction network, which must include all of the elementary processes that can contribute to the reactive fluxes, with rates computed at a suitable level of theory.