A Theoretical Study on Reactions of Alkylperoxy Radicals

Subscriber access provided by UNIV OF LOUISIANA

A: Kinetics, Dynamics, Photochemistry, and Excited States

A Theoretical Study on Reactions of Alkylperoxy Radicals Yanjiao Xu, Shuanghui Xi, Fan Wang, and Xiang-Yuan Li J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.9b01496 • Publication Date (Web): 16 Apr 2019 Downloaded from http://pubs.acs.org on April 16, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1 2 3

4

A Theoretical Study on Reactions of Alkylperoxy

5

Radicals

6 7

Yanjiao Xu 1, Shuanghui Xi 2, Fan Wang 1,* and Xiangyuan Li 2

8

9 10 11 12

13

1 Institute

of Atomic and Molecular Physics, Key Laboratory of High Energy Density Physics

and Technology, Ministry of Education, Sichuan University, Chengdu, P. R. China. 2 College

of Chemical Engineering, Sichuan University, Chengdu, P. R. China.

14 15 16 17 18 19

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 39

1

ABSTRACT

2

We carried out theoretical study on geometries, relative energies of stationary points and reaction

3

rate constants for ethyl + O2, propyl + O2 and butyl + O2 reactions, which are important reactions

4

in low temperature oxidation of corresponding alkanes. Geometries with CCSD(T)/aug-cc-pVTZ

5

for the ethyl + O2 system are adopted as benchmark to choose a proper exchange-correlation

6

functional for geometry optimization. Our results show that B3LYP with 6-311+G(d,p) can

7

provide reliable structures for this system and structures of the other two systems are determined

8

with this functional. Performance of the explicitly correlated CCSD(T)-F12a and the locally

9

correlated DLPNO-CCSD(T) methods on barrier heights and reaction energies are evaluated by

10

comparing their results with those of CCSD(T)/aug-cc-pVQZ for the ethyl + O2 system. Our

11

results indicate that reliable energy differences for this system are achieved with CCSD(T)-F12a

12

using the cc-pVDZ-F12 basis set and this method is employed in calculating single point

13

energies for the other two systems. The single-reference equation-of-motion spin-flip coupled-

14

cluster method is adopted to obtain potential energy surface of the barrierless reaction C2H5• +

15

O2  CH3CH2OO• and the results are compared with those using broken-symmetry DFT and the

16

Morse potential. Differences between energies with these methods are less than 1.6 kcal/mol, but

17

difference in the rate constants could be sizeable at temperatures lower than 500K and rate

18

constants obtained in this work are reliable only for temperatures higher than 500K. Pressure-

19

dependent rate constants for these reactions are determined using the Rice-Ramsperger-Kassel-

20

Marcus/Master equation method. The obtained reaction energies, barrier heights and rate

21

constants could be valuable for reactions between large alkane radical and O2, which is important

22

in low temperature combustion of fuels such as kerosene and gasoline.


2

Page 3 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1


I. INTRODUCTION

2

Detailed chemical kinetic models for combustion of alkanes have been extensively studied and

3

reasonable results at temperatures higher than 1000K can usually be achieved.1–5 Nevertheless,

4

there still exit some uncertainties in the reaction mechanisms at low-temperatures (500-1000K).6

5

Compared with high temperature mechanisms, much more elementary reactions and chemical

6

species are involved in the low-temperature combustion mechanisms of alkanes. The radical

7

species (•R) could react rapidly with oxygen molecules to form alkylperoxy radicals (ROO•),

8

which can undergo further reactions at low temperatures. ROO• could form carbon-centered

9

hydroperoxy alkyl radical (•QOOH) via intermolecular hydrogen abstraction. Both alkylperoxy

10

radicals and hydroperoxy alkyl radicals could produce bimolecular products such as hydroperoxy

11

radicals and oleﬁns or cyclic alkanes, hydroxyl radicals with aldehydes, ketones, or cyclic ethers.

12

On the other hand, a hydroperoxy alkyl radical could further react with an additional oxygen

13

molecule to produce hydroperoxyl-alkyl-peroxyl (•O2QOOH) radicals. These elementary

14

reactions become important in the low-temperature oxidation of alkanes.4,7

15

In recent years, some intermediates in the low-temperature combustion pathways of alkanes

16

have been detected experimentally and the low-temperature oxidation scheme of alkanes is

17

confirmed.8-13 However, it is still difficult to obtain reaction rate constants involving peroxy

18

radicals and hydroperoxy radicals experimentally.7,14,15 Therefore, rate constants of these low-

19

temperature combustion reactions are mainly estimated or calculated with theoretical

20

methods.7,15 To obtain reasonable reaction rate constants, reliable reaction energies, barrier

21

heights, geometric structures of reactants, products and transition states and the corresponding

22

vibration frequencies are required. High-accurate quantum chemistry calculations are thus


3


Page 4 of 39

1

needed. This paper focuses on quantum chemical study of reactions involving alkylperoxy

2

radicals and hydroperoxy alkyl radicals.

3

Reactions between ethyl + O2, propyl + O2 and butyl + O2 are important prototypes of the

4

alkyl+ O2 reaction class, and their rate constants provide insights on reaction pathways and

5

reaction rate of larger alkyl systems.16 In the overwhelming majority of cases, DFT method17 was

6

applied to the calculation of geometries and vibration frequencies for reactants, products and

