Reaction Kinetics of HBr with HO2: A New Channel for Isotope

Oct 10, 2016 - Reaction Kinetics of HBr with HO2: A New Channel for Isotope .... Mergers and acquisitions in the drug industry last year brought one ...
0 downloads 0 Views 1004KB Size
Download PDF
Article pubs.acs.org/JPCA

Reaction Kinetics of HBr with HO2: A New Channel for Isotope Scrambling Reactions Jonathan R. Church and Rex T. Skodje* Department of Chemistry and Biochemistry, University of Colorado, Boulder, Colorado 80309-0215, United States S Supporting Information *

ABSTRACT: The gas phase reaction kinetics of HBr with the HO2 radical are investigated over the temperature range of T = 200−1500 K using a theoretical approach based on transition state theory. The parameters for the potential energy surface are computed using density functional theory with the M11 exchange functional. The rate coefficient for the HBr + HO2 → Br + H2O2 abstraction channel is found to be somewhat larger than previous estimates at low temperatures due to quantum tunneling. The present study reveals the existence of a novel exchange pathway, HBr + H′O2 → H′Br + HO2, which exhibits a much lower reaction barrier than does the abstraction route. The transition state for this process is a symmetrical planar five-membered-ring-shaped structure. At low temperatures, this concerted double hydrogen transfer reaction is several orders of magnitude faster than the abstraction channel. The exchange process may be observed using isotope scrambling reactions; such reactions may contribute to observed isotope abundances in the atmosphere. The rate coefficients for the isotopically labeled reactions are computed.

I. INTRODUCTION The gas phase chemical reactions of small bromine-containing species continues to attract interest in atmospheric1 and combustion applications. In combustion chemistry,2,3 it is known that the introduction of halogen species can act as an inhibitor such as used in flame retardants.4,5 Mathieu et al.,6 e.g., have recently discussed the influence of radical reactions such as C2H3 + Br = C2H3Br and Br + HO2 = HBr + O2 on the laminar flame speed and ignition delay time for hydrocarbon combustion. The sensitivity of combustion observables to the rates of elementary bromine processes motivates the need to continue to improve the accuracy of rate coefficients used in chemical modeling. In atmospheric chemistry, bromine reactions such as HBr + O3 = BrO + HO2, Br + HO2 = HBr + O2, and HBr + HO2 = H2O2 + Br were studied in the past in connection with stratospheric ozone loss.7−11 More recently, bromine released from sea-salt aerosols in the troposphere is suspected to be of importance for reactions of secondary organic aerosols and in ozone consumption.12,13 Again, the accurate determination of elementary rate coefficients is important for the construction of reliable models of such atmospheric chemistry. In this work we consider the reactions of HBr with HO2. The HO2 radical is an essential species for both combustion and atmospheric chemistry, whereas HBr is one of the most important carriers of the bromine species. The most wellknown reaction channel is the endothermic abstraction process

the absence of a clear signal of the reaction, set an upper limit for k1 to be 3.0 × 10−17 (cm3 molecule)/s at 300 K. Leu,8 Posey et al.,9 and Toohey et al.14 separately studied the exothermic reverse reaction Br + H 2O2 → HBr + HO2

and set upper bounds for k1r in the range of (0.5−3.0) × 10−15 cm3/(s·molecule) at ∼300 K. Heneghan and Benson15 also studied reaction R1r in the range 300−350 K and found a somewhat higher rate, k1r(T=300 K) = (1.3 ± 0.45) × 10−14 cm3/(s·molecule). Using the equilibrium constant, the value of k1 inferred from k1r was also found to be less than 10−15 cm3/(s· molecule) at ambient temperatures. In an early experimental study of bromine inhibition of the second explosion limit of hydrogen, Clark et al.16 estimated that k1r = (1.7 ± 0.8) × 10−13 at the higher temperature of T = 773 K and hence k1 = 1.6 × 10−14 cm3/(s·molecule). Experimental results at higher combustion temperatures are not available at present. On the basis of the activation energy, it is estimated that the reaction barrier for reaction R1 is on the order of 10 kcal/mol. To provide a rate coefficient for reaction R1 appropriate for combustion temperatures, and to resolve some of the discrepancies reported at lower temperatures, Dixon-Lewis et al.17 have recently carried out a theoretical evaluation of k1(T) based on transition state theory (TST), which uses a quantum chemical determination of the transition state properties. The derived expression predicts a value of k1 = 3.2 × 10−20 cm3/ (molecule·s) at T = 300 K which is 3 orders of magnitude below the upper bound set by experiment. Hence, it is possible that barrier to reaction is significantly lower than calculated,

