Detailed benchmark ab initio mapping of the potential energy surfaces of the X + C 2 H 6 [X = F, Cl, Br, I] reactions

We investigate three reaction pathways of the X + C 2 H 6 [X = F, Cl, Br, I] reactions: H-abstraction, methyl-substitution, and H-substitution, with the latter two proceeding via either Walden-inversion or front-side-attack mechanisms. We report classical and adiabatic relative energies of unprecedented accuracy for the corresponding stationary points of the reaction potential energy surfaces (PESs) by augmenting the CCSD(T)-F12b/aug-cc-pVQZ energies by core-correlation, post-CCSD(T) and spin–orbit corrections. Taking these correction terms into account turns out to be essential to reach subchemical, i.e. o 0.5 kcal mol (cid:2) 1 , accuracy. Our new benchmark 0 K reaction enthalpies show excellent agreement with experimental data. Spin–orbit coupling in these open-shell systems is also monitored throughout the reaction paths and found to be non-negligible even in some transition-state geometries. Barrier heights corresponding to the different channels of the title reactions appear in the same order with increasing energy: H-abstraction, Walden-inversion methyl-substitution, Walden-inversion H-substitution, front-side-attack H-substitution and front-side-attack methyl-substitution, except for X = I where the latter two come in reverse order. Similarly, product channels follow the energy order of the corresponding barrier heights in all four cases. We find strongly reactant-like transition-state structures for the exothermic F + C 2 H 6 reaction paths, while more and more product-like transition states are observed along with increasing endothermicity as going from Cl to I. Several entrance and exit channel minima are also identified for the studied reactions with significant spin–orbit effects for the formers.


Introduction
Reactions of atoms (X) with methane [1][2][3][4][5][6] and ions (X À ) with methyl-halides (CH 3 Y) [7][8][9][10] have become benchmark systems to study the dynamics and mechanisms of polyatomic chemical reactions.The main reaction channel of the former systems is hydrogen-abstraction leading to HX + CH 3 , and for the latter bimolecular nucleophilic substitution (S N 2) resulting in Y À + CH 3 X.Recent studies, however, revealed that these reactions are not so simple and several other reaction channels may become available at higher collision energies. 3,7,8,11For X + methane reactions substitution leads to H + CH 3 X products via different stereospecific reaction pathways involving configuration inversion or retention around the carbon center. 11In the case of the X À + CH 3 Y systems the S N 2 reaction can proceed with Walden-inversion, front-side-attack, double-inversion, 12 roundabout, 10 etc. mechanisms and proton transfer leading to HX + CH 2 Y À may occur at higher collision energies.Furthermore, as we replace CH 3 Y with CH 3 CH 2 Y a new channel, called bimolecular elimination (E2), leading to Y À + HX + C 2 H 4 becomes possible.Recently we reported a benchmark ab initio characterization of the F À + CH 3 CH 2 Cl reaction. 13In the present work we plan to extend our previous work 1,11 on X + methane reactions by replacing methane with the ethane (C 2 H 6 ) molecule.Halogen + ethane reactions represent a more complex system with more of degrees of freedom and new possible reaction mechanisms.Besides the main H-abstraction reaction path and the higher-energy H-substitution studied also for halogen + methane reactions, in the X + C 2 H 6 [F, Cl, Br, I] reactive system methyl-substitution appears as a new mechanism and is expected to occur at lower energies then H-substitution.For these reactions so far only the H-abstraction channel (X = F, Cl) [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33] or its reverse reaction (X = Br, I) [34][35][36][37][38][39][40][41] have been investigated.
Stationary points for the F + C 2 H 6 -HF + C 2 H 5 reaction were identified in 2007 using the CCSD(T)/cc-pVDZ method along with energy computations at the CCSD(T)/cc-pVQZ level of theory, 14 proposing an early-barrier transition state (TS) for the highly exothermic reaction.Later state-to-state scattering experiments also proposed an early-barrier TS from the mostly vibrationally excited HF products, however, a significant amount of vibrationally non-excited HF molecules implied energy transfer from HF to the ethyl unit during the H-abstraction reaction. 15irect dynamics simulations were performed by Troya and coworkers using a reparametrized semiempirical Hamiltonian, based on MP2/aug-cc-pVDZ geometry optimization and CCSD(T)/ aug-cc-pVDZ energy computation on the stationary points, for the reactions of fluorine with small alkanes including ethane, where the effect of the length of the alkane chain on the dynamics was studied. 16Espinosa-Garcia et al. have very recently reported the first analytical valence-bond -force-field-based PES of the F + C 2 H 6 -HF + C 2 H 5 reaction based on CCSD(T)/triple-zeta ab initio data and depending on 63 adjustable parameters. 17he hydrogen-abstraction reaction between Cl and ethane has been studied more extensively both experimentally and theoretically.In the late '90s Zare and coworkers investigated this reaction by the state-selective detection of the HCl 18,19 or the ethyl radical product. 20They proposed a loose transition state and claimed that the ethyl unit was a spectator during the slightly exothermic reaction.Later crossed molecular beam experiments detecting HCl product states indicated significant internal excitation of the ethyl radical, which contradicts its spectator character. 21,23,24Crossed-beam scattering experiments of Suits and coworkers also strengthened the considerable ethyl internal excitation. 31,32In parallel, stationary points of the Cl + C 2 H 6 -HCl + C 2 H 5 reaction were theoretically characterized using the MP2 method with correlation-consistent basis sets combined with CCSD(T) energy computations suggesting a loose TS and a post-TS minimum. 25The reaction was also studied by direct dynamics simulations using semiempirical Hamiltonians. 22,257][28][29][30] In a recent paper Rangel and Espinosa-Garcia have published an analytical PES including valence-bond and molecular-mechanic terms for this abstraction reaction. 33They also characterized the stationary points using the CCSD(T) method with different basis sets and implied that high-level quantum chemical computations are needed to determine the exact energetics of the reaction.
5][36][37][38][39][40][41] For the HBr + C 2 H 5 -Br + C 2 H 6 reaction theoretical studies using the CCSD(T)/cc-pVTZ method as highest level of theory revealed a loosely-bound adduct in the entrance channel and a tightly-bound TS very close in energy and both being below the reactant asymptote. 35,38,39The HI + C 2 H 5 -I + C 2 H 6 reaction has also been studied both by experimental 34,40 and theoretical 40,41 techniques.Leplat et al.   reported CCSD/cc-pVDZ structures along with spin-orbit corrections for the stationary points of this reaction. 41n the present study we extend the investigations of the halogen + ethane reactions towards new, so far unstudied mechanisms, such as H-substitution and methyl-substitution.We also consider the Walden-inversion and front-side-attack pathways for both substitution channels.We identify the corresponding stationary points using the explicitly-correlated CCSD(T)-F12b/augcc-pVTZ level of theory, and provide the corresponding classical and adiabatic relative energies with unprecedented accuracy, i.e.CCSD(T)-F12b/aug-cc-pVQZ energies corrected by core-valence, post-CCSD(T), and spin-orbit effects.Spin-orbit coupling throughout the reaction paths along with structural changes are also thoroughly analyzed.This work is meant to be a first step towards developing global analytical ab initio potential energy surfaces thereby making possible to study the multi-channel dynamics of the title reactions in detail.