7

transition states, due to its compromise between accuracy and efficiency. In previous

8

studies14,15,18-28 on ethyl + O2, propyl + O2 and butyl + O2 reactions, geometries and frequencies

9

are almost always determined with DFT using B3LYP29 or M06-2X30. Despite its increasing

10

popularity, the employed exchange-correlation (XC) functional has a pronounced effect on

11

results. On the other hand, Schaefer et al.31,32 and Klippenstein et al. 16,33 optimized geometries of

12

stationary points in the ethyl+O2 reaction using the CCSD(T) method,34 which is the “gold

13

standard” of quantum chemistry. Schaefer et al. 32 only optimized geometries of molecules with

14

certain symmetry at the CCSD(T)/cc-pVQZ level in their work. Klippenstein et al. 16,33 employed

15

CCSD(T)/cc-pVTZ to optimize geometries in the ethyl + O2 system. However, diffuse basis

16

functions were not included in these previous CCSD(T) calculations16,32,33 for the ethyl + O2

17

system. Transition state structures were generally loose and diffuse basis functions may have

18

certain effect on the result. Furthermore, CCSD(T) is computationally intensive, and DFT

19

calculations are still required for structures of larger systems. Structures with CCSD(T) can be

20

employed to choose an appropriate XC functional for DFT calculations.

21

Reaction barrier height usually has a significant effect on reaction rate constants and accuracy

22

of DFT results on reaction barrier heights is generally insufficient. Single point energies are

23

mostly calculated using CCSD(T) or QCISD(T).14-16,18-26,31-33 In addition, a large basis set is


4

Page 5 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

needed in such calculation to achieve results with high accuracy. Scaling of the CCSD(T)

2

method is N7, where N represents system size, and it can therefore only be used for very small

3

molecules. For the ethyl + O2 system, Schaefer et al.32 employed a high-level focal point analysis

4

with basis set extrapolation to calculate single-point energies for structures with certain

5

symmetry. Klippenstein et al.16,33 adopted the ANL0 (CCSD(T)/CBS (aug-cc-pVQZ, aug-cc-

6

pV5Z)) to calculate energies of stationary point for the C2H5• + O2 system. QCISD(T)/CBS (cc-

7

pVTZ and cc-pVQZ) is used by Goldsmith et al.26 in computing energies of stationary points in

8

the C3H7• + O2 system. On the other hand, the explicitly correlated CCSD(T)-F12x(x=a,b)

9

method35,36 has been developed, where the wave function depends explicitly on the distance

10

between two electrons. Computational cost of the CCSD(T)-F12x method is only slightly larger

11

than that of the conventional CCSD(T) method when using the same basis set. However, the

12

CCSD(T)-F12x method is not as sensitive to the choice of the basis set, and highly-accurate

13

results can be obtained by using a small basis set.36 Therefore, the CCSD(T)-F12x method can

14

possibly be applied to somewhat larger systems. On the other hand, the recently developed

15

locally correlated coupled-cluster methods such as the DLPNO-CCSD(T) method37 have a

16

significantly lower computational cost than the conventional CCSD(T) method. DLPNO-

17

CCSD(T) has been reported to be applied to single-point calculations in transition metal

18

complexes with more than 100 atoms.38 However, its accuracy needs further investigation.

19

The reaction alkyl + O2  ROO• is barrierless, and a reasonable potential energy surface

20

(PES) is required to calculate its reaction rate constant using variational transition state theory39

21

(VTST). Multi-reference methods are usually adopted to calculate the PES. CASPT240 is one of

22

the most popular methods in such calculations. However, CASPT2 results depend on the chosen

23

active space. This method would be expensive when the active space is too large. Minimum


5


Page 6 of 39

1

active space was usually employed in previous work,16,26 which would reduce accuracy of the

2

obtained results. In practice, the broken-symmetry DFT method was adopted,41-43 but significant

3

spin contamination was introduced. On the other hand, the equation-of-motion spin-flip coupled-

4

cluster (EOM-SF-CC) method44 is a single-reference method, which effectively avoids the

5

selection of the active space. In EOM-SF-CC calculations, a quadruple-spin state is chosen as

6

reference, which can be described reliably by a single reference method. The doublet-spin state

7

can be obtained from the reference quadruple-spin state by a spin-flip transition using EOM-SF-

8

CCSD.44 A reliable PES can thus be achieved with the EOM-SF-CCSD method.

9

In the present work, we study the key chemical reactions of ROO• and •QOOH radicals

10

involved in the low-temperature combustion for ethyl + O2, propyl + O2 and butyl + O2 systems.

11

The following types of reactions are investigated:

12

R•+O2  ROO•,

(i)

13

R•+ O2  olefin + HO2,

(ii)

14

ROO•  •QOOH,

(iii)

15

ROO•  aldehydes, ketones, or cyRO + OH,

(iv)

16

ROO•  olefin + HO2,

(v)

17

•QOOH  cyRO + OH,

(vi)

18

•QOOH  olefin + HO2

(vii)

19

•QOOH  olefin + formaldehyde + OH

(viii)

20

where R• is an alkyl radical, cyRO is a cyclic ether, ROO• is an alkylperoxy radical and •QOOH

21

is an hydroperoxy alkyl radical. Geometries for stationary points in the C2H5• + O2 system will

22

