Theory Finally Agrees with Experiment for the Dynamics of the Cl + C2H6 Reaction

Since the pioneering reaction dynamics studies of H + H2 in the 1970s, theory increased the system size by one atom in every decade arriving to six-atom reactions in the early 2010s. Here, we take a significant step forward by reporting accurate dynamics simulations for the nine-atom Cl + ethane (C2H6) reaction using a new high-quality spin–orbit–ground-state ab initio potential energy surface. Quasi-classical trajectory simulations on this surface cool the rotational distribution of the HCl product molecules, thereby providing unprecedented agreement with experiment after several previous failed attempts of theory. Unlike Cl + CH4, the Cl + C2H6 reaction is exothermic with an adiabatically submerged transition state, allowing testing of the validity of the Polanyi rules for a negative-barrier reaction.

O ne of the main goals of modern chemistry is to understand how chemical reactions proceed step by step at the atomic and molecular level. Experimentally, the reactions of a chlorine atom (Cl) with organic molecules such as methane (CH 4 ), ethane (C 2 H 6 ), methanol (CH 3 OH), etc. have become benchmark systems to uncover the dynamics of hydrogen-abstraction processes forming hydrogen chloride (HCl) molecules. 1−11 These experiments provided deep insights into the state-to-state dynamics and mode-specific energy transfer of polyatomic reactions, thereby extending and modifying the fundamental rules 12 of chemical reactivity. Furthermore, the experimental findings have also challenged theory to provide explanations, predictions, confirmations, and sometimes contradictions, thereby moving the field forward. 13 −21 Theory of chemical reaction dynamics began with the X + H 2 reactions in the 1970s, 22 followed by initiative studies for X + H 2 O and X + CH 4 in the 1980s and 1990s, 23−25 where X is an atom (H, F, or Cl). The two key steps and challenges of the dynamics simulations are the description of the motion of the electrons and that of the nuclei. The former is done quantum mechanically resulting in a potential energy surface (PES), which governs the latter via classical or quantum methods. Following the pioneering three-dimensional quantum dynamics study of the H + H 2 reaction in 1975, 22 it took about two decades to achieve similar-quality description of the H + H 2 O reaction 24 and another two decades to perform accurate dynamics computations for the X + CH 4 reactions. 13−15,26 For a six-atom system, like X + CH 4 , the PES is a 12-dimensional function, whose accurate representation was a considerable challenge for theory. In 2011, one of us reported a high-quality full-dimensional ab initio PES for the Cl + CH 4 reaction, 13 allowing efficient dynamics simulations using the quasi-classical trajectory (QCT) or reduced-dimensional quantum methods. 14 The key experimental result of the Cl-atom hydrogenabstraction reactions is the rotational distribution of the HCl product. 6−9 In the 2000s several theoretical studies struggled to reproduce the measured data, [17][18][19]27 until 2011, when QCT simulations on the above-mentioned PES finally provided cold rotational distributions for the HCl product of the Cl + CH 4 reaction, 13 in excellent agreement with experiment. 6 For the Cl + C 2 H 6 reaction, experiment also shows the production of rotationally cold HCl, 7−9 whereas reaction dynamics simulations using on-the-fly electronic structure theories 17−19 or valence-bond/molecular-mechanics-type analytical PES 20,21 provide broad rotational distributions, in sharp contradiction with experiment. 8 Of course, this is not surprising, as history shows that two decades of theoretical research are needed to accurately describe 1−2 additional atoms in a chemical reaction. 22,24,26 In the present study, we aim to take a significant step forward by moving from six-atom systems to a nine-atom reaction within a decade using a similar level of theory for Cl + C 2 H 6 as was done in 2011 for Cl + CH 4 . 13 Can the new high-level dynamics cool the HCl rotational distribution for Cl + C 2 H 6 , thereby resolving the contradiction between theory and experiment? Is the reaction of ethane similar to that of methane with Cl? What can we say about the Polanyi rules-predicted 12 translational and vibrational energy effects on reactivity, which were challenged by Cl + CH 4 , 1 in the case of Cl + C 2 H 6 ? Below we answer these questions and describe our theoretical approach from the 21-dimensional PES development to the reaction dynamics simulations.
First we describe the motion of the electrons by numerically solving their Schrodinger equation at fixed nuclear configurations resulting in potential energies. For the Cl + C 2 H 6 system, we have to describe electron correlation and the coupling of the spin and orbital angular momenta, which lowers the reactant asymptote by 0.84 kcal/mol. 28 To achieve the best technically feasible accuracy we use a composite ab initio method based on explicitly correlated perturbation and coupled-cluster methods (MP2-F12 and CCSD(T)-F12b) 29,30 and basis sets up to aug-cc-pVTZ 31 (to obtain accurate correlated wave functions) as well as the interacting-states approach 32 via multireference configuration interaction 33 (to describe spin−orbit effects). The next step is to select the nuclear configurations and represent the potential energies by a 21-dimensional analytical function, which has been a challenge for decades. We recently developed a program system called ROBOSURFER, 34 which automates the construction of analytical PESs by iteratively improving the surface based on (1) selective addition of geometries extracted from classical trajectories, (2) ab initio computations, and (3) fitting the energy points by the monomial symmetrization approach 35 of the permutationally invariant polynomial method. 36 Using ROBOSURFER with the above-described composite ab initio method, we arrive at a high-quality spin−orbit−ground-state PES for the Cl( 2 P 3/2 ) + C 2 H 6 → HCl + C 2 H 5 reaction. With this analytical surface at hand we study the dynamics of the title reaction using the QCT method as the motion of the heavy (relative to the mass of the electron) nuclei can usually be well-described by classical mechanics, especially for the present barrier-less reaction, while the forces are obtained from the gradients of the quantum-mechanics-based PES. The interested reader can consult Computational Methods and the Supporting Information for additional computational details.
The topology of the PES showing the structures and relative energies along the H-abstraction pathway is shown in Figure 1. The classical and zero-point-energy-corrected adiabatic energies relative to Cl( 2 P 3/2 ) + C 2 H 6 can be compared to benchmark data 28 to assess the accuracy of the analytical PES. As seen, all the relative energies obtained on the PES agree with the benchmark values within 0.5 kcal/mol and the PES reaction enthalpy agrees with experiment 37 within only 0.14 kcal/mol, showing the subchemical accuracy of the new PES. In the entrance channel there is a shallow prereaction well, which is 0.57 kcal/mol deep for a reactive Cl···HCH 2 CH 3 configuration ( Figure 1). The PES describes the entrancechannel potential accurately as shown in Figure S1. The ground-state vibrationally adiabatic pathway is exothermic (ΔH 0 = −3.06 kcal/mol), featuring a submerged transition state (TS) and a minimum in the product channel (postmin) below the reactants by 2.12 and 5.19 kcal/mol, respectively. However, if we do not consider zero-point vibration, the classical PES is slightly endothermic (ΔE = 2.04 kcal/mol) and has a positive barrier (2.21 kcal/mol) as shown in Figure 1. In the case of the Cl + CH 4 PES, the classical(adiabatic) barrier is higher, 7.57(3.56) kcal/mol, and the reaction is endothermic, ΔE(ΔH 0 ) = 5.75 (0.73) kcal/mol. 13 At the Cl + C 2 H 6 TS, the C−H and H−Cl distances are stretched by 0.246 and 0.219 Å on the PES, relative to the corresponding bonds in C 2 H 6 and HCl, respectively, with a nearly collinear C−H−Cl arrangement (173°) featuring a slightly late, product-like TS, whereas the Cl + CH 4 TS is exactly collinear and clearly late as the corresponding stretching values are 0.314 and 0.161 Å, respectively. 13 The excellent agreement between the PES and benchmark data at the stationary points does not guarantee the global accuracy of the analytical surface; therefore, the HCl rotational distribution obtained on the new PES provides a critical test to confirm its good performance in dynamics simulations. This test is especially interesting considering the fact that all the previous theoretical attempts 17−21 failed to reproduce the cold experimental results. Considering the experimental−theoretical controversy, one may start to wonder that some quantum effects, which are not captured by QCT, cool HCl rotation. Here, we show that this is not the case, because the QCT simulations on the new PES provide cold HCl rotational distribution, in good and unprecedented agreement with experiment, 8 as shown in Figure 2. Now both the theoretical Figure 1. Schematic representation of the energetics of the Cl( 2 P 3/2 ) + C 2 H 6 → HCl + C 2 H 5 reaction. Classical and adiabatic relative energies obtained on the present PES are compared with benchmark ab initio data 28 (in parentheses) and, for the products, with experimental 0 K reaction enthalpy. 37 For details about premin, see Figure S1. Bond lengths and the C−H−Cl angles obtained on the PES compared to benchmark values 28 19,21 for J = 4−6, while our computations give 6, 3, and 1% for J = 4, 5, and 6, respectively, in excellent quantitative agreement with experiment. Thus, the present findings demonstrate that QCT is capable of correctly describing the dynamics of the Cl + C 2 H 6 reaction, similarly to the case of Cl + CH 4 , 13 if an accurate PES is used, whose postreaction valley plays a key role in hindering the rotation of the departing HCl product, 7 because the C−H−Cl arrangements are nearly collinear at the TS and also at the postmin structures as shown in Figure 1. Seeing the cold HCl rotational distributions and knowing that HCl is mainly formed in its vibrational ground state, the available energy must be transformed into internal excitation of the ethyl radical and relative translation of the products. Our simulations show that the main part of the available energy (∼63%) goes to product translation and about 34−36% channels to ethyl vibration and rotation, whereas the fraction of HCl rotational energy is only 3%, in almost quantitative agreement with the most recent experiment 11 as shown in Table 1. A previous theoretical study provided similar results for ethyl internal and translational energy fractions but overestimated the average HCl rotational energy by a factor of 2. 21 The measured J-resolved translational energy distributions at collision energy of 6.7 kcal/mol cover an energy range from 0 to 10 kcal/mol (translational energy limit was marked at ∼9.5 kcal/mol) and peak around 5−7 kcal/ mol, 11 in excellent agreement with the present computed results ( Figure S4).
Scattering angle distributions are isotropic at low collision energies (Figure 3), indicating indirect dynamics as expected in the case of a barrier-less exothermic reaction with a relatively deep minimum below the product asymptote. As collision energy increases, the angular distributions shift toward forward scattered products, in agreement with experiments of Suits and co-workers. 11 Forward scattering is a signature of direct stripping mechanism occurring at large impact parameters. These findings are in accord with the opacity functions ( Figure  S3), showing a slowly decaying shape at low collision energies and peaking at large impact parameters at high translational energies. The above-described translational energy dependence of the Cl + C 2 H 6 reaction dynamics is similar to Cl + CH 4 , 38 although, because of the more pronounced dispersion interaction between the reactants, the angular distributions of Cl + C 2 H 6 are more isotropic, especially at low collision energies, where Cl + CH 4 is clearly backward scattered. At collision energy of 6.7 kcal/mol J-resolved experimental angular distributions are available for Cl + C 2 H 6 , 11 which indicate a slight preference of forward scattering, in agreement with our computed result shown in Figure 3.
Finally, we investigate the effect of symmetric CH stretching excitation on the dynamics of the Cl + C 2 H 6 reaction, which challenged the rules of thumb of reaction dynamics, formulated by Nobel-laureate John Polanyi, 12 in the case of the reactions of F and Cl with CHD 3 . 1,13,14,39,40 The Polanyi rules predict, based on findings for atom−diatom systems, which form of energy is more efficient to activate a chemical reaction. 12 The rules say that one should invest in translational energy for early-barrier reactions, whereas vibrational energy accelerates late-barrier reactions more efficiently. On one hand, the Cl + C 2 H 6 reaction has a slightly late barrier, thus we expect that CH stretching excitation enhances the reactivity. On the other hand, the Cl + C 2 H 6 reaction has a submerged adiabatic barrier; thus, we do not really need to invest energy to activate the reaction, in contrast to Cl + CH 4 . Figure 4 shows the cross sections of the vibrationally ground-state and CH-stretchingexcited Cl + C 2 H 6 (v 1 = 0,1) reactions as a function of collision energy. The cross sections of the ground-state reaction are   . Integral cross sections for the Cl( 2 P 3/2 ) + C 2 H 6 (v = 0) → HCl + C 2 H 5 and the Cl( 2 P 3/2 ) + C 2 H 6 (v 1 = 1) → HCl + C 2 H 5 reactions as a function of collision energy. The results are obtained on the present PES, and the inset shows the ratio of the symmetric CHstretching-excited (v 1 = 1) and the ground-state (v = 0) cross sections as a function of collision energy.
The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter large (40−60 bohr 2 ), in accord with the substantial reaction probabilities (25−40%, Figure S3), and, apart from the fast decay at low energies, show only modest collision energy dependence (slight decrease with increasing translational energy). This finding is in sharp contrast to Cl + CH 4 , 13 whose excitation function has a threshold and increases with collision energy. The difference between the reactivity of CH 4 and C 2 H 6 with Cl can be explained by the different characteristics of the PESs, as Cl + CH 4 has a positive barrier, 13 whereas Cl + C 2 H 6 is adiabatically barrier-less; thus, additional translational energy slightly hinders the reactivity by allowing less time to find a favorable condition (orientation, vibrational phase, etc.) for reaction. CH stretching excitation increases the reactivity of Cl + C 2 H 6 by about 80% at a low collision energy of 1 kcal/mol and with a decreasing factor as collision energy increases. In the translational energy range of 10−20 kcal/mol, the vibrational enhancement is around 20%. Thus, on one hand, CH stretching excitation enhances the reactivity of the submerged late-barrier Cl + C 2 H 6 reaction more efficiently that the same amount of translational energy, in accord with the Polanyi rules. On the other hand, translational energy actually hinders the reactivity of Cl + C 2 H 6 and thus obviously cannot compete with vibrational excitation even if the stretching effect is modest, as a barrierless exothermic process requires interaction time rather than energy. We conclude that in 2020 one can develop a PES for a nineatom reaction of quality similar to what was possible for sixatom systems about a decade ago. 13 The keys to this achievement are the use of advanced explicitly correlated electronic structure theories, 29,30 a general fitting approach, 35 and automated PES development enabled by the ROBOSURFER program. 34 As a result, we can finally resolve a longstanding contradiction 17−21 between theory and experiment proving the accuracy of our approach for reaction dynamics simulations. We expect that similar dynamics studies of nine-atom reactions will become a new state-of-the-art in the next decade, thereby motivating future experiments.