HBr + HO2 → H 2O2 + Br ΔHrxn(T =0) = +2.8 kcal/mol

(R1)

Experimental work near room temperature has concluded that the rate of this process is quite slow. The only attempted direct measurement of reaction R1 was by Mellouki et al.10 who, in © XXXX American Chemical Society

(R1r)

Received: July 18, 2016 Revised: October 1, 2016

A

DOI: 10.1021/acs.jpca.6b07215 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 1. Mechanism for reaction R2. The double hydrogen exchange reaction occurs through the planar transition state TS2 depicted here.

the first to use ab initio quantum chemistry to explore the degree of stabilization in the hydrogen-bonded ring structures for SO3 + (H2O)n = H2SO4 + (H2O)n−1 reactions. Since then, there have been a number of studies of water catalyzed and carboxylic acid catalyzed gas phase reactions that have emphasized the stabilization of ring-shaped transition states.30−38 The degree of ring stabilization can be understood, to some degree, by appeal to concepts from molecular structure theory such as ring strain and electron delocalization. Of course, there are other factors that influence the rate of a reaction besides the barrier height. The formation of a cyclic transition state is expected to have a lower entropy of activation because the ring is more constrained than, e.g., a linear transition state. The degree of quantum tunneling may also be strongly affected by width of the reaction barrier. In our previous work on water catalyzed reactions, we have also shown that dynamical effects can be crucial in light activated reactions39−41 and may actually cancel the barrier lowering effect. In this work we shall compute the rate coefficients for reactions R1 and R2, and their isotopomers, using an approach based on ab initio transition state theory. In section II we present the theoretical methods used in this study. The barrier heights, structures, and vibrational frequencies are computed using density functional theory. The performances of several exchange functionals were compared. The rate coefficients were computed using TST with the harmonic oscillator/rigid rotor approximation for the partition function. It was found that variational modifications to the values were unimportant in the present case. The effects of quantum tunneling were incorporated into the rate coefficients using the small curvature tunneling, zero curvature tunneling, and Wigner methods. The results of the calculations are presented in section III. The rate coefficients were computed over a temperature range of 200− 1500 K. It is found that reaction R2 proceeds much faster at low temperatures than does reaction R1. At high temperatures, however, the entropic effects favor reaction R1 and thus the rate coefficients for reactions R1 and R2 become comparable. The role of quantum tunneling is discussed for both reactions. The role of isotopic substitution is explored for both reactions. In section IV, the significance of the results is discussed.

ΔE0 = 10.5 kcal/mol, or that quantum tunneling plays a larger role than predicted. One goal of the present work is to investigate these possibilities. In addition to reaction R1, it has been noted that there are other potential reaction channels for the HBr + HO2 system.19 Indeed, the reactions HBr + HO 2 → HOBr + OH (ΔH(T=298K) = −5.6 kcal/mol) and HBr + HO2 → BrO + H2O (ΔH(T=298K) = −28.3 kcal/mol) are both more exothermic than reaction R1. However, it does not appear that either of these channels has been actually observed in experiment. Furthermore, our theoretical calculations indicate that the saddlepoints leading to these channels lie much higher in energy than that for reaction R1. However, in the course of this work we have discovered another much lower barrier pathway for the HBr + HO2 system. This channel involves the simultaneous exchange of two H atoms, HBr + H′O2 → H′Br + HO2

(R2)

Naturally, this channel is unobservable unless the system is isotopically labeled, e.g. HBr + DO2 → DBr + HO2

(R2a)

DBr + HO2 → HBr + DO2

(R2b)

HBr + H16O18O → HBr + H18O16O

(R2c)