be obtained at the CCSD(T)/aug-cc-pVTZ level and these geometries will be used as benchmark

23

to choose an appropriate XC functional for propyl + O2 and butyl + O2 systems. Secondly, the


6

Page 7 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

corresponding single point energies will be calculated using the explicitly correlated CCSD(T)-

2

F12a method and accuracy of the DLPNO-CCSD(T) method on barrier heights will be

3

investigated. A proper approach will be chosen to calculate single point energies for these

4

systems. The PES of the barrierless reaction C2H5• + O2  CH3CH2OO• will be calculated with

5

the EOM-SF-CC method. However, the EOM-SF-CC method with a reasonable basis set is still

6

computationally expensive, and it is difficult to be applied to propyl + O2 and butyl + O2

7

systems. Results of EOM-SF-CC methods for C2H5• + O2  CH3CH2OO• can be employed to

8

evaluate performance of other approximate approaches such as a Morse potential and broken-

9

symmetry DFT method.

10

This paper is organized in the following manner: computational details are described in Sec. II.

11

In Sec. III, an appropriate XC functional is chosen for geometry optimization of stationary points

12

in the alkyl + O2 systems, as well as for calculation of frequencies. Method for high-accurate

13

single-point energy calculations will also be chosen. The PES for the barrierless reaction C2H5• +

14

O2  CH3CH2OO• will be obtained with the EOM-SF-CC method, and reliability of other

15

approximate approaches will be evaluated. PESs of C3H7• + O2  C3H7O2• and C4H9• + O2 

16

C4H9O2• systems will be determined with an approximate method. Finally the corresponding rate

17

constant of alkylperoxy radical reactions are presented. Conclusions will be given in Sec. IV.

18

II. COMPUTATIONAL DETAILS

19

The following five mechanisms, M1-M5, for the ethyl + O2 reaction are studied in detail and

20

most notations are consistent with those in an earlier work.31 Meanwhile, similar reactions were

21

studied for propyl + O2 and butyl + O2 systems, and the reactions of the specific studies are listed

22

in Scheme S1 of Supporting Information.

23

M 1 : C2 H 5  O2  TS1  Trans  C2 H 4  HO2


7


1

M 2 : C2 H 5  O2  CH 3CH 2OO  TS 2  CH 2CH 2OOH  TS 5  CH 2CH 2O  OH

2

M 3 : C2 H 5  O2  CH 3CH 2OO  TS 4  CH 3CHO  OH

3

M 4 : C2 H 5  O2  CH 3CH 2OO  TS 2  CH 2CH 2OOH  TS 3  C2 H 4  HO2

4

M 5 : C2 H 5  O2  CH 3CH 2OO  TS1  C2 H 4  HO2

5

1. Electronic Structure Calculations.

Page 8 of 39

6

CCSD(T)/aug-cc-pVTZ is adopted in geometry optimization of the reactants, transition states,

7

intermediates, and products in the reactions between ethyl radical and the oxygen molecule.

8

Molecular geometries were considered converged when the rms gradient fell below 10-6

9

hartree/bohr in these calculations. It is in many cases advantageous to perform geometry

10

optimizations using a pre-calculated force-constant matrix. Such a force-constant matrix is

11

obtained using CCSD/cc-pVTZ. CCSD(T) is too expensive to be applied to determine structures

12

for propyl + O2 and butyl + O2 systems, and DFT calculations are still required. Geometries for

13

ethyl + O2 system obtained from CCSD(T)/aug-cc-pVTZ is employed as the benchmark to

14

evaluate performance of the following six XC functionals: B3LYP,29 BHandHLYP,45 TPSSh,46

15

M06-2X,30 VSXC,47 PBE1PBE48 with 6-311++g(d,p) on structures in these reactions.

16

Differences between bond lengths and bond angles of stationary points with these XC

17

functionals and those with CCSD(T)/aug-cc-pVTZ are adopted as a criteria in choosing a proper

18

XC functional. In addition, the energy of the optimized configuration with these methods can

19

also be used to evaluate performance of these XC functionals. The CCSD(T)-F12a/cc-pVDZ-F12

20

method, which will be shown to be an economical and reliable method in energy calculations, is

21

adopted to calculate the single-point energies at structures obtained in DFT calculations. The

22

structure with the lowest energy is optimal for stable molecules. For transition states, the energy

23

difference between that at the DFT optimized structures and that at the CCSD(T) optimized


8

Page 9 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

structures is adopted as a criterion in evaluating these XC functionals. A proper XC functional

2

will be chosen based both on the structure and single-point energy to deal with propyl + O2 and

3

butyl + O2 systems.

4

In the calculation of reaction energy barriers and reaction energies, CCSD(T) method with the

5

aug-cc-pVQZ basis set is difficult to be applied to propyl + O2 and butyl + O2 systems. On the

6

other hand, DLPNO-CCSD(T)/aug-cc-pVQZ, CCSD(T)-F12a/aug-cc-pVDZ, and CCSD(T)-

7

F12a/cc-pVDZ-F12 methods are economical and possibly reliable methods. Reaction barrier

8

heights and reaction energies for the ethyl + O2 system with DLPNO-CCSD(T)/aug-cc-pVQZ,

9

CCSD(T)-F12a/aug-cc-pVDZ and CCSD(T)-F12a/cc-pVDZ-F12 will be compared with those of

10