Computational details
In this work three different channels of the X + C 2 H 6 [F, Cl, Br, I] reactions are studied: hydrogen-abstraction (HA), hydrogensubstitution (HS), and methyl-substitution (MS) leading to C 2 H 5 + HX, C 2 H 5 X + H, and CH 3 X + CH 3 , respectively.HS and MS can proceed via a Walden-inversion (W) transition state (TS) or via a front-side-attack (FS) TS.This nomenclature is used throughout the paper.In the case of HA post-TS minima, referred to as HA postmin, are also found.
The stationary-point geometries of the potential energy surfaces of the X + C 2 H 6 [F, Cl, Br, I] reactions are pre-optimized using the restricted open-shell second-order Møller-Plesset perturbation theory 42 (RMP2) using the correlation-consistent augcc-pVDZ basis set. 43A few geometries (HS FS TS for F + C 2 H 6 , and HA TS and postmin for Br + C 2 H 6 ) could only be identified with the unrestricted MP2 (UMP2) method, in these cases we use the UMP2 energies of the reactant halogen atom to determine relative energies.Then the restricted open-shell Hartree-Fock-based unrestricted explicitly-correlated coupled cluster singles, doubles, and perturbative triples method, ROHF-UCCSD(T)-F12b, 44 with the aug-cc-pVDZ basis set is applied for geometry optimization and harmonic frequency computations.Geometry optimizations are then also performed at the ROHF-UCCSD(T)-F12b/aug-cc-pVTZ 43 level of theory for all stationary points.To obtain benchmark relative energies further single-point energy computations are carried out at the most accurate ROHF-UCCSD(T)-F12b/augcc-pVTZ geometries: ROHF-UCCSD(T)-F12b/aug-cc-pVQZ, 43 UHF-UCCSD(T), 45 UHF-UCCSDT, 46 and UHF-UCCSDT(Q) 47 with the cc-pVDZ basis set, and ROHF-UCCSD(T)/aug-cc-pwCVTZ 48 with both frozen-core (FC) and all-electron (AE) approaches.For the Br and I atoms small-core relativistic effective core potentials (ECPs) 49 are used with the corresponding cc-pVDZ-PP, aug-cc-pVDZ-PP, aug-cc-pVTZ-PP, aug-cc-pVQZ-PP, and aug-cc-pwCVTZ-PP basis sets. 50ECPs replace the inner-core 1s 2 2s 2 2p 6 electrons of Br and the 1s 2 2s 2 2p 6 3s 2 3p 6 3d 10 electrons of I.In all cases the frozen-core approach is applied correlating the valence This journal is © the Owner Societies 2019 electrons only, except in AE computations, where the following electrons are also correlated: 1s 2 for C and F, 2s 2 2p 6 for Cl, 3s 2 3p 6 3d 10 for Br, and 4s 2 4p 6 4d 10 for I.
For spin-orbit (SO) computations a corrected multireference configuration interaction [MRCI+Q] method 51,52 is applied, within the frozen-core approach.The Davidson correction (+Q) 53 estimates higher-order correlation energy effects.In the multireference computations a minimal active space of 5 electrons in the 3 spatial nplike orbitals, where n = 2, 3, 4, and 5 for F, Cl, Br, and I, respectively, is used.The SO computations are based on the Breit-Pauli operator in the interacting-states approach. 54The SO eigenstates are obtained by diagonalizing a 6 Â 6 SO matrix, whose diagonal elements are replaced by the Davidson-corrected MRCI energies.SO computations have been carried out scanning the entrance channel of the X + C 2 H 6 [F, Cl, Br, I] reactions, representing three different directions of the halogen atom approaching the ethane molecule: (1) X approaching ethane from one of the methyl-groups maintaining a C 3v point-group symmetry, (2) X approaching one H atom of ethane (3) X approaching perpendicularly the C-C bond of ethane.The ethane molecule is kept fixed at its CCSD(T)-F12b/ aug-cc-pVTZ equilibrium geometry while changing the separation distance of X. SO computations are also performed at the CCSD(T)-F12b/aug-cc-pVTZ geometries of the TSs, postmins and products, where the products are separated by 15 Å.The DSO corrections to the TS, postmin and product relative energies are determined by subtracting the difference of the computed SO and non-SO groundstate energies from the 1/3 of the experimental SO splittings (e) of the halogen atoms, which splittings are the following: e = 1.16(F), 2.52(Cl), 10.54(Br), and 21.74(I) kcal mol À1 .
All ab initio computations are carried out with the Molpro program package, 55 except the UHF-UCCSD(T), UHF-UCCSDT, and UHF-UCCSDT(Q) computations, which are performed using the MRCC program 56 interfaced to Molpro.
The benchmark classical relative energies at the UCCSD(T)-F12b/aug-cc-pVTZ geometries of the stationary points are calculated as follows: where and The benchmark adiabatic relative energies are determined as follows: where DZPE is the zero-point-energy correction obtained from the harmonic frequency computations performed at the UCCSD(T)-F12b/aug-cc-pVDZ level of theory.
In determining the relative energies of the stationary points a few issues have emerged: (1) the HA TS of the F + C 2 H 6 reaction could not be identified by optimizing the geometry at the UCCSD(T)-F12b/aug-cc-pVDZ level of theory, probably due to the flatness of the entrance channel of the corresponding potential energy surface and/or the large multireference character of this TS, thus, instead, an UCCSD(T)-F12b/aug-cc-pVDZ single point energy computation is performed on the RMP2/ aug-cc-pVDZ geometry and for ZPE correction the RMP2/aug-cc-pVDZ harmonic frequencies are used.Due probably also to the larger multireference character, at the UCCSD(T)-F12b/aug-cc-pVTZ geometry of this HA TS neither the RHF, nor the UHF computation have converged during the correction-term computations.(2) Since in the I + C 2 H 6 HA reaction the TS and postmin geometries could only be identified at the CCSD/ cc-pVDZ-PP level of theory, reproducing previous results, 41 UCCSD(T)-F12b/aug-cc-pVDZ-PP, UCCSD(T)-F12b/aug-cc-pVTZ-PP single point energies and all the above defined correction terms are computed for these structures (and also for ethane in these cases) at the UCCSD/cc-pVDZ-PP geometries.To calculate the ZPE correction for these TS and postmin the CCSD/cc-pVDZ-PP harmonic frequencies are used.Attempts to optimize these geometries have been made at the following levels of theory: RMP2/aug-cc-pVDZ-PP, UMP2/aug-cc-pVDZ-PP, UCCSD-F12b/aug-cc-pVDZ-PP, UCCSD(T)-F12b/aug-cc-pVDZ-PP, UCCSD(T)-F12b/aug-cc-pVTZ-PP, using different geometry optimization algorithms available in Molpro.