Such isotope scrambling reactions may in fact prove to be important in understanding isotope abundances in atmospheric trace gases.18,19 Our interest here, however, is in understanding the basic chemistry behind the double exchange process, reaction R2. The HO2 radical is known to form multiply hydrogenbonded ring adducts with other species.20−22 For example, in the disproportionation self-reaction HO2 + HO2 → H2O2 + O2 the radicals come together on the triplet surface to produce a metastable six-membered-ring-shaped intermediate that is bound by approximately 10 kcal/mol relative to HO2 + HO2.23−27 Toohley et al.14 noted that reactions such as Br + HO2 → HBr + O2 may proceed through a four-center intermediate and transition state. Computational studies have likewise demonstrated that HO2 readily forms ring-shaped intermediates and transition states in gas phase reactions for reactions such as HO2 + OH. We propose that reaction R2 proceeds through the five-membered transition state, as shown in Figure 1. In this view, the reaction is a concerted double hydrogen exchange process. We also shall find that HBr---HO2 forms a ring-shaped prereactive complex that may play a role in the reaction kinetics at very low temperatures. However, at room temperature and above, the complex potential well is not sufficiently deep to significantly influence the reaction. The lowering of reaction barriers for hydrogen exchange reactions due to the formation of cyclic hydrogen-bonded transition states is a well explored idea dating back, at least, to the work of Eigen.28 Morokuma and Muguruma29 were among

II. THEORY In this section, we briefly describe the theoretical methods used to obtain the rate coefficients for the HBr + HO2 reactions (reactions R1 and R2), and their isotopomers. The results of the quantum chemical calculations (QM) are presented in this section. This includes the optimized structures, vibrational frequencies, and energetics of the transition states and the reactive intermediate. We also review the procedure for computing the TST rate coefficients with tunneling corrections. IIA. Quantum Chemistry. The QM calculations are challenging due to the presence of the bromine atom, which requires the use of a large basis set and relativistic corrections B

DOI: 10.1021/acs.jpca.6b07215 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Table 1. Zero-point-Corrected Energy of the Stationary Points for the HBr + HO2 System Computed at Various Levels of Theory M11/cc-pVnZ ΔE0(kcal/mol)a reactants intermediate TS1 TS2 products ΔE0(kcal/mol)

n=D

n=T

0.00 0.00 −5.42 −3.86 7.34 9.12 −1.58 2.00 7.10 5.05 M08-HX/cc-pVnZ//M11

wB97X/cc-pVnZ//M11 n=Q

CBS

ΔE0(kcal/mol)

0.00 −3.36 9.80 2.95 4.23

0.00 −3.11 10.04 3.56 3.95

reactants intermediate TS1 TS2 products

n=Q

CBS

ΔE0 (kcal/mol)

0.00 0.00 −6.24 −4.70 7.50 8.71 −0.67 2.54 7.19 4.47 M06-2X/cc-pVnZ//M11

0.00 −4.04 9.50 2.95 4.61

0.00 −3.85 9.58 3.63 3.90

reactants intermediate TS1 TS2 products

n=D

n=Q

CBS

ΔE0(kcal/mol)

0.00 0.00 0.00 −5.62 −4.39 −3.84 6.02 7.06 7.78 −1.11 1.66 2.24 2.44 −0.28 −0.04 wB97X-D/cc-pVnZ//M11

0.00 −3.69 7.83 2.76 −0.78

reactants intermediate TS1 TS2 products

ΔE0(kcal/mol)

n=D

n=T

n=Q

CBS

ΔE0 (kcal/mol)

reactants intermediate TS1 TS2 products

0.00 −5.17 7.09 1.87 8.00

0.00 −3.77 7.90 3.93 5.70

0.00 −3.20 8.49 4.62 5.62

0.00 −3.02 8.52 4.93 5.08

reactants intermediate TS1 TS2 products ΔE0(kcal/mol) reactants intermediate TS1 TS2 products

n=D

n=T

n=T

reactants TS1 ΔE0 (kcal/mol) reactants products

n=D

n=Q

CBS

0.00 0.00 −6.20 −4.77 9.16 9.94 1.71 4.05 8.04 5.49 wB97/cc-pVnZ//M11

n=T

0.00 −4.14 10.61 4.88 5.46

0.00 −3.97 10.61 5.22 4.84

n=D

n=Q

CBS

0.00 0.00 −6.63 −5.04 10.88 11.78 2.37 5.09 8.22 5.67 B3LYP/cc-pVnZ//M11

n=T

0.00 −4.38 12.55 6.06 5.61

0.00 −4.18 12.56 6.45 5.00

n=D

n=Q

CBS

0.00 0.00 0.00 −5.42 −2.93 −2.62 2.83 4.59 4.80 0.05 4.03 4.50 8.78 6.25 5.65 CCSD(T)/cc-pVnZ//B3YLP