CCSD(T)/aug-cc-pVQZ. The most reliable method will be applied to propyl + O2 and butyl + O2

11

systems.

12

Optimized geometries, vibrational frequencies and hindered internal rotor potentials are

13

computed with B3LYP, which will be shown to be among the most reliable functional to

14

determine structures. A scale factor of 0.963 recommended by Kashinski49 was applied to the

15

B3LYP/6-311++g(d,p) frequencies. All the transition-state structures were conﬁrmed with one

16

and only one single imaginary frequency. And intrinsic reaction coordinate50 (IRC) calculations

17

are carried out to verify the transition states structure that connect the desired reactants and

18

products. For low-frequency torsional motions, the one-dimensional (1-D) hindered internal rotor

19

model51 is adopted to handle these internal rotations. The potentials of each internal rotation

20

molecules are calculated at the same level as that for geometry optimization using a relaxed

21

energy scan of the dihedral angle. For transition states, internal rotor scans are performed by

22

freezing the atoms participated in the reaction coordinate.


9


Page 10 of 39

1

For the barrierless alkyl + O2  ROO• reactions, relaxed energy scans are performed along the

2

C–O bond with a 0.1 Å interval using the broken-symmetry B3LYP/6-311++g(d,p) to achieve

3

structures along the minimum energy path (MEP). In the B3LYP calculations, correctness of the

4

dissociation curve is verified by ensuring that the total energy of alkyl + O2 at an infinite distance

5

is equal to the sum of the energies of separate alkyl radical and O2 molecule. Vibrational

6

frequencies for points along the MEP are also obtained at the same level with low-frequency

7

torsional modes treated using the hindered internal rotor model. Energies on the MEP for

8

CH3CH2OO•  C2H5• + O2 reaction are further obtained with the EOM-SF-CCSD(dT)52 and

9

EOM-SF-CCSD methods at geometries obtained with B3LYP/6-311++g(d,p). In these

10

calculations, CCSD(T) is adopted in calculating the energy of the quartet state with MS=3/2

11

along the MEP. EOM-SF-CCSD(dT) or EOM-SF-CCSD with aug-cc-pVTZ is employed to

12

calculate energy difference between the doublet and quartet states. In EOM-SF-CC calculations,

13

the quartet state with MS=3/2 and the quartet state with MS=1/2 are not degenerate due to spin

14

contamination. To achieve reliable total energy for the double state, this energy difference is

15

obtained from the quartet state with Ms=1/2 and the doublet state whose energies are determined

16

at the same level in EOM-SF-CC calculations.

17

The EOM-SF-CC method is still too expensive for propyl + O2 and butyl + O2, and lower level

18

methods are required to obtain energies on the MEP. DFT with broken-symmetry is usually

19

employed to determine the PESs of barrierless reactions, i.e. the same method to determine the

20

structures on the MEP. On the other hand, the PESs can also be represented approximately using

21

the Morse potential with the form: V(R)=De{1-exp[-β(R-Re)]}2, where De is the dissociation

22

energy excluding zero point energy, Re is the equilibrium bond length, β is related to the force

23

constant of the corresponding bond. Reliability of energies along the MEP with broken-


10

Page 11 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

symmetry DFT and those with the Morse potential is evaluated by comparing the relative

2

energies and rate constants for the reaction CH3CH2OO•  C2H5• + O2 with those from EOM-

3

SF-CCSD(dT). In broken-symmetry DFT calculations, energies at the MEP with respect to

4

CH3CH2OO• are scaled by the ratio between reaction energy from CCSD(T)-F12a and that with

5

DFT. On the other hand, only β needs to be determined for the Morse potential, while De is

6

taken to be the dissociation energy with CCSD(T)-F12a, and Re is the optimized bond length

7

with B3LYP. In the present work, β is determined by fitting five points along the MEP with the

8

C-O bond distance starting from the optimized bond length and an interval of 0.05Å using

9

CCSD(T)-F12a/cc-pVDZ-F12.

10

2. Rate constant calculations

11

In this study, rates constants with tight transition state at high pressure limit are calculated with

12

the canonical transition state theory53 (TST), and the variational transition state theory39 (VTST)

13

is employed for barrierless reactions. Pressure-dependent rate constants are calculated with the

14

RRKM/ME theory at pressures varying from 0.1 to 100 atm. Our calculations are carried out

15

using the MESS program.54 The Eckart tunneling corrections55 were included for channels with

16

tight transition states. The collision frequency is estimated by using the Lennard-Jones (L-J)

17

potential, where the following L-J parameters are employed:  =4.87 Å and  =403.37 K for

18

C2H5O2,  =5.21 Å and  =454.77 K for C3H7O2,  =5.52 Å and  =502.26 K for C4H9O2. The

19

L-J parameters for the wells are estimated based on critical temperature and pressure.56 The L-J

20

parameters for the collider N2 with  =3.62 Å and  =97.53 K are taken from the JetSurF2.0

21

database.57 A single-exponential down model is applied to energy transfer per downward

22

collision,  E down , and a value of 200 (T/300)0.85 cm-1 is used to account for the temperature


11


Page 12 of 39

1

dependence of the energy transfer parameter. This approximation has been used in previous

2

related study. 26,33

3

To facilitate application of these reaction rate constants in the practical modeling, the modified

Ea ) is adopted to describe the temperature-dependent rate RT