Results and discussion
The schematic representation of the potential energy surfaces (PESs) and different pathways of the X + C 2 H 6 [F, Cl, Br, I] reactions are shown in Fig. 1, with the benchmark classical (adiabatic) energies of the stationary points relative to the reactants, obtained from eqn (1) (eqn ( 5)).The relative energies determined at different ab initio levels of theory (RMP2/aug-cc-pVDZ and UCCSD(T)-F12b with aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ bases) and the post-CCSD(T), core-correlation, spin-orbit and ZPE corrections are presented in Table 1.
For all four X + C 2 H 6 reactions the 0 K reaction enthalpies of the different reaction channels follow the same order: H-abstraction (HA), H-substitution (HS), and methyl-substitution (MS) with increasing relative energies.The energy difference between the products of HS and MS is around 10 kcal mol À1 in all cases, whereas, as the atomic number of the halogen increases the energy gap between the MS and the HA product channels decrease significantly.The different reaction paths become more and more endothermic from X = F to X = I, along with increasing barrier heights, which usually also follow the same energy order regarding all four reactions.The lowest-barriers appear for HA, the main reaction path at low collision energies, followed by the higher-energy TSs of the Walden-inversion pathway of MS and then that of HS.At high energies first the front-side-attack (FS) mechanism of HS becomes available and finally the FS MS channel features higher barriers, except for X = I, where FS HS has higher barrier than FS MS.Thus, with increasing collision energy, first the HA, then the MS and finally the HS channels open.For X = F and Cl the global minima of the PESs correspond to the HA postmin geometries, whereas for X = Br and I the lowest energy stationary points are located in the entrance channel, see discussion below.
In the F + C 2 H 6 reaction all channels are exothermic with 0 K reaction enthalpies À36.3, À21.1 and À12.6 kcal mol À1 , for HA, MS, and HS, respectively.For the Cl + C 2 H 6 reaction only HA has negative 0 K reaction enthalpy (À3.0 kcal mol À1 ), however, if ZPE correction is not taken into account, HA is an endothermic (2.3 kcal mol À1 ) reaction path.For X = Br and I all reaction channels are endothermic, even the main HA channel has 10.8 and 26.4 kcal mol À1 barrier heights, respectively, with the highest channel-opening at 71.7 kcal mol À1 for the front-side I + C 2 H 6 -C 2 H 5 I + H substitution reaction.Endothermicity and the corresponding barrier heights increase significantly, with 16-33 kcal mol À1 by changing F with Cl (except for HA TS), and somewhat more moderately (with 8-17 kcal mol À1 ) from X = Cl to Br, or from X = Br to I. Regarding H-abstraction, between the TS and postmin geometries an energy gain of 38.1 kcal mol À1 occurs for X = F, which decreases to 2.8 kcal mol À1 for X = Cl, and almost vanishes in the case of X = Br and I, where the close relative energies of the TS and the postmin indicate a flat potential energy surface region and make hard to identify these stationary points.(Note that the relative energies at the optimized postmin geometries without corrections are lower than at the optimized TS geometries both for X = Br and I).H-abstraction is a barrierless channel for the F + C 2 H 6 reaction and also for  1, obtained from eqn (1) (eqn ( 5)).
This journal is © the Owner Societies 2019 Cl + C 2 H 6 when the adiabatic energies are considered.Interestingly, the reverse HA reactions for Br and I, i.e. the widelystudied [34][35][36][37][38][39][40][41] HX + C 2 H 5 [X = Br, I] reactions also have submerged barriers, thus can occur spontaneously.The higher-energy HS and MS reactions proceed either via Walden or front-side TSs, where the FS mechanism has about 11-17 kcal mol À1 (17-23 kcal mol À1 ) higher barriers than Walden inversion in the case of H-substitution (methyl-substitution) for all four reactions.Two HS FS TSs are found close in energy (with 0.5-2.8kcal mol À1 differences) in all four cases, one with C s and the other with C 1 (having always the lower relative energy) point-group symmetry.
All ZPE corrections to the relative energies, as shown in Fig. 2, turn out to be negative covering the 0.5-6 kcal mol À1 energy range, so these terms always lower the barrier heights and reaction energies.The smallest corrections (excluding the very reactant-like HA TS of the F + C 2 H 6 reaction) correspond to  methyl-substitution TSs, and product channels, where the breaking C-C bonds are replaced by newly forming C-X bonds, resulting in similar harmonic vibrational frequencies, however, the number of vibrational degrees of freedom (dofs) decreases by 3 during product formation, which causes the negative ZPE correction.Larger corrections appear for the HS TSs and products, due to the formation of C-X bonds instead of the breaking C-H bonds, where C-X bonds have smaller vibrational frequencies.Slightly larger ZPE corrections are observed for the HA TSs and products, where the C-H bonds are replaced by H-X bonds, resulting in smaller vibrational frequencies along with the reduction of dofs by 2. As the C-X and H-X frequencies decrease from X = F to I, in Fig. 2 an increase in the magnitude of ZPE corrections can be observed for the corresponding stationary points from left to right.
In Table 1 the MP2/aug-cc-pVDZ (MP2/DZ) relative energies and ZPE corrections are also listed.This level of theory provides relative energies far from chemical accuracy, with average deviations of 2.8 kcal mol À1 (with a maximum value of 5.4 kcal mol À1 ) from the CCSD(T)-F12b/aug-cc-pVTZ energies.MP2/DZ energies are even poorer when compared to the benchmark classical values, with an average difference of 3.2 kcal mol À1 and a maximum deviation of 6.9 kcal mol À1 .The harmonic ZPE corrections obtained at the MP2/DZ level are of better quality, reproducing the CCSD(T)-F12b/aug-cc-pVDZ ZPE contributions with an average deviation of 0.37 kcal mol À1 (the maximum difference is 1.62 kcal mol À1 ).Therefore, despite being widely used, the MP2/aug-cc-pVDZ level of theory is not recommended for studying such systems if chemical accuracy is desired.
The different correction terms summarized in Table 1 and presented graphically in Fig. 3 and 4 allow estimating the accuracy of the new benchmark relative energies.The remarkably fast basis set convergence of the UCCSD(T)-F12b relative energies is highlighted in Fig. 3.The aug-cc-pVTZ (TZ) and the aug-cc-pVQZ (QZ) energies are of similar quality, with an average deviation of only 0.1 kcal mol À1 , and with the largest deviations observed for the X = I TSs.The TZ-DZ increments show an average magnitude of 0.6 kcal mol À1 , which, together with the small QZ-TZ corrections, indicates that even the DZ energies are within chemical accuracy for most of the stationary points.
As shown in Table 1 and Fig. 4 the post-CCSD(T) and corevalence corrections are usually in the same order of magnitude or even larger than the basis set convergence errors, especially for X = I.The core-correlation terms have an average value of 0.4 kcal mol À1 , and are usually positive and relatively small (average: 0.2 kcal mol À1 ) for X = F and Cl.For X = Br and I, this effect becomes more significant (average: 0.7 kcal mol À1 ) and always lowers the barrier heights and reaction energies.For the C 2 H 5 I + H products this correction reaches even À1.3 kcal mol À1 .d[CCSDT] and d[CCSDT(Q)] increments also represent a significant contribution to the relative energies and cannot be neglected  This journal is © the Owner Societies 2019 if chemical accuracy is desired.][59] However, unexpectedly for X = Br and I both the core-valence and post-CCSD(T) corrections are negative, thus the cancelling effect does not occur.In such open-shell systems the energy lowering due to the relativistic spin-orbit (SO) interaction in the halogen atoms should also be taken into account.The non-relativistic ground state of a halogen atom is 2 P, which splits into two energy levels, 2 P J=3/2 and 2 P J=1/2 , due to spin-orbit coupling, where J is the total angular momentum quantum number.Because of the (2J + 1)-fold degeneracy of the SO states, the ground 2 P 3/2 state is lower in energy by 1/3 of the SO-splitting, while the excited 2 P 1/2 state is higher by 2/3 of the splitting energy with respect to the non-SO 2 P state.In Fig. 5, where all the data are shown relative to the ground non-SO 2 P asymptotic limit, we monitored what happens if the isolated halogen atom approaches the ethane molecule.As seen in Fig. 5 the interaction with the ethane molecule splits the fourfold degenerate 2 P 3/2 state into a twofold degenerate SO ground state (SO 1 ) and an (also twofold degenerate) excited (SO 2 ) state.The higher-energy 2 P 1/2 state (SO 3 ) is not further split.The sixfold degenerate non-relativistic 2 P state is also split due to the interaction into a non-SO ground state (non-SO 1 ), and two non-SO excited states (non-SO 2 and non-SO 3 ), all three states having twofold degeneracy.Among these only the SO 1 and non-SO 1 ground states are reactive within the Born-Oppenheimer approximation and therefore of interest concerning the title reactions.We consider three different directions of the halogen atom (X) nearing to the C 2 H 6 unit: (1) X approaching one of the methyl-groups of ethane maintaining a C 3v point-group symmetry, (2) X approaching one H atom of ethane, and (3) X approaching perpendicularly the C-C bond of ethane.In the first case, the C 3v symmetry causes an exact degeneracy of the non-SO 2 and non-SO 3 states, as can be seen in the first column of Fig. 5.In our MRCI+(Q)(5,3)/augcc-pVDZ(-PP) (PPs for Br and I) computations the experimental atomic splittings, which are 1.16(F), 2.52(Cl), 10.54(Br), and 21.74(I) kcal mol À1 , are reproduced qualitatively well, within a few 0.1 kcal mol À1 , calculated as the difference between the asymptotic limits of the SO potential energy curves of Fig. 5.As seen from the energy scales of Fig. 5 and the experimental splitting values, spin-orbit coupling has a significant energylowering effect in the entrance channel of the title reactions, especially when the reactants are far from each other.This effect is around a few 0.1 kcal mol À1 for X = F, almost 1 kcal mol À1 for X = Cl, a few kcal mol À1 for Br, and can exceed 7 kcal mol À1 for X = I.The reactive non-SO 1 and SO 1 potentials feature van der Waals (vdW) minima along the CÁ Á ÁX coordinate, that is entrance-channel valleys of the PESs of the X + C 2 H 6 reactions, which may play an important role in their dynamics.In the case of X = Br and I these minima turn out to be the global minima of the corresponding PESs.
To visualize the relative depths of these vdW valleys, in Fig. 6 scaled non-SO 1 and SO 1 potential energy curves are shown with the SO 1 curves shifted by their asymptotic limits and all data being relative to the non-SO 1 potential energy asymptote of the corresponding halogen.From X = F to I the entrance channel minima become deeper and deeper, and the energy gap between the non-SO 1 and the less deep SO 1 minima increases.These minima are the deepest in the C 3v -symmetry case, when the halogen approaches one methyl-group (0.5-1.5 kcal mol À1 ), and are slightly shallower when the halogen atom comes perpendicularly (0.5-1.3 kcal mol À1 ).The smallest depths, 0.3-1.1 kcal mol À1 , appear when the halogen accesses one of the H-atoms, similar to the X + CH 4 systems. 1 It is also seen that the CÁ Á ÁX distance corresponding to the minima of both the non-SO 1 and the SO 1 potentials is shifted towards larger values as one goes from F to I, and the SO 1 minima always appear at larger CÁ Á ÁX distances than the non-SO 1 ones.The CÁ Á ÁX distances corresponding to the minimum values of the potentials are the smallest (3.0-3.7 Å) in the C 3v -symmetry case, whereas in the other two cases they are considerably larger (3.5-4.3Å).Such entrance-channel minima may play a crucial role in the reactivity of the low-barrier HA pathways of the F and Cl + C 2 H 6 reactions.SO-coupling is usually considered to play a role only in the entrance channel of the reactions of halogen atoms, where the halogen is ''far'' from its reaction partner.Nonetheless, we carried out SO computations at the TS and even at the product geometries to check the validity of this assumption.It turned out that if the TS is reactant-like (which is mostly the case for X = F, see discussion below) or if the SO splitting is large (thus especially for X = I), the lowering of relative energies due to SO coupling at the TS geometries can reach even the 20% of the SO corrections of the corresponding halogen atoms, as shown in Fig. 7.For X = F, the HA TS has a strongly reactant-like structure with the F atom being 1.94 Å from the H atom of ethane, which gives rise to considerable SO-effect (0.06 kcal mol À1 , 15% of the SO correction of F).In the case of X = Cl SO-coupling is not very significant (0.01-0.02 kcal mol À1 ) for the TSs, while for X = Br even a 0.3 kcal mol À1 (9% for MS W TS) energy Fig. 5 Potential energy curves obtained at the MRCI+Q(5,3)/aug-cc-pVDZ(-PP) (PPs for Br and I) level as a function of the C 2 H 6 Á Á ÁX separation, for X = F (1st row), Cl (2nd row), Br (3rd row), I (4th row).We consider three different separation directions: (left column) X approaching ethane from one of the methyl-groups maintaining a C 3v point-group symmetry, (middle column) X approaching one H atom of one of the methyl groups, and (right column) X approaching perpendicularly the C-C bond of ethane, see inserted geometries.The C 2 H 6 unit is kept frozen at its CCSD(T)-F12b/aug-cc-pVTZ equilibrium structure while changing the separation distance of X. SO 1 and non-SO 1 (red curves) denote the spin-orbit (SO) and non-SO ground electronic states, respectively, whereas the yellow and purple curves correspond to first and second excited states, respectively.lowering can be observed.The largest effect, both in the change in relative energy and in percent, emerges for X = I, where most of the TSs feature 15-20% SO corrections compared to that of the I atom, with a maximum value of 1.4 kcal mol À1 in the case of MS FS TS.It is apparent that if subchemical accuracy is aimed, the SO-coupling cannot be neglected in determining the relative energies of the TSs of the title reactions.In the case of the product geometries SO coupling is not significant, however, for X = I it has a À0.02 and À0.04 kcal mol À1 contribution to the MS and HA reaction energies.
Scalar relativistic effects are taken into account for Br and I through the effective core potentials used in our computations.For all other atoms the scalar relativistic corrections are expected to be small (10-50% of the core correlation effects), 59 thus these are neglected in this study.
The above determined corrections make possible to estimate the uncertainty of our computed benchmark relative energies of the stationary points of the X + C 2 H 6 [X = F, Cl, Br, I] reactions.Due to the fast basis-set convergence of the explicitly correlated CCSD(T)-F12b method, the basis-set error might only contribute to the overall uncertainty with about 0.1 kcal mol À1 .Somewhat larger uncertainties may come from the post-CCSD(T) corrections, whose error could reach 50% of the correction value, giving rise   uncertainty.Taken together, the total uncertainty of our computed data is estimated to be within 0.5 kcal mol À1 , far superseding chemical accuracy.
In Table 2 we compare the benchmark 0 K reaction enthalpies determined in this work to those obtained from the Active Thermochemical Tables (ATcT) along with their uncertainties.A very good agreement is observed as the differences of the best available experimental and our new theoretical data are between 0.05-0.65 kcal mol À1 regarding all the twelve cases.As expected, these small deviations are mostly within the predicted uncertainty of our new benchmark relative energies.
Deviations from a few tenth to even 3 kcal mol À1 can be observed between previously published relative energies and ours, with small differences found mostly when authors carried out some kind of extrapolation of low-level energies.Note that, as these barrier heights (for X = Br and I provided for the reverse reactions in the literature) have small absolute values, the above deviations can mean very large relative errors with respect to the accurate energies.
The three possible mechanisms of the title reactions are also studied in terms of structural changes; the corresponding data are shown Fig. 8. Together with energetics, the structural parameters of the stationary points allow for making predictions concerning the dynamics of the title reactions and the validity of the Polanyi rules, 60 which were established for the dynamics of the simple A + BC reactions.As seen in Fig. 8, the UCCSD(T)-F12b geometries obtained with the aug-cc-pVDZ and with the aug-cc-pVTZ basis sets are basically the same usually within a few 0.001 Å in bond lengths and a few 0.1 degrees in bond angles, confirming the very fast basis set convergence observed for the relative energies (Table 1 and Fig. 3).For the H-substitution channel in all three TS geometries both the breaking C-H and the forming C-X bond lengths increase as going from F to I, with the largest increments between X = F and Cl.These C-H and C-X bonds are usually the shortest in the Walden-inversion TS, as it is the sterically most favorable mechanism due to back-side attack.The breaking C-H bond lengths in the FS as well as in the W TSs show increasing differences (26 to 35% for FS TSs and 20 to 35% for W TS) from the C-H bond in ethane, accompanied by the increasing resemblance of the C-X bond lengths of the TSs (19 to 13% for FS TSs and 17 to 8% for W TS) to those observed in the corresponding ethyl-halides as the atomic number of the halogen increases, thus one can talk about more and more product-like TSs.
In the case of methyl-substitution both the C-C and the C-X distances increase in the TS geometries from X = F to I with the C-C bond lengths being slightly shorter and the C-X bond lengths changing in a narrower scale in W TS. For the MS FS TS structures also the jump in the bond-length values between the X = F (1.882 Å for C-C and 1.733 Å for C-F) and Cl (2.080 Å for C-C 2.183 Å for C-Cl) cases is the largest.The X-C-C bond angles rise slightly from X = F to I as the size of the halogen atom increases.The C-C bond lengths in the TS structures increase (with 18-32% for both FS and W TSs) compared to ethane from X = F to I, while the C-X bond lengths are more and more similar, with differences decreasing by 20-16% for FS TS and by 22-11% for W TS to those observed in the methylhalides, indicating more and more product-like TSs.
As to H-abstraction, in the TS geometry the breaking C-H bond lengths increase from F to I, whereas the X-H distances decrease significantly from the extremely large value for X = F (1.941 Å). (Note that, for X = I the HA TS and postmin geometries are obtained at the less accurate CCSD/cc-pVDZ level of theory.)The X-H-C bond angle is also unusually bent for X = F (140.71), whereas for the other halogens it is around 1751 indicating a nearly collinear TS.In the postmin structure the H-X bonds are shorter (in the case of X = F by 50%) and the breaking C-H bonds are much longer compared to those of the TS geometries.H-X bond lengths in the HA TSs differ less and less from those of the HX products, by 53(F), 14(Cl), 4(Br), and 3(I)%, whereas in the postmin geometries these differences are only between 0-2% for all four cases.The breaking C-H bond in the TS geometries is getting more and more elongated from X = F to I (1-42%), while a 45-55% increase of the breaking C-H bond lengths can be observed in the postmin structures compared to ethane.Accordingly, for the HA channel the TS and postmin structures are more and more product-like, with a very reactant-like TS in the case of X = F, similar to the X + CH 4 reactions. 1ur findings regarding the geometrical changes during the title reactions are in accord with their energetics (Fig. 1) and thus with Hammond's postulate 61 stating that in exothermic-(endothermic) reactions the TS has a geometry similar to the reactants(products).Knowing the position of the TSs along the reaction coordinate also allows for making some wary predictions regarding the dynamical behavior of the title reactions based on the Polanyi rules. 60According to these basic rules, one may expect for the early-barrier F + C 2 H 6 reaction paths that translational energy will be the most efficient in promoting the reaction, while e.g. the late-barrier I + C 2 H 6 reactions will require strong vibrational excitation in order to increase reactivity.However, only dynamics simulations on accurate PESs can reveal the real effect of the complexity of these 9-atomic systems and the role of e.g. the entrance-channel van der Waals minima in these reactions.