0.00 −2.10 5.18 5.35 5.20

n=T

n=T

n=Q

CBS

0.00 0.00 10.92 10.69 experimental

0.00 10.50 M11/CBS

0.00 2.80

0.00 3.95

a

The harmonic zero-point-corrected energies of stationary points for the reagents, products, transition states, and intermeadiate are computed with the zero set to the reagent energies.

vibrational frequencies that incorporate the scaling factors are provided for the stationary points. In Figure 2, we show the optimized geometries of the two transition states TS1 and TS2 along with that of the intermediate complex. The rotational constants for the transition states and the reagent species are reported in the Supporting Information. To put the quantum chemistry results into context, in Table 1 we also present the zero-point-corrected (adiabatic) energies obtained using several other methods. The M11/CBS adiabatic barrier energy for TS1, corresponding to HBr + HO2 → Br + H2O2, is 10.04 kcal/mol whereas the intermediate complex lies at −3.11 kcal/mol, below the reactant asymptote. The prediction of M08-HX/CBS//M11/cc-pVQZ-PP for the adiabatic barrier is 9.58 kcal/mol and the intermediate is at −3.85 kcal/mol, which are in good agreement with the M11 results. The barrier energy for M06-2X/CBS//M11/cc-pVQZPP is 7.83 kcal/mol, which is probably in greater error than M11. The dispersion-corrected ωB97X-D/CBS//M11/ccpVQZ-PP method gives a barrier of 8.52 kcal/mol whereas the form without dispersion ωB97X/CBS//M11/cc-pVQZ-PP gives 10.61 kcal/mol. The conventional B3LYP functional seems to be an outlier with a low barrier of 5.18 kcal/mol. It is also useful to compare the predictions for the reaction endothermicity of the various methods to the experimental result of ΔH(T=0) = +2.8 kcal/mol. It is seen that M11 (and M08-HX) does the best with an error of about +1 kcal/mol. The other methods have clearly larger errors for this quantity.

for the core electrons. To simplify the problem, the QM has been studied using density functional theory (DFT) with basis sets that employ an effective core potential to model the inner bromine electrons, viz. the cc-pVnZ-PP basis sets with n = D, T, or Q. A number of recent studies has verified that DFT can yield reliable geometries, frequencies, and energetics for small bromine-containing systems.42−44 The “Minnesota” functionals have proven to be among the most dependable. The M11 and M11-L functionals are the most recent and highly optimized of these expressions.45 It appears that the commonly used B3LYP functional fares poorly in a number of benchmark calculation for bromine-containing systems.42,43 We shall employ the M11 functional in this work. The geometries of the saddlepoints for reactions R1 and R2, which we designate as TS1 and TS2, were optimized using the M11/cc-pVQZ-PP method. The optimized structures of the reactant, product, and the intermediate species were likewise determined. The nuclear geometries are reported in the Supporting Information. The vibrational frequencies were computed at the same level of theory using a scaling factor of 0.967.46 The final energies for the stationary points, i.e., the saddlepoints and the reactant and product potential minima, were obtained through extrapolation to the complete basis set (CBS) limit based on M11/cc-pVnZ (n = D, T, Q). All QM calculations were carried out using the GAMESS-US software package.47 The results for the zero-point-corrected energies of the stationary points are presented in Table 1. In Table 2, the C

DOI: 10.1021/acs.jpca.6b07215 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 2. Frequencies (cm−1) for the Stationary Points of the HBr + HO2 System Computed Using M11/cc-pVQZ-PPa HBr DBr HO2 DO2 intermediate HO2) intermediate HO2) intermediate DO2) intermediate DO2) TS1

