Experimental and Modeling Studies of the Pressure and Temperature

Jan 14, 2011 - Electronic structure calculations of the stationary points of the potential energy surface were used with a master equation analysis to...
0 downloads 10 Views 2MB Size
Download PDF
ARTICLE pubs.acs.org/JPCA

Experimental and Modeling Studies of the Pressure and Temperature Dependences of the Kinetics and the OH Yields in the Acetyl þ O2 Reaction Scott A. Carr,† David R. Glowacki,*,§ Chi-Hsiu Liang,† M. Teresa Baeza-Romero,‡ Mark A. Blitz,† Michael J. Pilling,*,† and Paul W. Seakins† †

School of Chemistry and ‡School of Earth and Environment, University of Leeds, Leeds LS2 9JT, United Kingdom

bS Supporting Information ABSTRACT: The acetyl þ O2 reaction has been studied by observing the time dependence of OH by laser-induced fluorescence (LIF) and by electronic structure/master equation analysis. The experimental OH time profiles were analyzed to obtain the kinetics of the acetyl þ O2 reaction and the relative OH yields over the temperature range of 213-500 K in helium at pressures in the range of 5-600 Torr. More limited measurements were made in N2 and for CD3CO þ O2. The relative OH yields were converted into absolute yields by assuming that the OH yield at zero pressure is unity. Electronic structure calculations of the stationary points of the potential energy surface were used with a master equation analysis to fit the experimental data in He using the high-pressure limiting rate coefficient for the reaction, k¥(T), and the energy transfer parameter, ÆΔEdæ, as variable parameters. The best-fit parameters obtained are k¥ = 6.2  10-12 cm-3 molecule-1 s-1, independent of temperature over the experimental range, and ÆΔEdæ(He) = 160(T/298 K) cm-1. The fits in N2, using the same k¥(T), gave ÆΔEdæ(N2) = 270(T/298 K) cm-1. The rate coefficients for formation of OH and CH3C(O)O2 are provided in parametrized form, based on modified Troe expressions, from the best-fit master equation calculations, over the pressure and temperature ranges of 1 e p/Torr e 1.5  105 and 200 e T/K e 1000 for He and N2 as the bath gas. The minor channels, leading to HO2 þ CH2CO and CH2C(O)OOH, generally have yields 99.9%), acetone-d6 (Euriso-Top % D > 99.8%), and TBH (70% aq) were degassed prior to use. He (BOC, CP grade,99.999%), N2 (air products, 99.999%), and O2 (Air products, high purity, 99.999%) were used straight from the cylinder.

3. ELECTRONIC STRUCTURE CALCULATIONS 3.1. single-reference Methods. The goal of the G3X family of methods is to accurately approximate the QCISD(T)/G3XL model chemistry, where G3XL refers to the G3XLarge basis set.27,28 This method has been shown to give accurate energies on those molecules featured in the G3/99 test set. In the G3X methods, geometry optimizations are carried out at the B3LYP/6-311þG(3df,2p) model chemistry (alternatively, a smaller 6-31G(2df,p) Pople basis has been found to give energies that are essentially identical to those obtained with the larger Pople basis27). 1071

dx.doi.org/10.1021/jp1099199 |J. Phys. Chem. A 2011, 115, 1069–1085

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Calculated stationary points (G3X(MP2)-RAD) for acetyl þ O2. ZPE corrections have been made.

In the G3X method, a MP4/6-31(d) single-point energy calculation is carried out on the DFT determined geometry, and a series of corrective single-point energy calculations are subsequently performed. For molecules that feature atoms no larger than those in the second row, the G3X energy is determined as follows E0 ðG3XÞ ¼ MP4=6-31GðdÞ þ ΔEðCCÞ þ ΔEð þ Þ þ ΔEð2df ; pÞ þ ΔEðG3XLargeÞ þ EðHLCÞ þ EðZPEÞ ðE1Þ where ΔE(CC)29 is a correction for correlation effects beyond fourth-order perturbation theory, ΔE(þ) is a correction for diffuse functions, ΔE(2df,p) is a correction for higher polarization functions, and ΔE(G3XLarge) is a correction for larger basis set effects than 6-31G(d). E(ZPE) is the ZPE correction, and E(HLC) is a higher-level correction, which is essentially a semiempirical parameter30 that has been optimized to give the best possible agreement between the calculated molecular properties and the experimentally measured values featured in the G3/99 test set. The G3X(MP2) method is a less computationally expensive implementation of the G3X method, wherein the additivity equation describing the energy of a molecular geometry is E0 ½G3XðMP2Þ ¼ MP4=6-31GðdÞ þ ΔEðCCÞ þ ΔEðMP2Þ þ ΔEðG3XLargeÞ þ EðHLCÞ þ EðZPEÞ