Summary and conclusions
In this study we have thoroughly mapped the complex potential energy surfaces of the X + C 2 H 6 [F, Cl, Br, I] reactions using the explicitly-correlated UCCSD(T)-F12b method with different basis sets, concerning three possible reaction mechanisms: the already studied main H-abstraction, which is the only available path at low collision energies, and the so far uninvestigated higher-energy H-substitution and methyl-substitution reaction paths.For the substitution channels both the Walden-inversion and the front-side-attack mechanisms have been under consideration.We have identified the corresponding stationary points and report their benchmark-quality classical and adiabatic relative energies, including core-valence, post-CCSD(T) and spin-orbit corrections to the highly accurate UCCSD(T)-F12b/ aug-cc-pVQZ energies.It turns out, that methyl-substitution is the next available channel above the main H-abstraction pathway for all the title reactions, and at somewhat higher energies H-substitution also opens.Therefore, at high collision energies all the three reactions can occur, which could be experimentally investigated as was done in the high-energy crossed-beam study of the O( 3 P) + CH 4 multi-channel reaction. 62For both of the substitution channels the Walden-inversion pathway features lower barrier heights then front-side attack.Transition states follow the same order in energy for all four halogens (HA TS, MS W TS, HS W TS, HS FS TS, MS FS TS), except the X = I case, where the barrier of the front-side-attack pathway of H-substitution is higher than that of methyl-substitution.Product channels are also in the same increasing energy order in all four cases: H-abstraction, methyl-substitution and H-substitution.Our new benchmark 0 K reaction enthalpies are in excellent agreement with the best available experimental results.We have compared the relative energies obtained at the MP2/aug-cc-pVDZ level of theory, which is widely used for similar-sized systems, to the benchmark relative energies and we find MP2 results being far from chemical accuracy, with errors of several kcal mol À1 .Our results show that if subchemical accuracy, i.e. o0.5 kcal mol À1 , is desired for such chemical systems, one needs to take all of the above correction terms into account.Spin-orbit couplings have also been studied throughout the reaction paths revealing a nonnegligible effect even in the transition-state structures.We have also identified several minima in the entrance and exit channels along the different reaction paths, which may have a significant effect on the dynamics of the title reactions.We further investigated the structural changes during the course of the reactions and find fully reactant-like transition states for the X = F case, whereas more and more product-like TS geometries along with increasing endothermicity can be observed as the atomic number of the halogen increases, in accord with Hammond's postulate.Based on the basic dynamical rules proposed by Polanyi, the location of the TS along the reaction coordinate can give us hints concerning the effect of translational or vibrational excitation on reactivity, however, only accurate dynamics simulations can reveal the validity of these rules in systems of such complexity.In line with that, this study is a first step towards developing global analytical ab initio PESs, which provide opportunity to investigate the dynamics of the different channels of the title reactions in detail.Such a PES for the Cl + C 2 H 6 reaction is already being developed in our group.