(HBr + (DBr + (HBr + (DBr +

TS1 (DBr + HO2) TS2 TS2 (HBr + DO2) TS2 (DBr + HO2) TS2 (DBr+ DO2) H2O2 DHO2

We have also considered the possibility that a multireference character of the wave function may affect the results. To test this, the T2-diagnostic48 for multireference character was computed at the TS1 and TS2 structures. Our largest values for the T2 amplitudes of TS1 and TS2 were, −0.0308 and −0.0401, respectively, which were well below the value 0.2 value generally accepted for the validity of single reference calculations. It is particularly useful to compare with the recent results of Dixon-Lewis et al.17 for reaction R1 because those calculations employ a higher level ab initio calculation of the barrier energy. Using a full-electron DFT calculation, B3YLP/6-311(d,p), the geometries of the reagents and the saddlepoint were optimized. The vibrational frequencies at this level of theory (with a scaling factor) were used by Dixon-Lewis et al. for zero-point energies and the partition functions. The barrier energy was then calculated at a higher level of theory using the CCSD(T) method and taking a complete basis set extrapolation (CBS) of the cc-pVnZ-PP basis sets. The resulting zero-point-corrected barrier is 10.5 kcal/mol in good agreement with our M11 results. We note, however, that there are some differences between the optimized geometries obtained with B3YLP and those we calculated with M11/cc-pVQZ-PP. The largest differences were seen for the bond lengths to the abstracted hydrogen atom. The B3LYP calculation gave rOH = 1.230 Å and rHBr = 1.637 Å whereas M11 gave rOH = 1.156 Å and rHBr = 1.672 Å. To test whether the agreement between CCSD(T) and M11in the barrier heights was the result of a fortuitous cancelation of errors arising from the differences in geometries, we recomputed the barrier energy from M11 using the geometries obtained using B3LYP. The M11/CBS//B3YLP/ 6-311(d,p) result for the zero-point-corrected barrier was ΔE0 = 10.32 kcal/mol, which agrees well with both M11 and CCSD(T). Our conclusion is that the M11 functional probably provides a reasonably accurate choice for the DFT calculations. The optimized transition state structure for the exchange reaction (R2), i.e., TS2, is found to be the planar symmetrical five-membered ring shown in Figure 2b. The adiabatic barrier for this reaction is predicted by M11/CBS to be 3.56 kcal/mol, which is less than half the barrier height of TS1. The other functionals shown in Table 1 likewise show a dramatic lowering of the barrier, with the exception of B3LYP which, again, is an outlier. The frequencies obtained for TS2 using M11/CBS// M11/cc-pVQZ-PP are presented in Table 2. As might be expected, many of the vibrational frequencies for TS2 are larger than their counterparts in the TS1 structure. In particular, the constrained cyclic geometry of TS2 serves to significantly increase the low-frequency high-amplitude torsional modes and bending modes found in the linear TS1 structure. In addition to TS1 and TS2, we have computed the properties of a doubly hydrogen-bonded intermediate complex. The optimized planar structure for this species is shown in Figure 2c. It is seen that the intermediate has the proper geometry to access TS2, and thus under some circumstances, e.g., very low temperatures, it may serve as a precursor to reaction R2. A binding energy of 3.11 kcal/mol (relative to the reactant asymptote) is predicted by the M11/CBS method. The other functionals listed in Table 1 show comparable results, with B3LYP again yielding the largest deviation. To summarize the influence of isotopic variation, in Figure 3 we show the energy diagram for the various isotopomers of reactions R1 and R2. The zero point effects are seen to have a significant effect on the dynamical barriers for both reactions.

2545.3 1811.8 3585.6, 1407.7, 1214.1 2612.0, 1223.1, 1038.8 3426.4, 2337.1, 1458.75, 1247.3, 462.5, 405.6, 345.9, 197.1, 119.16 3425.8, 1665.7, 1457.6, 1247.1, 404.6, 351.5, 253.3, 186.1, 118.8 2499.6, 2333.4, 1250.6, 1079.8, 462.0, 352.6, 294.0, 188.8, 117.8 2495.5, 1664.5, 1250.3, 1079.3, 347.6, 300.8, 250.1, 178.9, 117.4 3627.3, 1415.7, 1144.2, 1099.4, 821.9, 354.7, 298.0, 139.5, 1867.7i 3626.8, 1410.5, 1119.6, 816.2, 620.8, 357.6, 276.5, 136.8, 1368.6i 2105.7, 1780.8, 1574.1, 1331.1, 1062.3, 994.7, 776.2, 337.8, 127.6i 1951.7, 1519.9, 1331.2, 1262.3, 954.1, 922.5, 609.7, 334.7, 119.4i 1951.7, 1519.9, 1331.2, 1262.3, 954.1, 922.5, 609.7, 334.7, 119.4i 1529.3, 1332.4, 1281.7, 1115.9, 915.6, 737.6, 555.6, 330.8, 109.0i 3709.6, 3707.6, 1417.1, 1311.8, 1004.0, 374.1 3708.9, 2702.4, 1369.7, 1011.0, 994.5, 326.2

a

The symmetry number is 2 for the TS2 structure of HBr + HO2 and DBr + DO2, and is 1 for all other reagents and TS’s.

Figure 2. Configurations of the stationary point for the HBr + HO2 system. Panel a shows the nonplanar TS1 structure for the abstraction reaction (R1). Panel b shows the symmetric planar TS2 for the exchange reaction (R2). Panel c shows the planar intermediate. The bond lengths are given in angstroms and the bond angles are in degrees.

D

DOI: 10.1021/acs.jpca.6b07215 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

the predicted rate coefficient very substantially. The influence of the tunneling corrections is strongest at lower temperatures and can lead to a curved Arrhenius plot. Therefore, the rate coefficients were fit to a sum of two generalized Arrhenius forms, i.e. k(T ) = A1(T /1K )n1 exp( −E A1/kBT ) + A 2 (T /1K )n2 exp( −E A 2 /kBT )

(3.3)

The parameters were set by a nonlinear least-squares fit using a uniform sampling over the full temperature range of 200−1500 K.

III. RESULTS IIIA. HBr + HO2 → H2O2 + Br. We have computed the TST rate coefficient for the abstraction reaction (R1), and its isotopomers, using the M11/CBS level of quantum chemistry for the energetics and M11/cc-pVQZ-PP for the frequencies and geometries. The reaction proceeds through the transition state TS1 depicted in Figure 2a. The adiabatic barrier lies 10.04 kcal/mol above the reactants and 6.09 kcal/mol above the products; the imaginary frequency of the classical barrier is 1867.7 cm−1. These values are in good agreement with the CCSD(T)/CBS//B3LYP results of Dixon-Lewis et al.17 who predicted an adiabatic barrier height of 10.50 kcal/mol. It is not surprising, therefore, that the resulting kTST/W(T) values (conventional TST combined with the Wigner tunneling coefficient) are likewise quite similar. This is illustrated in Figure 4 where we present the Arrhenius plots for k1(T) of

Figure 3. Schematic diagram showing the energetics of the two reaction pathways, HO2 + HBr → H2O2 + Br (TS1) and HO2 + H′Br → H′O2 + HBr (TS2). The zero-point-corrected results for M11/CBS are shown in units of kcal/mol. The isotopic shifts of the energetics due to zero-point energy are displayed.

The isotope effects will be discussed in more detail in section III. IIB. Transition State Theory. The rate coefficients for reactions R1 and R2 were computed using the standard harmonic-oscillator/rigid-rotor approximation to transition state theory (TST). Conventional TST is obtained by assuming the TS lies at the saddlepoint for the computation of the partition functions. To assess the contribution of variational effects and quantum mechanical tunneling to the reaction, it proved useful to compute the minimum energy path (MEP) emanating from the saddlepoint. The MEP was generated using GAMESS47 with the GS2 algorithm with a step size of 0.035 bohr at the M11/cc-pVQZ level of theory. The TST calculations were carried out using the POLYRATE package.49 It was found that that variational effects, due to the shift of the TS away from the saddlepoint, were small and are ignored in this work. The influence of quantum mechanical tunneling is incorporated by multiplying the TST rate coefficient by a quantum transmission coefficient, i.e. k(T ) = k(T ) ·kTST(T )

Figure 4. Rate coefficients for the abstraction reaction HBr + HO2 → H2O2 + Br. The solid lines are theoretical calculations using TST. The red curve is the CCSD(T)/CBS+Wigner tunneling result of DixonLewis.17 The black curve is the M11/CBS+SCT result, the green curve is the M11/CBS result with no tunneling, the aqua curve is M11/CBS +Wigner tunneling, and the purple curve is M11/CBS+ZCT. The experimental results, shown with symbols, represent the various “upper limits” obtained for the rate coefficient.8−10,15,16

(3.1)

The simplest tunneling model is the Wigner expression, i.e. κ (T ) = 1 +

2 1 ⎛ ℏω ⎞ ⎜ ⎟ 24 ⎝ kBT ⎠

(3.2)

where ω is the magnitude of the imaginary barrier frequency evaluated at the saddlepoint. The Wigner tunneling expression is popular because of its simplicity. Because eq 3.2 requires merely a normal-mode analysis of the TS, the complexities of locating tunneling paths can be ignored. However, it is widely appreciated that the Wigner tunneling coefficient often shows large error because it does not include effects such vibrational adiabaticity and dynamical corner cutting.50 The zero curvature tunneling method49 (ZCT) is a more accurate theory because it incorporates tunneling through the exact adiabatic potential barrier. The SCT (small curvature tunneling) theory,51 which is a modification of the earlier SCSAG concept,52 is still more reliable because it takes into account the corner cutting effect. We shall find that more accurate tunneling models can increase

reaction R1 obtained in various ways. The CCSD(T)/CBS// B3LYP result (red curve) and the present M11/CBS//M11/ccpVQZ-PP with Wigner tunneling result (aqua curve) are in good agreement down to about 300 K. The results at high temperatures are close, e.g., kTST/W(T=1000K) = 1.31 × 10−14 for M11 and 9.05 × 10−15 for CCSD(T). The ZCT results (purple curve) are seen to be similar to Wigner tunneling down to 300 K and then become somewhat larger than Wigner tunneling. As we noted earlier, a more accurate representation of the tunneling coefficient is provided by the SCT method. The k1(T) rate coefficient for SCT computed with the M11/ccpVQZ-PP reaction path is depicted with the black curve. It is E

DOI: 10.1021/acs.jpca.6b07215 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A seen that at room temperatures the rate coefficient is roughly 1 order of magnitude larger when the SCT method was used than was obtained with Wigner tunneling. At high temperatures, the tunneling contribution to the reaction is small and the results for SCT are close to those obtained using ZCT. The best fit parameters obtained for a nonlinear least-squares optimization of a two exponential form for k1(T) are presented in Table 4. Table 3. Energies for the Stationary Points of the Isotopes of the HBr + HO2 System Computed Using M11/CBS including Harmonic Zero-Point Energy ΔE0 (kcal/mol) HO2 + HBr

HO2 + DBr

DO2 + HBr

DO2 + DBr

0 −3.11 10.04 3.56 3.95 0

−1.05 −4.39 9.28 2.02 1.93 −1.91

−1.91 −5.15 N/A 2.02 N/A −1.05

−2.96 −6.43 N/A 0.47 N/A −2.96

reactants intermediate TS1 TS2 products R1 products R2

Figure 5. Rate coefficient for two isotopomers of reaction R1. The k1(T) is computed using M11/CBS barrier energies and the SCT tunneling coefficient obtained from the M11/cc-pVQZ-PP reaction path.

presented in Table 1 that the adiabatic barrier energy is much lower for reaction R2 than it is for reaction R1. Using M11/ CBS, we find ΔE0 = 3.56 kcal/mol for TS2 compared with ΔE0 = 10.04 kcal/mol for TS1. The TS2 structure is seen to bear a strong resemblance to the planar intermediate shown in Figure 2c. The hydrogen bonding in the intermediate is seen to be fairly weak, ΔE0 = 3.11 kcal/mol. This can be compared with other doubly hydrogen-bonded systems such as ΔE0 = 6.1 kcal/ mol for H2O·HO2,53 and 9.8 kcal/mol for (HO2)2.26 Hence, we suspect that the TS stabilization likewise may not be as strong for HBr + HO2 as seen in other reactions. Indeed, the barrier for the HO2 + HO2 reaction is actually submerged; i.e., it lies below the energy of the reactants. Furthermore, the fivemembered ring of TS2 is clearly strained with a bond angle of only 53.1°, which may also tend to lessen the stabilization. Because the complex intermediate species is only stabilized by 3.11 kcal/mol, we have computed the rate coefficient using only positive collision energies, i.e., corresponding to unbound HBr and HO2 reagents. The TST rate coefficient for reaction R2 was computed using ZCT theory and was fit to the two exponential form with parameters given in Table 4. In Figure 6 we plot the rate coefficient k2 versus temperature for reaction R2 along with k1 obtained previously for reaction R1. We see that near T = 300 K the reaction R2 rate coefficient is roughly 2−3 orders of magnitude larger than that for reaction R1 even including the quantum tunneling effects. The difference is due in most part to the large disparity in the barrier heights. However, at higher combustion temperatures near T = 1000 K, the rate coefficients for reactions R1 and R2 actually become roughly comparable. To understand this behavior, we need to consider factors other than the activation energy that may influence them.

We also show in Figure 4 some of the experimental results that have been obtained for reaction R1. As noted above, most experiments were carried out at lower temperatures (