ðE2Þ

where all of the elements of this equation are as described above and ΔE(MP2) is a correction for the differences introduced by the much larger G3MP2L basis set with respect to the 6-31G(d) basis. For open-shell species, the G3X(MP2) protocol indicates that the post Hartree-Fock (HF) calculations in eq E2 feature an unrestricted HF (UHF) calculation. While UHF wave functions allow electron delocalization with respect to the ROHF

wave function, they do not rigorously conserve spin symmetry and may be contaminated by higher spin states, which can affect the calculated energies,31 as noted by Hou et al.22 in their calculations on this system. The G3X(MP2)-RAD model chemistry is a variant of G3X(MP2) theory designed for open-shell species in those cases where spin contamination in the UHF G3X(MP2) wave function suggests that the calculated energies may not be reliable.29 For all post-SCF calculations, the G3X(MP2)-RAD method replaces the UHF wave function with a restricted open-shell HF (ROHF) calculation. The G3X(MP2)RAD method replaces the QCISD(T) calculation with a UCCSDROHF calculation, implemented herein using MOLPRO.32 All other single-reference calculations were carried out using G03.28 Figure 1 shows the stationary points calculated using the G3X and G3X(MP2)-RAD methods, and Table 1 gives the corresponding energies and spin expectation values, ÆS2æ, for these stationary points. Apart from the transition states TS3 and TS4 for dissociation of Int2, all of the stationary points have G3X and G3X(MP2)-RAD energies that agree within ∼3 kJ mol-1. Given that the T1 and D1 diagnostics for TS3 and TS4 are larger than those for any of the other species shown in Figure 1, the difference in energy between the two model chemistries likely derives from multireference effects in this region of the PES. The spin expectation value for TS3 indicates contamination of a lowlying quartet state, and multireference calculations suggest that the energy of the TS4 ground-state wave function is affected by a low-lying excited state, which is observed in both C1 and Cs(A00 ) symmetry. Despite the differences in calculated energies for these two species, they do not have a significant impact on the kinetics and yields simulated in the ME model described below. As discussed below, stabilization tends to occur only in the acetyl peroxy radical (Int1) well. Even if the TS3 energy is set to its maximum calculated value at the G3X model chemistry, stabilization 1072

dx.doi.org/10.1021/jp1099199 |J. Phys. Chem. A 2011, 115, 1069–1085

The Journal of Physical Chemistry A

ARTICLE

Table 2. Lennard-Jones Parameters Used in the Calculation of the Collision Frequency in the ME Calculations species Int1 and Int2 He N2

Figure 2. Geometry and atomic labeling of TS1.

of the hydroperoxide isomer (Int2) is unimportant; at the lowest temperatures and highest pressures considered in this study, peroxy species that cross TS1 react via TS3. Similarly, reaction via TS4 is unable to compete with reaction via TS3; even if TS4 is set to its lowest computed energy at the G3X(MP2)-RAD model chemistry, the ketene yield is less than 0.01, which is in agreement with previous experimental measurements that found minimal ketene yields.15,25 Also shown in Table 1 are the stationary point energies calculated by both Hou et al.24 and Maranzana et al.25 Our calculated DFT energies are in good agreement with those calculated by Maranzana et al., although there are some significant differences between the DFT and G3X calculations. Typical of the B3LYP functional, several of the transition-state energies are underestimated.33 Our G3X energies generally agree well with those found by Hou et al. The biggest difference between their G3X energies and ours concerns TS4 and is likely due to the low-lying excited state discussed above. With respect to the ME kinetics and yields that we describe below, the most significant discrepancy between the DFT and wave function energies pertains to the relative energy of acetyl þ O2 with respect to that of the acetyl peroxy radical. Whereas the G3X methodology calculates the acetyl-O2 binding energy to be more than 140 kJ mol-1, the DFT calculations give a result of 130 kJ mol-1. The consequences of this difference are discussed further below. The DFT geometries obtained in this work are given in the Supporting Information. 3.2. Multireference Methods. In addition to the singlereference electronic structure calculations described above, we utilized two different active spaces to carry out multireference calculations on two different parts of the PES shown in Figure 1, (1) Int1, TS1, Int2, TS3, and 3C2H2O2 þ OH and (2) the acetyl þ O2 entrance channel. In order to guarantee a consistent active space for the CAS calculations with respect to different molecular geometries, isosurfaces of the aforementioned molecular orbitals were visualized with MOLDEN.34 Figure 2 shows the geometry of TS1, and the following discussion of active spaces used in the multireference calculations refers to the numbering scheme shown in that figure. For Int1, TS1, Int2, TS3, and 3C2H2O2 þ OH, the CASSCF reference wave function was computed in order to investigate static electron correlation on the SCF energies using an 11 electron, 11 orbital (11,11) active space. Precaution was taken to maintain the consistency of the active space for all isomers.