Fig. 2
Fig. 2 Harmonic zero-point-energy corrections, determined at the UCCSD(T)-F12b/aug-cc-pVDZ level of theory, to the relative energies corresponding to the different stationary points of the X + C 2 H 6 [X = F, Cl, Br, I] reactions.In the case of the HA TS of the F + C 2 H 6 reaction the MP2/aug-cc-pVDZ ZPE correction, and for the HA TS and postmin of the I + C 2 H 6 reaction the UCCSD/cc-pVDZ ZPE corrections are shown.

Fig. 3
Fig. 3 Convergence of the UCCSD(T)-F12b relative energies with the aug-cc-pVDZ (DZ), aug-cc-pVTZ (TZ), and aug-cc-pVQZ (QZ) basis sets for the different stationary points of the X + C 2 H 6 [X = F, Cl, Br, I] reactions.In the case of the HA TS of the F + C 2 H 6 reaction neither ROHF, nor UHF converged at the UCCSD(T)-F12b/aug-cc-pVTZ geometry.
d[CCSDT] corrections are larger in magnitude than the d[CCSDT(Q)] terms, with average values of 0.3 and 0.2 kcal mol À1 , respectively.

Fig. 6
Fig. 6 Potential energy curves obtained at the MRCI+Q(5,3)/aug-cc-pVDZ(-PP) (PPs for Br and I) level as a function of the C 2 H 6 Á Á ÁX [X = F (light green), Cl (dark green), Br (red), I (purple)] separation, considering three different separation directions: (left panel) X approaching ethane from one of the methylgroups maintaining a C 3v point-group symmetry, (middle panel) X approaching one H atom of one of the methyl groups, and (right panel) X approaching perpendicularly the C-C bond of ethane, see inserted geometries.The C 2 H 6 unit is kept frozen at its CCSD(T)-F12b/aug-cc-pVTZ equilibrium structure while changing the separation distance of X. Solid symbols represent non-spin-orbit ground states (non-SO 1 ), while open symbols refer to SO ground states (SO 1 ).SO ground states are shifted with their asymptotic (10 Å separation) limits in each case.