■ COMPUTATIONAL METHODS
The construction of the PES starts from randomly displaced geometries of the stationary points 28 of the Cl + C 2 H 6 reaction and utilizes the ROBOSURFER program system, 34 which enables the automated development of the surface by the iterative and selective addition of new geometries extracted from QCT simulations, then subjected to ab initio quantum chemical computations, and fitted by using the monomial symmetrization approach 35 of the permutationally invariant polynomial method. 36 The fitting function is a polynomial expansion of the y ij = exp(−r ij /a) Morse-like variables of the r ij interatomic distances, where a = 1.5 bohr is applied and the highest total order of the polynomials is 5. The energy points are fitted using a least-squares fit with an E 0 /(E + E 0 ) weighting factor, where E is the actual energy relative to the global minimum of the fitting set and E 0 = 0.04 hartree. The fifth-order expansion requires 3234 fitting coefficients. In the first round of the PES development we apply the RMP2/aug-cc-pVDZ level of theory for quantum chemical computations and perform 368 ROBOSURFER iterations. The final energy points are then recomputed at the following composite level of theory: where "SO corr " is the spin−orbit (SO) correction to each energy point and refers to the difference between the relativistic spin−orbit and nonrelativistic ground state of the system. The SO correction is determined via the interactingstates approach 32  Cross sections are calculated by a b-weighted numerical integration of the P(b) opacity functions at each collision energy. To obtain the rotational state distribution of the HCl(v = 0) product molecules, only those trajectories are used where the vibrational energy of the C 2 H 5 product is greater than its zero-point energy (ZPE) and the internal energy of the HCl product is greater than its ZPE corresponding to its actual rotational state. Additional