ε/K

σ/A

473

5.1

10.2 48.0

2.6 3.9

In Int1, the active space included the σ and σ* C3-H4 orbitals, the σ and σ* C1-O2 orbitals, the σ and σ* C1-O7 orbitals, the σ and σ* O7-O8 orbitals, the π and π* C1-O2 orbitals, and the singly occupied molecular orbital (SOMO), which is localized on C8 in Int1. Geometry optimizations and frequencies for each stationary point were calculated at the CAS(11,11)/cc-PVTZ model chemistry. Inspection of the orbital occupation numbers obtained upon diagonalization of the electron density matrix showed that all of the orbitals within the (11,11) active space had occupation numbers between 0.02 and 1.98. The effects of dynamic electron correlation on the CASSCF energies were computed using second-order multireference Rayleigh-Schrodinger perturbation theory (CASPT2), as implemented in MOLPRO.35,36 Correlation within the active space was calculated by considering all possible excitations within the active space and single and double excitations within the valence orbitals. Core electrons were frozen and not correlated with the other electrons in the molecule. Following previous multireference studies on O2 addition to radicals,25,37 we treated the acetyl þ O2 entrance channel with a seven electron, five orbital active space that included all of the π and π* electrons on O2 and the radical orbital on the acetyl radical. A relaxed PES scan along the C1-O7 bond distance was carried out at the CAS(7,5)/cc-pVTZ model chemistry, and CASPT2 calculations using the cc-pVTZ basis were subsequently performed. The CASPT2/cc-pVTZ results showed a small barrier at 2.5 Å; however, previous work by Klippenstein and co-workers37 on ethyl þ O2 association has shown that CASPT2 maxima which are observed using smaller basis sets disappear with a complete basis set (CBS) extrapolation scheme. To check whether or not this was similarly the case for the acetyl þ O2 PES, we performed subsequent CASPT2/cc-pVQZ calculations at each CAS(7,5)/cc-pVTZ optimized geometry. Following from the work of Martin38 as well as Feller and Dixon,39 the energies at the infinite basis set limit, ECBS, were extrapolated using the relation ECBS ¼ Elmax -

B ðlmax þ1Þ4

ðE3Þ

where lmax is the highest angular momentum component in the basis set (lmax = 3 and 4 for cc-pVTZ and cc-pVQZ, respectively), Elmax is the corresponding energy, and B is a constant determined by fitting. Figure 3 shows the results of CASPT2 relaxed PES scans on the acetyl þ O2 entrance channel. At the CASPT2/cc-pVTZ model chemistry, there is a van der Waals well with a depth of 1.86 kJ mol-1 at a separation of ∼3.1 Å and a barrier of 0.3 kJ mol-1 above the reaction asymptote at 2.5 Å. With a cc-pVQZ basis, the van der Waals well depth is 1.3 kJ mol-1; the barrier is 0.3 kJ mol-1 and occurs at a separation of 2.55 Å. The CBS surface shows the initial well to be shallower, with a depth of 0.95 kJ mol-1, and the activation barrier, which occurs around 2.6 Å, is submerged by 0.24 kJ mol-1 below the reaction asymptote. Transition-state searches initialized using a CASSCF and CASPT2 model chemistry with a cc-pVTZ basis set were unable to locate 1073

dx.doi.org/10.1021/jp1099199 |J. Phys. Chem. A 2011, 115, 1069–1085

The Journal of Physical Chemistry A

ARTICLE

Table 3. Molecular Parameters Corresponding to Figure 1 (ZPE(H) - ZPE(D))/

rotational constantsc/cm-1

frequencies (nondeuterated)c/cm-1

kJ mol-1

CH3CO

0.3155, 0.333, 2.85

111, 470, 853, 956, 1050, 1359, 1454, 1460, 1930, 3014, 3105, 3111

O2

1.456

1645

Int1

0.108, 0.161, 0.311

