Cangtao
Yin
*,
Viktor
Tajti
and
Gábor
Czakó
*
MTA-SZTE Lendület Computational Reaction Dynamics Research Group, Interdisciplinary Excellence Centre and Department of Physical Chemistry and Materials Science, Institute of Chemistry, University of Szeged, Rerrich Béla tér 1, Szeged H-6720, Hungary. E-mail: yincangtao@dicp.ac.cn; gczako@chem.u-szeged.hu
First published on 21st September 2022
We report a full-dimensional spin–orbit-corrected analytical potential energy surface (PES) for the HBr + C2H5 → Br + C2H6 reaction and a quasi-classical dynamics study on the new PES. For the PES development, the ROBOSURFER program package is applied and the ManyHF-based UCCSD(T)-F12a/cc-pVDZ-F12(-PP) energy points are fitted using the permutationally-invariant monomial symmetrization approach. The spin–orbit coupling at the level of MRCI-F12+Q(5,3)/cc-pVDZ-F12(-PP) is taken into account, since it has a significant effect in the exit channel of this reaction. Our simulations show that in the 1–40 kcal mol−1 collision energy (Ecoll) range the b = 0 reaction probability increases first and then decreases with increasing Ecoll, reaching around 15% at the medium Ecoll. No significant Ecoll dependence is observed in the range of 5–20 kcal mol−1. The reaction probabilities decrease monotonically with increasing b and the maximum b where reactivity vanishes is smaller and smaller as Ecoll increases. Unlike in the case of HBr + CH3, the integral cross-section decays sharply as Ecoll changes from 5 to 1 kcal mol−1. Scattering angle distributions usually show forward scattering preference, indicating the dominance of the direct stripping mechanism. The reaction clearly favors H-side attack over side-on HBr and the least-preferred Br-side approach, and favors side-on CH3CH2 attack over the CH2-side and the least-preferred CH3-side approach. The initial translational energy turns out to convert mostly into product recoil, whereas the reaction energy excites the C2H6 vibration. The vibrational and rotational distributions of the C2H6 product slightly blue-shift as Ecoll increases, and very few reactive trajectories violate zero-point energy.
The dynamics of two other similar chemical reactions have been theoretically investigated in our group: the study on the Cl + C2H6 reaction provided unprecedented agreement with the experiment for the rotational state distribution of the HCl product.7 In the meanwhile the vibrationally-resolved rotational state distributions of the HF product obtained from the computations agree well with the single-collision experimental data for the v = 1, 2, and 3 states, in the case of the F + C2H6 reaction.6 The HBr + C2H5 reaction has been the subject of both theoretical and experimental investigations due to their supposedly non-Arrhenius behavior.12–17 Experimentally, Seetula12 investigated the kinetics of the title reaction under pseudo-first-order conditions in a heatable tubular reactor and fitted the pressure-independent rate constants determined by the Arrhenius expression: k = (1.87 ± 0.14) × 10−12 exp[+(3.7 ± 0.2) kJ mol−1/RT]. They also reported the formation values of entropy (244 ± 6 J K−1 mol−1) and enthalpy (120.7 ± 2.1 kJ mol−1) for the radical studied at 298 K using the above kinetic data in a second-law procedure. Leplat et al.13 performed a reinvestigation of the absolute rate constants and led to the following Arrhenius expressions: k = 3.69(±0.95) × 10−11 exp(−10.62(±0.66) kJ mol−1/RT) in the temperature range of 293–623 K using a Knudsen reactor coupled to a single-photon photoionization mass spectrometer. They recommended the standard heat of formation of C2H5 to be 117.3 ± 3.1 kJ mol−1 resulting from an average of “third law” evaluations using S298°(C2H5) = 242.9 ± 4.6 J K−1 mol−1. Theoretical investigations began with the evaluation of absolute rate constants and kinetic isotope effects by transition state theory and RRKM theory as applied to the dissociation of the intermediate complex at the (U)MP4/6-311G**//(U)MP2/6-31G* level.14 Seetula15 calculated the enthalpy of formation of C2H5 at 298 K to be 120.4 ± 2.7 kJ mol−1 by optimizing the transition state of the reaction HBr + C2H5 → Br + C2H6 at the MP2(fc)/6-31G(d,p) level of theory, and studied the chemical nature of the transition state by localizing it along the minimum energy path of the reaction. Sheng et al.16 studied the reaction over the temperature range from 200 to 1400 K. They calculated the electronic structure information at the BHLYP/6-311+G(d,p) and QCISD/6-31+G(d) levels, and further refined the energies along the minimum energy paths at both levels by performing the single-point calculations at the PMP4(SDTQ)/6-311+G(3df,2p)//BHLYP and QCISD(T)/6-311++G(2df,2pd)//QCISD levels. Their computed results further indicated that the title reaction has a negative temperature dependence at T < 850 K, but clearly shows a positive temperature dependence at T > 850 K, and predicted that the kinetic isotope effect for the title reaction is inverse in the temperature range from 200 to 482 K. Golden et al.17 employed RRKM theory to analyze the kinetics of the title reaction. They characterized the stationary points along the reaction coordinate with coupled cluster theory combined with basis set extrapolation to the complete basis set limit, and located a shallow minimum bound by 9.7 kJ mol−1 relative to the reactants, with a very small energy barrier to dissociation to the products. They reported that the transition state is tight compared to the adduct, explored the influence of vibrational anharmonicity on the kinetics and thermochemistry of the title reaction, and suggested that with adjustment of the adduct binding energy by around 4 kJ mol−1, the computed rate constants may be brought into agreement with most experimental data in the literature. They indicated that at temperatures above those studied experimentally, the activation energy may switch from negative to positive. As an initial step for the dynamics study reported here, in 2019 our group18 determined benchmark geometries and energies for the stationary points of the backward reaction Br + C2H6, also considering the H-substitution and the methyl-substitution reaction pathways, 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. 0.5 kcal mol−1, accuracy.
In 2013 one of us worked out the potential energy function of a similar but smaller reactive system, Br + CH4.19 Later, using the quasi-classical trajectory (QCT) method in combination with an improved version of the high-level ab initio PES, Góger et al.20 also provided reliable rate coefficients for the HBr + CH3 → Br + CH4 reaction at combustion temperatures where no experimental data are available. The QCT calculations had been validated by reproducing the experimental rate coefficients at room temperature. At temperatures between 600 and 3200 K, the QCT rate coefficients display positive activation energies. Very recently, Gao et al.21 performed quantum and quasi-classical dynamics calculations for this reaction and very good agreement was found between the two approaches. The cross-sections were found to diverge when Ecoll decreases, indicating that the reactant attraction is responsible for the dynamics at low Ecoll. The quantum mechanical and the quasi-classical rate constants also agree very well and almost exactly reproduce the experimental results at low temperatures up to 540 K. The negative activation energy observed experimentally is confirmed by the calculations and is a consequence of the long-range attraction between the reactants. They found that at very low Ecoll, the reactive system crosses the potential barrier because the forces within the complex guide them. When Ecoll increases, the system does not follow the most favorable path and the reactants are, with increasing probability, reflected from the repulsive walls of the nonreactive parts of the reactants, providing a picture beyond the decreasing excitation function.
In this study, we aim to give a better understanding of how the HBr + C2H5 → Br + C2H6 reaction proceeds step-by-step at the atomic and molecular level. A full-dimensional ab initio PES for the title reaction is developed to provide a correct description. The spin–orbit coupling is taken explicitly into account in our computations, since it has a significant effect in the exit channel of this reaction, especially when the products are far from each other. Using this PES, we carry out QCT simulations, and discuss detailed dynamics results of the title reaction.
For the fitting of the energy points of the PES we utilize the Monomial Symmetrization Approach (MSA).29 Within this approach the PES is fitted using a full-dimensional analytical function which is inherently invariant under the permutations of like atoms. This function is an expansion of polynomials of the yij = exp(−rij/a) Morse-like variables, where rij are the inter-atomic distances and the a parameter, set to 2.0 bohr based on careful tests, controls the asymptotic behavior of the PES. The highest order of the polynomials is 5. The energy points are fitted using a least-squares fit with an E0/(E + E0) weighing 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 ROBOSURFER program a hard energy limit of 150 kcal mol−1 relative to the energy of the free reactants is applied, above which no energy point is added. A hard energy limit of 30 kcal mol−1 below the reactants is set to avoid spurious minima. A 0.5 kcal mol−1 target accuracy of the fitting is set. These parameter values are chosen to obtain the best possible description of the HBr + C2H5 reaction. The ManyHF-based RMP2/aug-cc-pVDZ (aug-cc-pVDZ-PP for Br atom) level of theory is used for the initial PES development, which consists of 128 ROBOSURFER iterations in total. During the PES development QCT computations30 are run to obtain new geometries, where the Ecoll is set from 1 to 60 kcal mol−1, which is enough to cover the energies in the PES we are interested in. We use both sides of the reaction, i.e., HBr + C2H5 and Br + C2H6, as the starting points of the QCTs. At the end we have 12819 geometries and the energies of them are recomputed at the following composite level of theory: ManyHF-UCCSD(T)-F12a/cc-pVDZ-F12 + SOcorr(MRCI-F12+Q(5,3)/cc-pVDZ-F12) (cc-pVDZ-PP-F12 for Br atom), where “SOcorr” is the spin–orbit (SO) correction to each energy point. In free halogen atoms the non-relativistic 2P ground electronic state splits due to the relativistic spin–orbit interaction, and the relativistic 2P3/2 ground state (where the total angular momentum quantum number, J = 3/2) is lower by one third of the splitting energy with respect to the non-relativistic ground state. Thus, this energy-lowering effect is especially relevant in the exit (Br + C2H6) channel of the reactive PES, where the Br atom is far from and thus unbound to the ethane molecule. The SO correction is determined at the MRCI-F12+Q(5,3)/cc-pVDZ-F12(-PP)31 level of theory for each geometry. The multireference computations utilize a minimal active space of 5 electrons on 3 spatial 4p-like orbitals, and the Q Davidson-correction32 estimates higher-order correlation energy effects. The SO computations make use of the Breit–Pauli operator in the interacting-states approach,33 where the SO eigenstates are determined by diagonalizing the 6 × 6 SO matrix whose diagonal elements are replaced by the Davidson-corrected MRCI energies.
The MSA is used again for the fitting of the geometries with new energies of the PES. The final PES is built from 11364 geometries and the corresponding composite energies, whose distribution is shown in Fig. 1, and features small root mean square (RMS) errors in the chemically interesting energy ranges as shown in Table 1.
![]() | ||
Fig. 1 Distribution of the final sampling points. Energies are given with respect to the asymptote of reactants. |
Energy range (kcal mol−1) | 0–20 | 20–40 | 40–100 |
---|---|---|---|
RMS (kcal mol−1) | 0.248 | 0.487 | 1.015 |
Integral cross-sections (σ) for the title reaction are calculated by a b-weighted numerical integration of the P(b) opacity functions at each Ecoll. The product ZPE-constraint, i.e., the classical vibrational energy of C2H6 has to be greater than its ZPE on the present PES, only rules out 8 trajectories from the 5793 reactive trajectories. So from now on we can safely ignore the differences between features with and without ZPE constraint.
The scattering angle distributions are obtained by binning the cosine of the angle (θ) of the relative velocity vectors of the center of masses of the products and those of the reactants into 10 equidistant bins from −1 to 1. cos(θ) = 1 (θ = 0°) corresponds to forward scattering and cos(θ) = −1 (θ = 180°) corresponds to backward scattering. The initial attack angle distributions for the reactants are calculated by binning the cosine of the angle (α for HBr and β for C2H5) of the velocity vector of center of mass of the examined reactant and an interatomic vector that is considered as the Br–H bond for HBr and the C–C bond for C2H5. We also use 10 equidistant bins between −1 to 1 like in the case of scattering angle distributions. For HBr cos(α) = −1 means that HBr approaches with its Br atom side and in the situation of cos(α) = 1 HBr goes with its H atom towards C2H5. While for C2H5 cos(β) = −1 means that C2H5 approaches HBr with its CH3 side and in the situation of cos(β) = 1 C2H5 goes with its CH2 side towards the HBr. Rotational quantum numbers of C2H6 are obtained by rounding the lengths of classical rotational angular momentum vectors to the nearest integer values.
![]() | ||
Fig. 2 Schematic potential energy diagram of the HBr + C2H5 → Br + C2H6 reaction comparing the classical relative energies obtained on the present PES, ManyHF-UCCSD(T)-F12a/cc-pVDZ-F12 + SOcorr(MRCI-F12+Q(5,3)/cc-pVDZ-F12) (cc-pVDZ-PP-F12 for Br atom) energies at geometries optimized on the PES, and the relativistic all-electron CCSDT(Q)/complete-basis-set-quality benchmark relative energies18 of the stationary points, in kcal mol−1. |
The imaginary frequency of the TS is only 162 cm−1 so it is a soft TS, which means the energy is a little higher than nearby. The lowest normal mode of premin has a wavenumber of 61 cm−1 so it is a shallow minimum, its energy is a little lower than nearby. It is possible that the TS lies below the premin, as shown at the benchmark level and on the PES. Note that at the ab initio level used for the PES development this trend is reversed, but with only 0.1 kcal mol−1 difference. Actually in this reaction the TS and the premin have very similar energies and their relative energy is within our error bar.
Besides the comparison of the stationary-point energies, the one-dimensional potential energy curve as a function of the C–Br distance obtained on the PES, compared to that obtained using the ab initio method, is another good tool to evaluate the accuracy of the PES. The PES describes the entrance channel potential accurately as shown in Fig. 3.
To capture the properties near TS, we took the TS geometry, elongated and shortened the C–Br distance, while freezing the other degrees of the two parts, CH3CH2 and HBr. The asymptote in Fig. 3 is not the reactants in equilibrium, but bent over a little, so its energy is higher than that of the reactants, leading to a deeper well with −5.1 and −4.8 kcal mol−1 relative energies with respect to the asymptote in the case of the PES and the direct ab initio method, respectively. The comparison of the two minima indicates low (0.3 kcal mol−1) fitting errors of the full-dimensional PES, again consistent with the small RMS values of Table 1. The positions of the two minima located at C–Br distances of 3.10 and 3.22 Å, only differ by 0.12 Å, showing good behavior of the analytical PES.
![]() | ||
Fig. 4 Reaction probabilities as a function of the b impact parameter for the title reaction at different collision energies (given in kcal mol−1). |
As seen in Fig. 4, no threshold energy above 1.0 kcal mol−1 can be observed for the title reaction to proceed, in accordance with the negative barrier height relative to the reactants. It is clear that the opacity function at 1 kcal mol−1Ecoll has a much smaller reaction probability and much further distance where the reaction can still occur. The b = 0 reaction probability increases first and then decreases with increasing Ecoll, reaching around 15% in the medium Ecoll. No significant Ecoll dependence is observed in the range of 5–20 kcal mol−1. The reaction probabilities decrease monotonically with increasing b and the maximum b where reactivity vanishes decreases as Ecoll increases. The bmax value is the largest (14 bohr) at 1 kcal mol−1Ecoll, due presumably to the fact that the dipole–dipole interaction between the reactants is the least counteracted by translational momenta. Considering ZPE (41 kcal mol−1 for reactants and 43 kcal mol−1 for TS), the sub-barrier energy is −1.7 kcal mol−1. In shallow sub-barrier reactions, the reaction probability could be small at low Ecoll.6,10
The integral cross-section reaches its highest value at around 5 to 10 kcal mol−1Ecoll. At the high Ecoll part it decreases slowly, consistent with the HBr + CH3 → Br + CH4 reaction.19 However when Ecoll is low it decreases sharply with decreasing Ecoll, due probably to the fact that the reactants’ interaction is more anisotropic in the case of C2H5.
It is clear that forward scattering is favored at high Ecoll (40 and 20 kcal mol−1) indicating the dominance of the direct stripping mechanism. As Ecoll decreases (10 and 5 kcal mol−1), the distributions are less anisotropic. The DCS is almost isotropic for 2 kcal mol−1. However, it is quite surprising that the DCS at Ecoll = 1 kcal mol−1 shows a forward peak, which may indicate a different mechanism at low Ecoll. Nevertheless, this forward scattering is consistent with the significantly larger bmax value at low Ecoll (see Fig. 4), because at large b usually stripping occurs resulting in forward scattering. The title reaction favors H-side attack over side-on HBr and the least-preferred Br-side approach, as expected, because an H–C bond forms in the abstraction process, and favors side-on CH3CH2 attack over CH2-side and the least-preferred CH3-side approach. None of these attack characters shows any significant Ecoll dependence, except that the HBr attack angle distribution is more isotropic at the lowest Ecoll.
In order to provide more insight into the scattering dynamics, we present the correlation diagrams of the impact parameters and scattering angles in Fig. 7. At higher collision energies the reaction with direct forward scattering is favored for large b while backward scatting is favored for small b as expected. At low Ecoll some correlation between b and cos(θ) still exists showing that the isotropic scattering angle distribution at Ecoll = 2 kcal mol−1 rather originates from the balance of direct rebound and stripping than from indirect dynamics. At Ecoll = 1 kcal mol−1 the b − cos(θ) correlation is less significant showing some indirect feature. Furthermore, the forward-scattered trajectories at b > 9.5 bohr, which are not present at Ecoll = 2 kcal mol−1, boost the forward-scattered peak of the scattering angle distribution at Ecoll = 1 kcal mol−1 due to the b-weighted integration over the impact parameters.
![]() | ||
Fig. 7 Correlation diagrams of the impact parameters and scattering angles for the title reaction at different collision energies (given in kcal mol−1). |
![]() | ||
Fig. 8 Normalized product relative translational energy distributions for the title reaction at different collision energies (given in kcal mol−1). |
As shown in Fig. 8, the distributions become broader as the Ecoll increases, and their maxima are shifted by almost the total increment of the Ecoll, indicating that the major part of the initial translational energy ends up in translational recoil in all cases.
Consistent with the above observation, the Ecoll dependence of the internal energy distribution of ethane, plotted in Fig. 9, also suggests that only a small portion of the collision energy is transferred into the vibrational and rotational degrees of freedom of ethane, since the peaks of the distributions are less affected by the change in the initial translational energy than those of the product relative translational energy distributions. The internal energy excitations of the product ethane mainly come from the reaction energy.
As shown in Fig. 9, considering that the ZPE of ethane is 47 kcal mol−1, very few trajectories are ZPE-violating, consistent with our previous points arguing that ZPE-constraint only rules out 8 reactive trajectories.
This journal is © the Owner Societies 2022 |