4

Arrhenius form k (T )  AT n exp(

5

constants at different pressures. All of the rate constants of investigated reactions for the

6

temperature ranging from 500-2500K with an interval of 100K are calculated and fitted to the

7

modified Arrhenius expression.

8

In this work, all the DFT calculations were performed with the G09 program package,58

9

whereas the CFOUR program package59 was used to carry out CCSD(T)/aug-cc-pVTZ

10

optimization for geometries of ethyl system. The DLPNO-CCSD(T) calculations were performed

11

with the ORCA program,60 while the other CCSD(T) and CCSD(T)-F12a energies were obtained

12

with the MOLPRO program.61 The EOM-SF-CC calculations were performed with the Q-Chem

13

program.62 The rate constants are obtained with the MESS software.54 All calculations are

14

carried out on National Supercomputing Center of ShenZhen.

15

III. RESULTS AND DISCUSSION

16

1. Structures of TSs and stable molecules

17

Optimized geometries of the reactants, transition states and products for the ethyl + O2

18

system using the six XC functionals with 6-311++g(d,p), CCSD(T)/ aug-cc-PVTZ and available

19

experimental data are illustrated in Figures S1-10 of Supporting Information. Experimental

20

geometries are available for the following species: C2H5•，O2, OH，C2H4, HO2, CH2CH2O and

21

CH3CHO. The average and maximum absolute deviation on geometric parameters between the

22

experimental values and those with CCSD(T)/aug-cc-pVTZ are listed in Table S1 of Supporting

23

Information. It can be seen that CCSD(T)/aug-cc-pVTZ geometries agree very well with


12

Page 13 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

experimental values, and the maximum absolute deviations of bond lengths and bond angles are

2

0.0076Å and 1.327°, respectively. Furthermore, geometries of CH3CH2OO• and TS1 using

3

CCSD(T)/aug-cc-pVTZ are in good agreement with previous CCSD(T)/cc-pVQZ results,32 with

4

a maximum difference on the bond lengths and bond angles of 0.006 Å and 0.3°. According to

5

the work by Schaefer et. al.31 optimized CH3CH•OOH was located with CCSD/DZP and

6

CCSD(T)/DZP. However, the CH3CH•OOH configuration could not be obtained with

7

CCSD(T)/aug-cc-pVTZ, and aldehydes and OH are formed during optimization for

8

CH3CH•OOH. In fact, a four-membered ring-like TS is required to produce CH3CH•OOH from

9

CH3CH2OO• and such a TS is difficult to be formed. Furthermore, this structure is not reported

10

by Schaefer et. al.32 in their work on geometry optimization of this system at the CCSD(T)/cc-

11

pVQZ level. On the other hand, CH3CH•OOH can be obtained only using BHandHLYP and

12

M06-2X with 6-311++g(d,p). In addition, the configuration CH3CH•OOH is not reported in

13

Refs. 16 and 27 either.

14

To evaluate performance of these XC functionals on geometries of transition states and stable

15

molecules, absolute mean deviation (AMD) and maximum absolute deviation (MAD) on bond

16

lengths and bond angles of each XC functionals with respect to those of CCSD(T) are listed in

17

Table 1. It should be noted that TS1-Trans cannot be located using the M06-2X functional. It can

18

be seen from this table that the bond lengths and bond angles of transition states with PBE1PBE

19

agree best with those of CCSD(T) and BHandHLYP also provides reasonable bond lengths.

20

AMD on bond length with B3LYP is the smallest, but the MAD of bond lengths with B3LYP is

21

somewhat larger. In fact structures with TPSSh are also reasonable, while error of VSXC on TS

22

structures is more pronounced. On the other hand, error in the bond angles of TS1-Trans with

23

these XC functionals is rather larger, but effect of bond angles on total energy is usually


13


Page 14 of 39

1

insignificant. As for stable molecules, B3LYP and TPSSh afford structures that are in good

2

agreement with CCSD(T) results. The other XC functionals also give reliable structures for

3

stable species except that the error of BHandHLYP is somewhat larger.

4

Table 1. The absolute mean deviation (AMD) and maximum absolute deviation (MAD) on bond

5

lengths and bond angles of each XC functionals with respect to those of CCSD(T) for the

6

transition states and the stable molecules in the C2H5• + O2 system. B3LYP

BHandHLYP

TPSSh

M06-2X

VSXC

PBE1PBE

AMD

0.0097

0.0152

0.0119

(0.0107)

0.0199

0.0109

MAD

0.0996

0.0519

0.1148

(0.0632)

0.2522

0.0372

AMD

1.26

1.30

1.32

(0.85)

2.27

1.14

MAD

10.97

9.64

12.51

(3.25)

13.07

9.71

AMD

0.0037

0.0137

0.0033

0.0071

0.0030

0.0082

MAD

0.0092

0.0339

0.0098

0.0305

0.0244

0.0280

AMD

0.47

0.45

0.37

0.32

0.48

0.46

MAD

1.71

2.16

1.43

2.15

1.31

1.95

　 Transition states

d/Å

A / (°)

Stable molecules

d/Å

A / (°) 7 8

In order to re-confirm the optimal XC functional for geometries of the transition states and the

9

stable molecules, energies with CCSD(T)-F12a/cc-pVDZ-F12 at structures obtained from DFT

10

calculations are compared with those at the CCSD(T)/aug-cc-pVTZ structures. CCSD(T)-

11

F12a/cc-pVDZ-F12 will be shown later to provide reaction energies and barrier heights that

12

agree best with CCSD(T)/aug-cc-pVQZ. Stable molecules are local minimum on PES and the

13

optimal XC functional should provide structures of stable species with the lowest energies. We


14

Page 15 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

note that CCSD(T)-F12a/cc-pVDZ-F12 energies at CCSD(T)/aug-cc-pVTZ structures are the

2

lowest on average for stable species. Mean deviation in energy difference between DFT

3

structures and CCSD(T) structures of the stable molecules for the six XC functionals is listed in

4

Table 2. It can be seen from this table that the average total energy with B3LYP is the lowest,

5

which is consistent with the fact that B3LYP structures agree best with CCSD(T) structures for

6

stable molecules. Things are more complicated for TSs since they are saddle points on PES. The

7

absolute mean deviation in energy difference between DFT structures and CCSD(T) structures is

8

employed to choose the optimal XC functional for TS structures and they are listed in Table 2.

9

One can see from this table that single point energies at B3LYP structures or TPSSh structures

10

agree best with those at the CCSD(T) structures.

11

Table 2. Energy difference using CCSD(T)-F12a/cc-pVDZ-F12 at DFT structures and those at

12

CCSD(T) structures for stable molecules and transition states for the C2H5• + O2 system. (unit:

13

kcal/mol) B3LYP

BHandHLYP TPSSh M06-2X VSXC

PBE1PBE

Stable moleculesa

0.026

0.301

0.030

0.125

0.066

0.103

Transition statesb

0.167

0.373

0.165

0.272

0.525

0.209

14

15

aMean

16

bAbsolute

deviation in energy difference between DFT structures and CCSD(T) structures; mean deviation in energy difference between DFT structures and CCSD(T) structures.

17

Based on the above results, we conclude that the B3LYP functional with the 6-311++g(d,p)

18

basis set is superior for geometries in the C2H5• + O2 system from the aspect of bond lengths,

19

bond angles and energies of both transition states and stable molecules. Meanwhile, our


15


Page 16 of 39

1

conclusion are consistent with Ignatyev's results27 that B3LYP functional is the optimal

2

functional on structures of this system. We believe this functional is also promising for the

3

similar C3H7• + O2 and C4H9• + O2 systems. All the stationary points in the C3H7• + O2 and

4

C4H9• + O2 systems are optimized using B3LYP/6-311++g(d,p) and their Cartesian coordinates

5

are given in the Supporting Information.

6

2. Reaction Energies and Barrier Heights

7

To achieve reaction energies and barrier heights with high accuracy for these systems, the

8

explicitly correlated CCSD(T)-F12a or the locally correlated DLPNO-CCSD(T) approach will be

9

adopted. Accuracy of DLPNO-CCSD(T)/aug-cc-pVQZ, CCSD(T)-F12a/aug-cc-pVDZ, as well

10

as CCSD(T)-F12a/cc-pVDZ-F12 on the energy differences will be evaluated by comparing their

11

results with those of CCSD(T)/aug-cc-pVQZ for the C2H5• + O2 system. Energies for transition

12

states, intermediates, and products in the C2H5• + O2 system with respect to energy of the

13

reactant C2H5• + O2 using these methods are given in Table S2 of Supporting Information.

14

Differences between forward and reverse barrier heights as well as reaction energies in the C2H5•

15

+ O2 system with these methods and those with CCSD(T)/aug-cc-pVQZ are listed in Table 3.

16

Our results show that reaction energies with DLPNO-CCSD(T) agree rather well with those of

17

CCSD(T), while their difference in barrier height is somewhat larger. Barrier heights for M1:

18

C2H5• + O2  Ts1-Trans  C2H4 + HO2 with DLPNO-CCSD(T)/aug-cc-pVQZ are overestimated

19

by about 3.8 kcal/mol compared with those of CCSD(T)/aug-cc-pVQZ. Agreement between the

20

energy differences with CCSD(T)-F12a/aug-cc-pVDZ and those of CCSD(T)/aug-cc-pVQZ is

21

always better than DLPNO-CCSD(T), but their difference in barrier heights for M1 is still larger

22

than 2.5 kcal/mol. Reaction energies and barrier heights with CCSD(T)-F12a/cc-pVDZ-F12

23

agree best with those of CCSD(T)/aug-cc-pVQZ and their difference is less than 0.4 kcal/mol


16

Page 17 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

except for the barrier heights of M1. Based on these results, we conclude that energy differences

2

for the C2H5• + O2 system with CCSD(T)-F12a/cc-pVDZ-F12 is the best among these three

3

methods.

4

Table 3. Differences in energy barriers and reaction energies of the C2H5• + O2 system from

5

those of CCSD(T)/aug-cc-pVQZ using different methods. (unit: kcal/mol) DLPNO-CCSD(T) CCSD(T)-F12a /aug-cc-pVQZ /aug-cc-pVDZ

CCSD(T)-F12a /cc-pVDZ-F12

CH3CH2OO•  TS1  C2H4+HO2  E1

0.79

0.10

0.14

 E-1

1.24

-0.91

0.09

H

-0.45

1.02

0.05

 E1

3.76

-2.74

-1.39

 E-1

3.78

-2.51

-1.28

H

-0.02

-0.23

-0.11

 E1

0.44

-0.11

0.07

 E-1

0.48

-0.27

-0.16

H

-0.04

0.16

0.22

 E1

0.60

0.17

-0.13

 E-1

1.01

-0.69

0.04

H

-0.41

0.86

-0.17

 E1

0.40

0.37

0.18

 E-1

1.10

-0.16

0.37

C2H5•+O2  TS1-Trans  C2H4+HO2

CH3CH2OO•  TS2  •CH2CH2OOH

•CH2CH2OOH  TS3  C2H4+HO2

CH3CH2OO•  TS4  CH3CHO+OH


17


H

Page 18 of 39

-0.70

0.53

-0.19

 E1

0.28

0.27

0.05

 E-1

1.02

-0.02

0.23

H

-0.75

0.29

-0.18

•CH2CH2OOH  TS5  CH2CH2O+OH

1 2

 E1 is the difference of forward energy barriers between these methods and CCSD(T)/ aug-cc-

3

pVQZ;  E-1 is the difference of backward energy barriers between these methods and CCSD(T)/

4

aug-cc-pVQZ;  H is the difference of reaction energies between these methods and CCSD(T)/

5

aug-cc-pVQZ.

6

To further evaluate accuracy of reaction energies and barrier heights with CCSD(T)-F12a/cc-

7

pVDZ-F12 for these systems, relative energies of stationary points in the C2H5• + O2 system

8

obtained in this work together with those using ANL016 and FPA 32 methods are listed in Table 4

9

when available. In ANL0, CCSD(T) results at CBS is obtained using aug-cc-pVQZ and aug-cc-

10

pV5Z basis set. Besides basis set extrapolation, contribution of higher order excitations to

11

correlation energy is considered in FPA. It can be seen from Table 4 that difference between the

12

present results and those with ANL0 or FPA is generally less than 0.4 kcal/mol with a maximum

13

difference of about 0.7 kcal/mol. In fact, our results usually lie between those of ANL0 and FPA

14

and the largest difference between results of ANL0 and FPA reaches 0.86 kal/mol. This again

15

shows reliability of CCSD(T)-F12a/cc-pVDZ-F12 on energy differences of this systems and

16

reaction energies and barrier heights for all the systems studied in this work are obtained using

17

this method. A schematic plot of PES for the C2H5• + O2 system at the CCSD(T)-F12a/cc-

18

pVDZ-F12 level is provided in Figure 1. Schematic plots for PESs of C3H7• + O2 and C4H9• + O2

19

systems are provided in Figure S11-14 of Supporting Information. We note that relative energies


18

Page 19 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

in the C3H7• + O2 systems also agree well with those by Goldsmith et al.26 with a maximum

2

difference of 1 kcal/mol.

3

Table 4. Relative energies of stationary point in the C2H5+O2 system.a

4 5

a

Stationary point

Present work

ANL0

FPA

Energy16

Energy32

CH3CH2• + O2

0

0

0

CH3CH2OO•

-32.86

-32.77

-32.99

•CH2CH2OOH

-15.98

-15.74

CH3CH2OO•  •CH2CH2OOH

3.86

4.10

CH3CH2OO•  •CH3CHO + OH

9.08

9.32

•CH2CH2OOH  CH2CH2 + HO2

0.48

-0.13

•CH2CH2OOH  CH2CH2O + OH

-1.78

-2.06

CH2CH2 + HO2

-13.02

-13.38

-13.74

CH3CHO + OH

-60.00

-59.33

-60.10

CH2CH2O + OH

-32.44

-31.87

-32.73

All energies include vibrational zero-point energy(ZPE) corrections and are in kcal/mol relative

to C2H5•+O2.

6


19


Page 20 of 39

1

Figure 1. A schematic plot of the C2H5•+O2 PES calculated at the CCSD(T)-F12a/cc-pVDZ-F12

2

level of theory at 0 K, including zero-point energy correction from B3LYP/6-311++G(d,p). (unit:

3

kcal/mol)

4

According to our results, the C-O bond energy in n-C3H7O2 (33.21kcal/mol) or n-C4H9O2

5

(32.93kcal/mol) is slightly smaller than that in i-C3H7O2 (34.76kcal/mol) or s-C4H9O2

6

(34.99kcal/mol). This trend does not agree with relative stability of alky radicals: secondary alky

7

radical is more stable than primary radical. The energy of ROO• is lower than that of R• + O2,

8

and one would expect a smaller C-O bond energy for a more stable R radical. In fact the bond

9

energy of R-X is also closely related to electronegativity of X and it has been found that trend in

10

R-X bond energy will be reversed to relative stability of alky radicals for X with large

11

electronegativity,63 which is consistent with the present results. In fact, similar trend on this bond

12

energy has also been observed in Ref.19. In addition, difference in the C-O bond energies from

13

primary alkyls or from secondary alkyls is less than 0.35 kcal/mol regardless of the size of the

14

alkyl radical.

15

Isomerization reactions ROO•  •QOOH via a ring-like TS play an important role in the alkyl

16

+ O2 systems since the hydroperoxy alkyl radicals are believed to be the key to chain

17

branching.64 The forward and reverse barrier heights of the RO2 isomerization reactions

18

investigated in this work are listed in Table 5. It can be seen from this table that the size of the

19

ring in the transition state structure has a pronounced effect on both the forward and the reverse

20

barrier heights. The carbon site of the H-transfer also affects the barrier heights to some extent.

21

Barrier heights for reactions with a 6-member-ring TS are similar to those with a 7-member-ring

22

TS and they are about 10 kcal/mol smaller than reactions with a 5-member-ring TS for H-

23

transfer from the same type of carbon site. In addition, forward barrier heights for H-transfer


20

Page 21 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

from a primary site are about 3 kcal/mol higher than those from a secondary site with the same

2

ring size of TS. Similar conclusions have been reached in previous works.14,20 On the other hand,

3

the reverse barrier heights do not depend on the carbon site of H-transfer. Reaction energies also

4

depend on the carbon site of H-transfer and they are about 3 kcal/mol larger for H-transfer from

5

a primary site than that from a secondary site. This is consistent with relative stability of alky

6

radicals. One can also see from Table 5 that difference in the forward and reverse barrier heights

7

as well as reaction energies are around 1 kcal/mol for isomerization reactions of different

8

alkylperoxy radicals with the same carbon site for H-transfer and the same ring size of transition

9

states.

10

Table 5. Forward and reverse barrier heights and the reaction energies (include vibrational zero-

11

point energy (ZPE) corrections) for RO2 isomerization reactions. E1a(kcal/mol)

E-1b(kcal/mol)

ΔHc(kcal/mol)

CCOO•  •CCOOH

36.72

19.85

16.87

CC(C)OO•  •CC(C)OOH

36.15

18.71

17.44

CCOO•CC  •CCOOHCC

36.08

18.64

17.44

CCCOO•  C•CCOOH

33.26

19.64

13.62

CCCCOO•  CC•CCOOH

33.00

19.19

13.81

CCOO•CC  CCOOHC•C

32.96

18.82

14.14

CCCOO•  •CCCOOH

25.01

9.13

15.88

CCOO•CC  CCOOHCC•

24.83

8.44

16.38

CCCCOO•  C•CCCOOH

22.42

9.39

13.02

Reaction 1,4-H(p)

1,4-H(s)

1,5-H(p)

1,5-H(s)


21


Page 22 of 39

1,6-H(p) CCCCOO•  •CCCCOOH

25.30

9.53

15.77

1 2

aE 1

is the forward energy barrier; bE−1 is the reverse energy barrier; cΔH is the reaction energy.

3

3. Reaction Rate Constants

4

The dissociation reaction of alkylperoxy radicals to form alkyl + O2 is a barrierless reaction.

5

The alkylperoxy radicals are lower in energy than the dissociation product alkyl + O2, but the

6

dissociation reaction becomes important at elevated temperatures. The VTST method is

7

employed to calculate the rate constant for this reaction. The PES are obtained using EOM-SF-

8

CCSD, EOM-SF-CCSD(dT), broken-symmetry DFT and the Morse potential for the

9

CH3CH2OO•  C2H5• + O2 reaction. PES with EOM-SF-CCSD(dT) is the most accurate one

10

and our results show that difference in the PES with these approaches is less than 1.6 kcal/mol

11

for the CH3CH2OO•  C2H5• + O2 reaction at the points employed in VTST calculations.

12

However, reaction rate constants with these different PESs could be pronounced particularly at

13

low temperatures. Our results shows that the rate constant with PES from broken-symmetry DFT

14

is 7.6 times that from EOM-SF-CCSD(dT) at 300K and it is 2.88 times at 500K. Difference

15

between rate constants with these PESs decreases as temperature rises. This indicates that rate

16

constants obtained in this work at low temperatures are unreliable and they are reasonable for

17

temperatures higher than 500K. On the other hand, the ratio between rate constants from the PES

18

with EOM-SF-CCSD and those with EOM-SF-CCSD(dT) is 1.34 at 500K and this ratio is 1.12

19

for the PES with the Morse potential. These results show that the rate constants with different

20

PESs agree reasonably with each other. For dissociation of C3H7O2 and C4H9O2, their PESs are

21

obtained with the Morse potential in this work.


22

Page 23 of 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

The complete set of rate constants with modified Arrhenius parameters at all studied T/P for

2

the ethyl + O2, n-propyl + O2, isopropyl + O2, n-butyl + O2, and sec-butyl + O2 reactions are

3

provided in the Table S4 of Supporting Information. The rate constants obtained in this work

4

would be valuable to improve the accuracy of chemical kinetic models for low temperature

5

oxidation of large straight-chain alkanes which are important components in kerosene and

6

gasoline. It should be noted that there are many other reactions that have sizeable effects on low

7

temperature combustion of alkanes. Rate constants of those reactions are generally not the same

8

in different mechanisms. In additional, there still exist some uncertainties in the rate constants

9

obtained in this work. These rate constants need to be adjusted to some extent if they are to be

10

employed in a specific reaction mechanism.

11

The rate contants at the high-pressure limit rate for C2H5O2 dissociate reaction are shown in

12

Figure 2(a), together with those from Guo et al.,65 Sheng et al.,28 Villano et al.,20 and

13

Klippenstein et al.33 It can be seen from this figure that rate constants in the present work agree

14

reasonably well with those in previous works and the present rate constants are somewhat

15

smaller for T

A Theoretical Study on Reactions of Alkylperoxy Radicals

Recommend Documents