139, 167, 325, 507, 543, 550, 734, 982, 1048, 1130, 1195, 1402, 1463, 1468, 1894, 3060, 3120, 3161

26.77

TS1

0.113, 0.163, 0.343

1688i, 154, 307, 490, 545, 604, 685, 844, 947, 1035, 1038, 1106, 1149, 1420, 1723, 1865, 3109, 3192

21.98

TS2

0.090, 0.125, 0.305

1063i, 59, 211, 351, 472, 493, 550, 657, 886, 1028, 1090, 1276, 1311, 1389, 1702, 2182, 3091, 3165

22.31

Int2

0.109, 0.153, 0.377

206, 323, 337, 426, 500, 639, 684, 768, 873, 1008, 1030, 1248, 1451, 1502, 1710, 3170, 3290, 3474

26.14

TS3 TS4

0.100, 0.142, 0.316 0.327, 0.142, 0.099

656i, 157, 188, 296, 320, 378, 616, 634, 820, 973, 1009, 1013, 1183, 1430, 1860, 3156, 3276, 3683 369i, 150, 260, 398, 471, 550, 590, 685, 765, 980, 1128, 1174, 1436, 1513, 2009, 3136, 3173, 3286

25.34

OH lactone OH C2H2O2a

3

HO2 ketene

18.753

25.72 0.00

3724

0.842, 0.271, 0.214 18.753

494, 533, 719, 946, 1003, 1066, 1121, 1194, 1484, 2001, 3116, 3217 3724

0.406, 0.364, 0.192 20.648, 1.131, 1.072b 9.478, 0.344, 0.332b

441, 462, 602, 620, 706, 857, 1007, 1146, 1418, 1572, 3339, 3466 1184, 1449, 3580b 450, 552, 593, 994, 1176, 1414, 2246, 3194, 3290b

a

Rotational constants and frequencies obtained using a CASSCF(8,8) model chemistry. b UB3LYP/6-31G(2df,p). c UB3LYP/6-31G(3df,2p). Differences in ZPEs for the H and D systems are only given for the significant species in the calculations.

Figure 3. Multireference results for the acetyl þ O2 entrance channel. The results shown are from a relaxed PES scan along the nascent C-O bond length. The energies are relative to the acetyl þ O2 reaction asymptote.

a saddle point, presumably due to the extremely slight curvature shown in Figure 3.

4. MASTER EQUATION CALCULATIONS The form of the energy grained ME used in this work has been described in detail previously.40-42 The phase space for each isomer was divided into energy grains with a width of 100 cm-1. Collisional energy transfer in the grained phase space and interconversion between species were described with a set of coupled differential equations. Briefly, the form used in this work is d p ¼ Mp dt

ðE4Þ

where p is a vector containing the populations, ni(E), of the energy grains, where i refers to the ith isomer, and M is a matrix that determines grain population evolution due to collisional energy transfer and reaction. Collisional energy transfer was

described using an exponential down model, parametrized with the appropriate Lennard-Jones parameters given in Table 2, and the average energy transferred in a downward direction, ÆΔEdæ. The ÆΔEdæ parameter was determined using a least-squares fit to the experimental data, as discussed in the following section. As input for the ME calculations, we included stationary points for the entire PES shown in Figure 1, with the exception of the channel leading to 3C2H2O2, for which the yield is negligible over the conditions considered in this study. We investigated the sensitivity of the ME results to each of the stationary point energies depicted in Figure 1. Neither the kinetics nor the yields that we report below significantly depend on whether the G3X, G3X(MP2)-RAD, or CASPT2 results in Table 1 were used to model the energies of Int2, TS3, TS4, and the dissociation channels. The yields and kinetics did, however, show significant sensitivity to the method used for modeling the relative energies of acetyl þ O2, TS1, and TS2 with respect to Int1. To model the relative energies of TS1 and TS2, we used the CASPT2 energies given in Table 1. As we were unable to obtain a reliable CASPT2 energy for acetyl þ O2 with respect to Int1 due to difficulties in maintaining a consistent active space, we modeled the acetyl þ O2 energy by taking the average of the G3X and G3X(MP2) values in Table 1. Microcanonical rate coefficients for isomerization and association were calculated using vibrational harmonic oscillator approximations and by treating all molecules as either two- or three-dimensional classical free rotors. The relevant molecular parameters are given in Table 3. The k(E)s for the hydrogen-transfer reaction of TS1 were corrected for tunnelling effects by convoluting the classical RRKM k(E)s with tunnelling probabilities calculated assuming a one-dimensional asymmetric Eckart potential.43 All of the unimolecular and isomerization reactions shown in Figure 1 were treated using rovibronic state densities to calculate the RRKM k(E)s. For the barrierless acetyl þ O2 reaction, conventional TST and variational TST are not appropriate because the reaction coordinate is not well-defined, a situation which is typical of association reactions between openshell species.44 We have calculated the association k(E)s using a more pragmatic approach, by taking the inverse Laplace transform (ILT) of the modified Arrhenius form of the experimental 1074