Fig. 7
Fig. 7 Deviation in percent between the computed (MRCI+Q(5,3)/aug-cc-pVDZ(-PP), PPs for Br and I) energy difference of the spin-orbit and nonspin-orbit ground states regarding each stationary points of the X + C 2 H 6 [X = F (light green), Cl (dark green), Br (red), I (purple)] reactions and the 1/3 of the experimental SO splitting (energy difference between the 2 P 3/2 and the 2 P 1/2 states) of the halogen atoms.

Fig. 8
Fig. 8 Relevant bond lengths (in Å) and bond angles (in degrees) of the stationary-point structures of the X + C 2 H 6 [X = F (light green), Cl (dark green), Br (red), I (purple)] reactions obtained at the UCCSD(T)-F12b/aug-cc-pVTZ (UCCSD(T)-F12b/aug-cc-pVDZ) levels of theory.In the case of the HA TS of the F + C 2 H 6 reaction the UCCSD(T)-F12b/aug-cc-pVDZ geometry could not be converged (superscript a).The HA TS and postmin geometries of the I + C 2 H 6 reaction, and for comparison in this figure the HI geometries are determined at the UCCSD/cc-pVDZ level of theory (superscript b).

Table 1
Benchmark classical and adiabatic energies (kcal mol À1 ) of the stationary points obtained from eqn (1) and (5) relative to the reactants for the X + C 2 H 6 [X = F, Cl, Br, I] reactions Stationary points MP2/DZ a MP2 D ZPE

Table 2
Comparison between the best available experimental and our computed benchmark 0 K reaction enthalpies, given in kcal mol À1 , for the X + C 2 H 6 [X = F, Cl, Br, I] reactions -0.4 kcal mol À1 uncertainty.The harmonic ZPE corrections are supposed to account for approximately the 95% of the anharmonic correction resulting in about a 0.2 kcal mol À1 This journal is © the Owner Societies 2019 Phys.Chem.Chem.Phys., 2019, 21, 396--408 | 407