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.

One 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 (CH4), ethane (C2H6), methanol (CH3OH), 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 + H2 reactions in the 1970s,[22] followed by initiative studies for X + H2O and X + CH4 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 + H2 reaction in 1975,[22] it took about two decades to achieve similar-quality description of the H + H2O reaction[24] and another two decades to perform accurate dynamics computations for the X + CH4 reactions.[13−15,26] For a six-atom system, like X + CH4, 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 + CH4 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 hydrogen-abstraction reactions is the rotational distribution of the HCl product.[6−9] In the 2000s several theoretical studies struggled to reproduce the measured data,[17−19,27] until 2011, when QCT simulations on the above-mentioned PES finally provided cold rotational distributions for the HCl product of the Cl + CH4 reaction,[13] in excellent agreement with experiment.[6] For the Cl + C2H6 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 + C2H6 as was done in 2011 for Cl + CH4.[13] Can the new high-level dynamics cool the HCl rotational distribution for Cl + C2H6, 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 + CH4,[1] in the case of Cl + C2H6? 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 Schrödinger equation at fixed nuclear configurations resulting in potential energies. For the Cl + C2H6 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(2P3/2) + C2H6 → HCl + C2H5 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 . The classical and zero-point-energy-corrected adiabatic energies relative to Cl(2P3/2) + C2H6 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···HCH2CH3 configuration (Figure ). The PES describes the entrance-channel potential accurately as shown in Figure S1. The ground-state vibrationally adiabatic pathway is exothermic (ΔH0 = −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 . In the case of the Cl + CH4 PES, the classical(adiabatic) barrier is higher, 7.57(3.56) kcal/mol, and the reaction is endothermic, ΔE(ΔH0) = 5.75 (0.73) kcal/mol.[13] At the Cl + C2H6 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 C2H6 and HCl, respectively, with a nearly collinear C–H–Cl arrangement (173°) featuring a slightly late, product-like TS, whereas the Cl + CH4 TS is exactly collinear and clearly late as the corresponding stretching values are 0.314 and 0.161 Å, respectively.[13]

Figure 1

Schematic representation of the energetics of the Cl(2P3/2) + C2H6 → HCl + C2H5 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] (in parentheses) of the transition state (TS) and the post-TS minimum (postmin) geometries are given in angstroms and degrees, respectively.

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 . Now both the theoretical and experimental results peak at J = 1 and decay sharply with increasing J, vanishing at J = 6, whereas previous simulations[17−21] gave broad distributions between J values of 0 and 10, peaking around 2–4 and having 15–20% populations[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 + C2H6 reaction, similarly to the case of Cl + CH4,[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 .

Figure 2

Rotational state distribution of the HCl product at 5.5 kcal/mol collision energy of the Cl(2P3/2) + C2H6(v = 0) → HCl(v = 0, JHCl) + C2H5 reaction. The theoretical results, obtained on the present PES, are compared to experimental data taken from ref (8).

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 . 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).

Table 1

Fractions (%) of the Available Energy Corresponding to the Rotational Motion of the HCl Molecule (frot), the Internal Excitation of the Ethyl Radical (fint), and the Relative Translation of the Products (ftrans), Compared to Experimental Results[11] at Different Collision Energies (Ecoll)

	Ecoll = 5.5 kcal/mola		Ecoll = 6.7 kcal/mol
	theoryb	experimentc	theoryb	experimentc
frot(HCl)	3	3	3	2
fint(C2H5)	34	35	36	30
ftrans	63	63	62	68

The experimental Ecoll is 5.2 kcal/mol.

QCT results obtained on the present PES.

See Table 1 of ref (21).

Scattering angle distributions are isotropic at low collision energies (Figure ), 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 + C2H6 reaction dynamics is similar to Cl + CH4,[38] although, because of the more pronounced dispersion interaction between the reactants, the angular distributions of Cl + C2H6 are more isotropic, especially at low collision energies, where Cl + CH4 is clearly backward scattered. At collision energy of 6.7 kcal/mol J-resolved experimental angular distributions are available for Cl + C2H6,[11] which indicate a slight preference of forward scattering, in agreement with our computed result shown in Figure .

Figure 3

Normalized scattering angle distributions for the Cl(2P3/2) + C2H6(v = 0) → HCl + C2H5 reaction. The results are obtained on the present PES at different collision energies (Ecoll given in kcal/mol), where cos(θ) = −1 corresponds to backward scattering.

Finally, we investigate the effect of symmetric CH stretching excitation on the dynamics of the Cl + C2H6 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 CHD3.[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 + C2H6 reaction has a slightly late barrier, thus we expect that CH stretching excitation enhances the reactivity. On the other hand, the Cl + C2H6 reaction has a submerged adiabatic barrier; thus, we do not really need to invest energy to activate the reaction, in contrast to Cl + CH4. Figure shows the cross sections of the vibrationally ground-state and CH-stretching-excited Cl + C2H6(v1 = 0,1) reactions as a function of collision energy. The cross sections of the ground-state reaction are 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 + CH4,[13] whose excitation function has a threshold and increases with collision energy. The difference between the reactivity of CH4 and C2H6 with Cl can be explained by the different characteristics of the PESs, as Cl + CH4 has a positive barrier,[13] whereas Cl + C2H6 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 + C2H6 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 + C2H6 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 + C2H6 and thus obviously cannot compete with vibrational excitation even if the stretching effect is modest, as a barrier-less exothermic process requires interaction time rather than energy.

Figure 4

Integral cross sections for the Cl(2P3/2) + C2H6(v = 0) → HCl + C2H5 and the Cl(2P3/2) + C2H6(v1 = 1) → HCl + C2H5 reactions as a function of collision energy. The results are obtained on the present PES, and the inset shows the ratio of the symmetric CH-stretching-excited (v1 = 1) and the ground-state (v = 0) cross sections as a function of collision energy.

We conclude that in 2020 one can develop a PES for a nine-atom reaction of quality similar to what was possible for six-atom 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 + C2H6 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 = exp(−r/a) Morse-like variables of the r 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 E0/(E + E0) weighting factor, where E is the actual energy relative to the global minimum of the fitting set and E0= 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 “SOcorr” 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 interacting-states approach[32] at the MRCI+Q/aug-cc-pVDZ[31,33] level of theory for each geometry. The second (and final) round of the PES development uses the above-defined composite ab initio level of theory and consists of 140 Robosurfer iterations. The final PES is built from 11 701 geometries and the corresponding composite energies.

QCT simulations are run at seven collision energies: 1.0, 3.0, 5.5, 6.7, 10.0, 15.0, and 20.0 kcal/mol for the ground-state (v = 0) and symmetric CH-stretching-excited (v1 = 1) Cl + C2H6(v = 0, v1 = 1) reactions. The initial vibrational state of ethane is set by standard normal-mode sampling,[41] and the spatial orientation of the reactants is randomly sampled. The initial distance of the Cl atom and the center of mass of the ethane molecule is , where x = 16 bohr and the b impact parameter is varied between 0 and bmax with a step size of 0.5 bohr. 1000 trajectories are run at each b value.

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 C2H5 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 computational details are described in the Supporting Information.