dx.doi.org/10.1021/jp1099199 |J. Phys. Chem. A 2011, 115, 1069–1085

The Journal of Physical Chemistry A

ARTICLE

high-pressure association rate coefficient to derive the sum of states, W(E), for the association transition state.45 Briefly, if QR(β) represents the rovibronic partition function for the reactants in the association reaction, k1¥(β) represents the high-pressure limiting association rate coefficient and k-1¥(β) represents the high-pressure limiting dissociation rate ¥ (β)/k1¥(β), then the coefficient, where β = 1/kbT, and Keq(β) = k-1 mathematical relationship for deriving W(E) via ILT is as follows WðEÞ ¼ hL - 1 fk1 ¥ ðβÞKeq ðβÞQR ðβÞg

ðE5Þ

where L-1{f} indicates the ILT of the function f. The ILT in eq E5 has an exact analytical solution if Keq is rewritten in terms of the reactant and product partition functions and k1¥(T) can be represented using the modified Arrhenius form  n T k1 ¥ ðTÞ ¼ A expð - Ea =ðkb TÞÞ ðE6Þ T0 In this work, T0 was set to 300 K. Ea was set to 0 kJ mol-1 based on the CASPT2/CBS results in Figure 3, which suggest a barrierless association process. Both A and n were left as variable parameters and were determined according to a least-squares fit to the experimental data. The fitting methodology is discussed in further detail below. For each ME calculation performed in this work, the discretized matrix M was diagonalized, and the eigenpairs were determined to give a solution of the form p ¼ Ueλt U - 1 pð0Þ

Figure 4. A typical kinetic trace of the OH LIF signal following 248 nm acetone photolysis; [O2] = 1.1  1016 molecules cm-3, [acetone] = 6  1013 molecules cm-3, He bath gas pressure = 25 Torr, and T = 298 K. The line through the data is a nonlinear least-squares fit of eq E8 to the data. The lower trace shows the residuals of the fit.

ðE7Þ

where p(0) contains the initial conditions (i.e., t = 0) for each grain (i.e., niE(0)), U is a matrix of eigenvectors obtained from diagonalization of M, and λ is a vector of the corresponding eigenvalues. Within our formulation of the ME, bimolecular products are represented using the infinite sink approximation. To analyze the available experimental kinetic data, the microcanonical information contained in the ME solution of eq E5 must be transformed to phenomenological rate coefficients, and the procedure that we use for doing so is similar to that described by Bartis and Widom.40,46,47 All ME calculations and the subsequent eigenvalue-eigenvector analysis were carried out using our recently developed MESMER program, which is freely available on the Web.26

5. EXPERIMENTAL RESULTS 5.1. OH Kinetics. Following 248 nm acetone photolysis, the kinetics of OH production is controlled by the reaction of CH3CO þ O2 (reactions R1 and R2) and by other loss processes of CH3CO. At higher temperatures, the latter are dominated by dissociation to form CH3 þ CO (reaction R3). OH is also slowly lost from the system by reaction with acetone and impurities and diffusion out of the probed region. Solving the rate equations for this reaction set gives the timedependent OH concentration, eq E86

½OHt ¼

R½CH3 CO0 k 0 ðexpð - kl tÞ - expð - k 0 tÞ ðk 0 - kl Þ

ðE8Þ

where R = k2/(k1 P þ k2), k = k1 þ k2, k0 = k[O2] þ k3, kl ¼ k½acetone þ i ki ½impurityi þ kdiff , and [CH3CO]0 is the initial acetyl concentration. An analogous reaction scheme was used to model OD production following photolysis of acetone-d6. Equation E8 was fitted to the observed LIF signal

Figure 5. Bimolecular plots of k0 versus [O2]. 9: 298 K, 25 Torr He; 4: 100 Torr He. Solid lines are weighted least-squares regression fits. Dashed lines show the 95% confidence limits at 500 K.

with k0 and kl as fitting parameters; a typical fit and set of residuals are shown in Figure 4. Typically, partial pressures were