Gas-Phase Ozone Reactions with a Structurally Diverse Set of

Article Cite This: J. Phys. Chem. A XXXX, XXX, XXX−XXX

pubs.acs.org/JPCA

Gas-Phase Ozone Reactions with a Structurally Diverse Set of Molecules: Barrier Heights and Reaction Energies Evaluated by Coupled Cluster and Density Functional Theory Calculations Daniela Trogolo,† J. Samuel Arey,†,‡,§ and Peter R. Tentscher*,‡,∇ School of Architecture, Civil and Environmental Engineering, É cole Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland ‡ Eawag, Swiss Federal Institute of Aquatic Science and Technology, Ü berlandstrasse 133, 8600 Dübendorf, Switzerland Downloaded via IOWA STATE UNIV on January 4, 2019 at 13:27:38 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

S Supporting Information *

ABSTRACT: Reactions with ozone transform organic and inorganic molecules in water treatment systems as well as in atmospheric chemistry, either in the aqueous phase, at gas/ particle interfaces, or in the gas phase. Computed thermokinetic data can be used to estimate the reactivities of molecules toward ozone in cases where no experimental data are available. Although the gas-phase reactivity of olefins with ozone has been characterized extensively in the literature, this is not the case for the richer chemistry of ozone with polar molecules, which occurs in the aqueous phase or in microhydrated environments. Here, we selected a number of model reactions with small molecules (ethene, ethyne, hydrogen cyanide, hydrogen chloride, ammonia, bromide, and trimethylamine) to study the accuracy of different quantum chemical methods for describing the reactivities of these molecules with ozone. We calculated benchmark electronic energies of gas-phase reactions of these systems with singlereference coupled cluster (CC) theory. These benchmark results for the binding energy in the van der Waals complex, the energy of the transition structure, and the reaction energy were estimated to be accurate within 1−2 kcal mol−1. Singlet oxygen (1O2) is a common product of ozone reactions. Coupled cluster calculations with up to perturbative quadruples (CCSDT(Q)) were needed to obtain reaction energies accurate within 1 kcal mol−1 when this species was involved. In (micro)hydrated environments or at interfaces, coupled cluster methods are prohibitively expensive in most cases. We tested the suitability of some contemporary density functional theory (DFT) methods to reproduce the benchmark electronic energy differences. Range-separated functionals were found to be promising candidates to estimate forward barrier heights, with LC-ωPBE rivaling the accuracy of CCSD(T). For energies of reaction, however, DFT methods exhibited large systematic errors, depending on their fraction of orbital exchange. This was found to worsen when 1O2 is a product, and no safe recommendation can be given for DFT reaction energies in such cases.

1. INTRODUCTION

In both aqueous and atmospheric reactions, it is desirable to estimate second-order rate constants of the initial reaction of ozone with the organic molecule, as well as to elucidate the reaction mechanisms that lead to the observed products of ozonolysis. Knowledge of (estimated) second-order rate constants can be used to predict the atmospheric half-lives of organics6 or the efficiency of their removal in water treatment systems.7 In water treatment, the transformation products of such reactions are of toxicological relevance.8 To predict potential reaction routes, mechanistic insight is essential.9,10 Computational quantum chemistry can be used to investigate these reactions mechanisms, to estimate rate

1.1. Reactivity of Ozone. Ozone plays a crucial role as an oxidizing agent in atmospheric and aquatic chemistry. In the atmosphere, ozonolysis is involved in the formation of secondary organic aerosols from volatile organic compounds, a process with broad implications for the residence of organics in the atmosphere as well as for cloud formation.1 In water treatment, ozone is a widely applied oxidation agent for the disinfection of water, the removal of taste and odor compounds, and the abatement of organic micropollutants.2−4 In these aqueous-phase reactions, in which we are chiefly interested, a diverse set of functional groups can be expected to react with ozone. These include electron-rich aromatics, aliphatic amines, olefins, thiols, or thioethers, as well as organic and inorganic ions.5 © XXXX American Chemical Society

Received: October 23, 2018 Revised: December 6, 2018

A

DOI: 10.1021/acs.jpca.8b10323 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

renormalized CC theory, which is more appropriate for diradicals. They showed that, when including quadruples excitations in the CC treatment, the results obtained by conventional and completely renormalized CC theories are very similar. CCSD(T) is considered the “gold standard” for systems without multireference character, reaching chemical accuracy (1.0 kcal mol−1) for most problems. The authors report that RHF-based CCSD(T) results with an augmented triple-ζ basis differ from their best estimates of activation energy and reaction energy by 0.5 kcal mol−1 to 1.5 kcal mol−1. These studies show that closed-shell, single-reference CC theory can provide accurate thermokinetic data for O3 reactivity, at least toward these olefins. Ultimately, larger molecules require the use of protocols that are computationally less demanding than CCSD(T). DFT methods are the obvious choice from the standard toolbox of computational quantum chemistry. Zhao et al.25 tested a wide variety of density functionals for the thermokinetic energies of the reaction of O3 with C2H2 and C2H4. They reported energy differences between reactants and van der Waals complexes, transition structures, and products. For these systems, (comprehensive) mean unsigned errors (MUEs), with respect to the best estimates, were reported. The best performing DFT method was M05 (with a quadruple-ζ basis) exhibiting a MUE of 1.1 kcal mol−1, compared to 0.3 kcal mol−1 for CCSD(T) with a triple-ζ basis set. The error in M05 results stems chiefly from the reaction energies, and excluding those data lowers the MUE to 0.4 kcal mol−1. For comparison, the B3LYP functional exhibited MUEs of 2.7 and 3.3 kcal mol−1 with and without reaction energies, respectively. According to these data, thermokinetic calculations with the M05 functional should yield reliable (within ∼1 kcal mol−1 of the best estimates) results for the reactions between O3 and olefins. The reactivity of O3 is not limited to addition reactions to olefins, but high-level quantum chemical studies on other systems are rare. Additional gas-phase theoretical work included a hydrogen abstraction mechanism from NH3,26 H2S,27 alcohols,27 and even saturated hydrocarbons in polymers.28 Generally, these reactions seem to exhibit multiple possible reaction channels, several of which may be contributing to the observed reactivity. The employed model chemistries (closed-shell, single-reference CC energies with geometries obtained at a lower computational level) show that these systems have one favorable reaction channel in common: the insertion of O3, forming a hydrotrioxide (−OOOH) either through H-abstraction followed by recombination, or through a concerted mechanism. In one example, the H-abstraction reaction of O3 from methanol,29 post-CCSD(T) effects were considered. This led to a change in the computed barrier height of >1.2 kcal mol−1, which is slightly higher than the effect of the post-CCSD(T) treatment on the barrier heights for the Criegee addition to C2H2 and C2H4.25 O3 can also react with halide anions in both gas-phase and hydrated environments. Theoretical work has been performed for reactions with Br−30,31 and I−.32 In this case, no hydrogen can be abstracted, and the computed mechanism involves the addition of O3 to the halide, forming a trioxide. This initial step could be followed by the cleavage of O2. The gas-phase study by Gladich et al.31 on O3 + Br− employed exclusively unrestricted WF methods (U-MP2, U-CCSD, U-CCSD(T)), whereas Rayne and Forest30 used the composite G4 and G4MP2 protocols, presumably with a restricted reference wave

constants from computed barrier heights, and to evaluate the thermodynamic feasibility of the reactions of O3. Potential applications include the prediction of O3 regioselectivity, crucial for the prediction of toxic transformation products formed from the reactions of organic micropollutants with several O3-reactive sites. However, all applications necessitate the calculation of reliable barrier heights and reaction free energies. Therefore, appropriate quantum chemical methods for these calculations need to be identified. 1.2. Ozone as a Problematic Case in Theoretical Chemistry. The ozone molecule exhibits singlet diradical character,11 and it is considered to be a multireference system. Depending on differing definitions of diradical character, O3 has been found to bear 16%,12 18%,13 or even 44%.11,14 For such systems, the use of common, single-reference wave function (WF) or density functional theory (DFT) methods can be problematic.15−18 These problems have been studied in some detail for gas-phase reactions, where theoretical work has been mostly concerned with the addition of the ozone molecule to the double bond of unsaturated systems, discussed below. Other ozone-reactive functional groups render a molecule more polar, and they are, thus, rarely found in volatile organic compounds.6,19 High-level quantum chemical studies have focused on the (gas-phase) reaction of O3 with the prototypical olefins ethene (C2H4) and ethyne (C2H2). Not only are these molecules small and, thus, traceable with advanced wave function methods, but also experimental data for the reaction energies and barrier heights are available. These experimental thermodynamic and kinetic properties can be used as benchmarks to which computed data can be compared. For the addition of O3 to C2H4, two possible reaction channels exist on the CASPT220 potential energy surface (PES):21 a concerted cycloaddition (Criegee mechanism) and a stepwise addition with a biradical intermediate (DeMoore mechanism).21 Further single-point calculations with different variants of multireference perturbation theory, but also with the single reference “gold standard” CCSD(T), showed that the activation energy for the Criegee mechanism is much lower than that of the DeMoore mechanism. The author’s best estimates of the barrier heights are 56 kJ mol−1 (DeMoore) versus 12 kJ mol−1 (Criegee), which would make the Criegee mechanism the only relevant channel. Another approach to the diradical problem is to perform calculations with a single reference WF, but using an unrestricted, spin-contaminated UHF reference. The energies (and gradients) can then be corrected by approximate spin projection.22 On such a spin-projected PES, the DeMoore intermediate is depicted as an artifact of spin contamination.17 The predominance of the Criegee mechanism is slightly less clear for dienes: a study using single-reference DFT and WF methods as well as multireference-MP2 showed that the DeMoore mechanism could account for up to 30% (B2PLYP)23 of the total rate constant in trans-butadiene, whereas it reaches only 4% using multireference-MP2 energies.24 Using a similar methodology, these authors estimated that the contribution of the DeMoore mechanism to the overall rate constant of C2H2 is 1%−8%.18 Two studies by the groups of Truhlar25 and Houk16 aimed at providing highly accurate electronic energies for the barrier heights and the reaction energies of the Criegee addition of O3 to C2H2 and C2H4. Their best estimates are based on closedshell, single-reference CC theory, as well as completely B


Article


Figure 1. Reaction pathways considered. Displayed structures correspond to stationary structures optimized on the CCSD/(A-/M-)PVTZ PES. Reaction pathways were initially explored with the CCSD/PVDZ method and DFT methods (ωB97XD or M11).

function, supplemented with a continuum solvation model. A study on the reaction of O3 with I− also used restricted MP2 and CCSD(T) energies.32 Besides the above-mentioned studies that used high-level methods on small molecules, much computational work has

been performed with more approximate model chemistries. In the field of gas phase reactions, one common application is the estimation of absolute reaction rate constants, using either canonical or variational transition-state theory,33,34 or employing the Rice−Ramsperger−Kassel−Marcus (RRKM) model.35 C


Article


methods against a more varied set of thermokinetic data of O3 reactions, making it possible to recommend approximate quantum chemical methods for the study of ozone chemistry.

These studies considered almost exclusively the addition of O3 to olefins, using a wide variety of model chemistries. Often, the DFT geometries are combined with CCSD(T) or MP2 energies. In other cases, DFT methods have also been used to obtain the electronic energies: restricted Kohn−Sham methods included B3LYP36−40 (by far, the most used method), M11,28 MPW1PW91,37 MPW1K,37,39,41 BB1K,42 M06-2X,39 and BH&HLYP.40 Less frequently, unrestricted, spin-contaminated DFT methods were used (for example, UB3LYP43 or UM1128). Only a few gas-phase studies were concerned with molecules other than olefins. Some examples are NH3,26 the ether/thioether thioxane,39 aldehydes,44 benzene,45 catechol,46 and dioxines.47 Previous theoretical studies of aqueous phase (or microhydrated) ozone reactions covered a more diverse set of functional groups and ions: from inorganic ions (halide,48 nitrite,49 HOO−50) over small radicals (HOO•51 and others52), H2S53 and organic sulfides,54 aliphatic secondary amines,55 to larger organic molecules, as in aromatics56 including anilines57 and phenols,58 to complex organics36 as large as sulfamethoxazole.59,60 The majority of these studies36,48−50,53,56−58,60 employed closed-shell B3LYP, but also other closed- and open-shell methods (BH&HLYP,51 B97K,51 spin-projected UB3LYP,55 M06-2X,58 MPW1K,58 and UM0559) have been used. In those studies relying on electronic energies calculated with DFT methods, it appears that rather little attention has been paid to the selection and the accuracy of the chosen DFT methods. For comparison, the omnipresent B3LYP functional yielded errors (with respect to best estimates) of reaction energies and barrier heights of up to 4 kcal mol−1 for the Criegee additions of O3 to C2H2 and C2H4. The error in the B3LYP reaction energy (with respect to best estimates) for C2H2 and C2H4 was not systematic, amounting to ∼4 kcal mol−1 and 0.5 kcal mol−1, respectively.25 Therefore, error cancellation when comparing similar reactions cannot be counted upon. Even less is known about the appropriateness of different density functionals when other reactive moieties are involved, which is relevant to studies on aqueous ozone chemistry. One example is the reaction of O3 with HOO•,61 where DFT methods were compared to CASPT2 barriers for hydrogen and oxygen transfer. None of the few tested DFT methods were able to provide accurate results, with barrier heights computed by hybrid functionals differing by 3 kcal mol−1 from CASPT2 best estimates. Based on fitting of DFT imaginary frequencies to CASPT2 frequencies by varying the fraction of orbital exchange, the authors argued that methods with ∼40% of orbital exchange should be suitable for the problems studied. Nevertheless, these various DFT methods were used to derive thermokinetic data for aqueous-phase reactions (vide supra) but were lacking an estimation of the uncertainty in computed energies arising from the electronic structure method alone (that is, without the additional error contributed by implicitly or explicitly modeled solvation energies). In the present study, we extended the set of highly accurate benchmark thermokinetic data for gas phase ozone reactions to mechanisms other than Criegee additions. For a set of small model molecules (ethene, ethyne, HCN, HCl, NH3, Br−, trimethylamine, see Figure 1). we selected one possible reaction channel, for which we computed accurate thermochemical data based on single-reference CC methods up to CCSDT(Q). This allowed us to test many contemporary DFT

2. COMPUTATIONAL DETAILS 2.1. Coupled Cluster Calculations. 2.1.1. Basis Sets. Dunning’s correlation consistent basis sets of double-, triple-, quadruple-, pentuple-, hextuple-(X = 2, 3, 4, 5, 6)ζ quality were used throughout.62,62−70 X is the maximum angular momentum quantum number represented in the correlation-consistent basis set (for example, X = 4 for aug-cc-pVQZ). We abbreviated as PVXZ for the cc-pVXZ62−65 bases, APVXZ for the augmented aug-cc-pVXZ62,66 bases, PCVXZ for the core−valence cc-pCVXZ67−69 bases, and MAPVXZ for the minimally augmented maug-cc-pVXZ70 bases. 2.1.2. Geometry Optimization. All geometry optimizations were performed at the frozen-core (FC) CCSD71 level of theory, using a restricted Hartree−Fock (RHF) reference, with the Gaussian09 vD.0172 electronic structure suite, using standard settings. It was previously shown that CCSD yields reasonable structures for closed-shell73 and open-shell molecules,74 including multireference systems. Although it is a single-determinantal method, the inclusion of single excitations in CCSD allows for some inclusion of multireference effects in the coupled cluster wave function.73 Basis sets were chosen depending on the requirements of the system (adding diffuse functions in cases where this led to significant changes in the transition structure) and the computational cost: PVTZ for ethyne (C2H2), ethene (C2H4), and hydrogen cyanide (HCN); MAPVTZ for ammonia (NH3) and hydrochloric acid (HCl); APVTZ for bromide (Br−). For trimethylamine (TMA, (CH3)3N), we used the PVTZ basis on O and N atoms, and PVDZ on C and H. As C and H atoms are not involved in the studied TMA reaction channel, errors in computed energies arising from the use of PVDZ bases for the geometries were expected to cancel to a large extent. All stationary structures were characterized by harmonic frequency calculations. The connectivity of stationary points was determined by a relaxed scan of the bond lengths related to the reaction coordinate, performed at the CCSD/PVDZ level of theory. This information was confirmed by intrinsic reaction coordinate (IRC)75 calculations at the DFT level (ωB97XD76 or M1177 ), starting from correspondingly reoptimized transition structures. All species arising in a reaction sequence were optimized with consistent basis sets. O3 was optimized with PVTZ, MAPVTZ, and APVTZ basis sets, so that consistent sets of structures could be used for subsequent electronic energy calculations. 1O2 was treated in the same way. 2.1.3. Electronic Energy Calculations. Accurate thermochemical protocols based on single-reference CC methods need to achieve convergence in two domains: the description of the single electrons (basis set convergence), and the descriptions of the many-body problem (electron correlation).73 Common standardized approaches, such as W478 and HEAT,79 exploit the fact that higher-order (post-CCSD) correlation effects exhibit increasingly rapid basis-set convergence with each higher CC excitation level,80 and these contributions may be added to the total electronic energy as additive corrections with increasingly smaller basis sets. Single-reference CC methods form a hierarchy,73,81,82 starting from CC with single and double excitations from the reference determinant (CCSD), to perturbative (CCSD(T)) D


Article


in the correlation-consistent basis set (e.g., 3 for APVTZ and MAPVTZ, 4 for APVQZ and MAPVQZ). For most CBS extrapolations, we used the APVTZ/APVQZ basis sets; details are given in Table 1 and Section S2 in the Supporting Information (SI). The different correlation energy components were evaluated with different basis sets whenever this was affordable. Core electrons were kept frozen. We note that the default frozen core option for bromine differs between Gaussian09 and CFOUR v1; the CFOUR v1 scheme (freezing the 3d-orbitals) was used throughout. Otherwise, for the RHF-based CC calculations, these software packages produce identical results within the numerical convergence limits of the energy calculation. Post-CCSD(T) calculations with the MAPVDZ and APVDZ basis sets used the aug-cc-pVDZ-PP89 basis with Stuttgart−Koeln pseudopotentials on bromine, and MAPVDZ/APVDZ on all other elements. For all reactions except that of TMA, the energy contribution of core−valence correlation was assessed as an additive term, using the PCVTZ basis:

or full (CCSDT) inclusion of triples excitations, quadruples (CCSDT(Q)/CCSDTQ), pentuples (CCSDTQ(P)/ CCSDTQP), and so on,83 approaching the full configuration interaction (CI) limit that represents the exact wave function. In the present study, we have also included the CC3/CC4/ CC5 methods, 83,84 nonperturbative approximations to CCSDT/CCSDTQ/CCSDTQP, which are usually used for the computation of optical spectra. We included these methods because they should be less sensitive to shortcomings in the reference wave function compared to perturbative approaches. All “best estimates” of electronic energies E and energy differences ΔE are based on calculations that used a closedshell, RHF reference wave function. Additional UHF-based calculations are indicated with a “U”-prefix, whereas the absence of a prefix always refers to RHF-based calculations. The calculations on the triplet ground state of oxygen, O(3P), used a UHF reference, which was not spin-contaminated. All CC calculations were performed with the Gaussian09 or the CFOUR v185 software packages. Post-CCSD(T) calculations in this hierarchy (CC3, CCSDT, CCSDT(Q), CC4, CCSDTQ, CCSDTQ(P), CC5, CCSDTQP) were performed with CFOUR v1 or the MRCC86 code interfaced to CFOUR v1. The manuscript uses three types of energy differences, denoted by different symbols: ΔE denotes an energy difference between different species, δ̃E denotes the difference of the electronic energy of a single species due to an additive correction of correlation treatment (eq 1), and δE denotes the energy difference in ΔE due to an additive correction of correlation treatment (eq 7). In the present protocol, the total electronic energy of a species is given by an additivity scheme:

core‐valence δ Ẽ CCSD( T ) = EAE − CCSD(T )/PCVTZ − E FC − CCSD(T )/PCVTZ

(6)

where AE indicates all electrons, meaning that all orbitals have been included in the correlation treatment, in contrast with the frozen core (FC) approximation. We use ΔE to denote energy differences of various structures, with respect to isolated reactants, calculated with a given method, with subscripts (vdW, TS, rxn, 1O2) indicating the species shown in Figure 1. In analogy to the energy 2 level 2 difference δ̃Elevel level 1, we use δElevel 1 to denote the effect of correlation treatment on an energy difference among two systems. For example,

CCSD CCSD(T ) CC3 Eelec = E HF + δ Ẽ HF + δ Ẽ CCSD + δ Ẽ CCSD( T) CCSDT CCSDT(Q ) core−valence + δ Ẽ CC3 + δ Ẽ CCSDT + δ Ẽ CCSD( T)

(1)

CCSDT δECC3 /APVDZ = ΔE TS(CCSDT/APVDZ)

In this scheme, δ̃E corresponds to the difference in correlation energy, compared to the next lower level of correlation. For example, CCSD(T ) = ECCSD(T ) − ECCSD δ Ẽ CCSD

− ΔE TS(CC3/APVDZ) CCSDT − δ Ẽ C2H2 CCSDT = δ Ẽ transition structureCC3 CC3

− δ Ẽ O3 CCSDT

Subscript indications of the type of energy difference (vdW, TS, rxn, 1O2) are usually dropped from the δE notation. 2.1.4. Multireference Diagnostics and Open-Shell Calculations. For all structures, we computed two multireference diagnostics using both an RHF reference and a UHF reference, using Gaussian09: the T1 diagnostic90 and the %TAE[(T)] diagnostic,91 which is defined as

where the same basis set is used for both energy calculations, unless complete basis set extrapolation (CBS) is indicated. In 88 most cases, we used CBS extrapolation for EHF,87 δ̃ECCSD HF , 88 CCSD(T) and δ̃ECCSD : E HF /CBS = E HF /X +

E HF /X − E HF /X − 1 X exp(9 X − X − 1 ) X+1

−1

CCSD CCSD δ Ẽ HF /CBS = δ Ẽ HF /X CCSD CCSD ̃ δ E HF /X − δ Ẽ HF /(X − 1) + 3 X −1 X−1

(

)

(3)

%TAE[(T )] =

(4)

CCSD(T ) CCSD(T ) δ Ẽ CCSD /X − δ Ẽ CCSD /(X − 1) 3

( X X− 1 )

TAE(CCSD(T )) − TAE(CCSD) TAE(CCSD(T ))

(8)

where TAE is the total atomization energy excluding any vibrational and thermal contributions. A wave function stability analysis was performed for all UHF reference wave functions prior to the CCSD(T) calculation. Stability analysis verifies that the lowest (ground state) UHF reference wave function has been found. TAEs denoted “RHF” employed an ROHF reference for open-shell atoms. 2.2. Density Functional Theory Calculations. All DFT calculations were performed with the Gaussian09 electronic structure package. Two types of calculations were performed: (1) single point energy calculations on the CCSD-optimized geometries, which were also used for the best estimate energies

CCSD(T ) CCSD(T ) δ Ẽ CCSD /CBS = δ Ẽ CCSD /X

+

(7)

CC3

(2)

−1 (5)

where the consecutive cardinal numbers X − 1 and X are the maximum angular momentum quantum number X represented E


F

−32.42 −31.94 N/A N/A

19.68 22.01

15.48 19.94

N/A N/A

CBSc CBSd

7.29 7.72

6.48 7.19

6.78 7.83 5.08 7.41 6.00 0.63 6.76 13.12 −5.83

−3.35 −3.22 −4.46 −1.84 −4.12 0.24 −0.21

−0.35 −0.43 −0.16 −0.25 −0.27 −0.05 −2.12 −0.65

δECCSD(T) CCSD

N/A N/A

CBSc APVTZ

CBSc CBSc CBSc CBSc CBSc CBSc APVTZ APV5Z CBSc

CBSc CBSc CBSc CBSc CBSc CBSc APVTZ

CBSc CBSc CBSc CBSc CBSc CBSc CBSc APVTZ

basis

0.75 0.40

0.64 0.65

1.57 1.68 1.26 1.85 1.59 −0.79 1.00 1.35 −0.59

0.04 0.07 −0.06 −1.37 −0.56 −0.97 −1.58

−0.02 −0.02 −0.02 −0.02 −0.04 −0.11 −0.61 −0.04

δECC3 CCSD(T)

APVTZ mixed

APVTZ PVTZ

APVTZ APVTZ APVTZ APVTZ APVTZ APVTZ PVTZ APVTZ APVTZ

APVTZ APVTZ APVTZ APVTZ APVTZ APVTZ PVTZ

APVTZ APVTZ APVTZ APVTZ APVTZ APVTZ APVTZ PVTZ

basis

−4.21 −4.21

−3.74 −

−1.21 −1.43 −1.01 −1.81 −1.70 0.66 − −4.21 −

0.94 0.52 0.96 1.31 0.74 0.81 −

0.13 0.10 0.04 0.04 0.04 0.12 1.01 −

δECCSDT CC3

MAPVDZ MAPVDZ

MAPVDZ −

APVDZ MAPVDZ MAPVDZ MAPVDZ MAPVDZ APVDZ − MAPVDZ −

APVDZ MAPVDZ MAPVDZ MAPVDZ MAPVDZ APVDZ −

APVDZ MAPVDZ MAPVDZ MAPVDZ MAPVDZ APVDZ APVDZ −

basis

−0.62 −0.62

−1.19 −

1.92 2.11 1.74 2.06 1.85 0.50 − −0.62 −

−0.36 −0.43 −0.51 −1.64 −0.83 −0.05 −

−0.04 −0.06 0.01 0.00 −0.07 0.13 −0.38 −

δECCSDT(Q) CCSDT

MAPVDZ MAPVDZ

PVDZ −

PVDZ PVDZ PVDZ PVDZ PVDZ PVDZ − MAPVDZ −

PVDZ PVDZ PVDZ PVDZ PVDZ PVDZ −

PVDZ PVDZ PVDZ PVDZ PVDZ PVDZ PVDZ −

basis

0.16 0.16

−0.09 −

−0.19 −0.07 −0.07 −0.04 0.02 −0.19 − 0.16 −

0.09 0.08 0.16 0.01 0.08 −0.15 −

−0.01 −0.01 −0.02 −0.01 −0.01 −0.07 −0.07 −

δEcor.‑val. CCSD(T)

PCVTZ PCVTZ

PCVTZ −

PCVTZ PCVTZ PCVTZ PCVTZ PCVTZ PCVTZ − PCVTZ −

PCVTZ PCVTZ PCVTZ PCVTZ PCVTZ PCVTZ −

PCVTZ PCVTZ PCVTZ PCVTZ PCVTZ PCVTZ PCVTZ −

basis

−9.20 −12.09

Crrtd. est.

15.72 ± 1.09 −4.64

−63.01 ± 1.00 −56.19 ± 0.94 −19.39 ± 1.02 −12.12 ± 0.93 −4.99 ± 1.58 5.15 ± 0.78 −7.52 49.83 ± 0.36 −59.03

7.65 ± 0.47 3.01 ± 0.51 17.92 ± 0.53 23.13 ± 0.92 20.23 ± 0.86 9.73 ± 0.75 10.33

−1.83 ± 0.07 −2.01 ± 0.11 −1.96 ± 0.06 −2.75 ± 0.05 −2.23 ± 0.13 −4.87 ± 0.14 −6.13 ± 0.41 −3.80

best est.b

1 b Total ΔE values are given by the energy difference at the HF level, ΔEHF, and different correlation energy contributions, with respect to that difference, δElevel level 2. The uncertainty is given as ΔE ± 2σtot (a 95% probability interval), where σtot is a simplistic estimate of the total standard deviation of ΔE, based on the focal point scheme. See Section S3 in the SI. cAPVTZ and APVQZ basis sets are used for the extrapolations. dMAPVTZ and MAPVQZ basis sets are used for the extrapolations.

a

TMA: 1O2 + TMA-O TMA product

corrected values (see text)

CBSc CBSd

−1.87 −32.42

Br−: 1O2 + BrO− TMA: 1O2 + TMA-O

CBSc CBSc CBSc CBSc CBSc CBSc CBSd APV6Z CBSc

CBSc CBSc CBSc CBSc CBSc CBSc CBSd APV6Z CBSc

−87.15 −82.06 −38.05 −38.02 −25.80 4.49 −32.90 −32.45 0.03

C2H2 product C2H4 product HCN product NH3 product HCl product Br− product TMA product O3 → 1O2 + O(3P) TMA + O(3P) → TMA-O

15.26 15.76 11.66 16.43 13.05 −0.14 17.61 72.32 −52.63

CBSc CBSc CBSc CBSc CBSc CBSc CBSd

−12.53 −11.84 −13.23 −1.83 −9.26 −0.68 0.82

CBSc CBSc CBSc CBSc CBSc CBSc CBSd

22.84 17.82 35.07 28.49 34.19 10.52 11.30

C2H2 TS C2H4 TS HCN TS NH3 TS HCl TS Br− TS TMA TS

basis CBSc CBSc CBSc CBSc CBSc CBSc CBSc CBSd

δECCSD HF −1.74 −2.16 −0.97 −1.51 −1.16 −0.83 −6.56 −3.10

basis CBSc CBSc CBSc CBSc CBSc CBSc CBSc CBSd

0.20 0.57 −0.85 −1.00 −0.71 −4.06 2.59 −0.01

ΔEHF

Best estimates C2H2 vdW C2H4 vdW HCN vdW NH3 vdW HCl vdW Br− vdW Br− vdW (2) TMA vdW

species

Table 1. Best Estimates of ΔEvdW, ΔETS, ΔErxn, ΔE1O2 [kcal mol−1]a

The Journal of Physical Chemistry A Article


Article


Figure 2. Convergence of the absolute electronic energies within the CC hierarchy using the PVDZ basis set. (a) Ozone using RHF and UHF reference wave functions; zero was set to the RHF-CC5 result. (b) 1O2 using a RHF reference wave function; zero was set to the CCSDTQP result. (c) C2H2 using an RHF reference wave function; zero was set to the CCSDTQ(P) result. Note the 2-fold change in scale of the y-axis. (d) O(3P) using a UHF reference wave function; zero was set to the U-CC5 result. Note the 10-fold change in scale of the y-axis. Source data in Table S18 in the Supporting Information.

Figure 3. Convergence of electronic reaction energies to isolated reactants (O3, C2H2) within the CC hierarchy using RHF and UHF reference wave functions with the PVDZ basis set. (a) O3 + C2H2 → C2H2O3, where for each energy (ΔEvdW, ΔETS, and ΔErxn), zero was set to the mean of the (R/U)-CCSDT(Q) results. (b) O3 → 1O2 + O, where absolute reaction energies are given, the red line corresponds to the CC5 result. Source data are given in Table S18 in the Supporting Information.

D3,108,109 ωB97XD,76 LC-ωPBE,110 and M1177 functionals, as implemented in Gaussian09.

described above; (2) a reoptimization of the geometries with the respective DFT method, where the nature of the obtained stationary points was verified through a vibrational frequency analysis. All DFT calculations employed the APVTZ basis set. Single-point energy calculations used a pruned (99 590) grid (int = grid = ultrafine keyword in Gaussian09). Calculations used the PBE,92,93 OLYP,94,95 HCTH407,96 O3LYP,97 B3LYP,98 B97-2,99 PBE0,100 M06,101 M05,102 BMK,103 B2PLYP-D3, 23 BH&HLYP, 104 M06-2X, 101 N12SX, 105 MN12SX,105 HSEh1PBE,106 LC-PBE,92,93,107 CAM-B3LYP-

3. RESULTS AND DISCUSSION 3.1. Selection of the Types of Reactions Studied. We conducted calculations for seven molecules reacting with ozone. We considered a single reaction channel for each species (Figure 1): (1) concerted cycloadditions to a double/ triple bond (Criegee mechanism): ethyne, ethene, hydrogen cyanide, (2) monodentate addition with simultaneous hydroG


Article


Table 2. Comparison of Thermokinetic Data for the O3 + C2H2 and O3 + C2H4 Reactions Reported Here and Previously Published by Zhao et al.25 and Wheeler et al.16 Thermokinetic Data [kcal mol−1] C2H2 + O3 study this study Wheeler et al.16 this study Wheeler et al.16 this study this study Wheeler et al.16 Zhao et al.25 Zhao et al.25

model a

CCSD(T)/APVQZ CCSD(T)/APVQZa CCSD(T)/CBSb CCSD(T)/CBSb best estimate omitting CC3c best estimated CCSDT+(2) eQ CCSD(TQ)e

C2H4 + O3

ΔEvdW

ΔETS

ΔErxn

ΔEvdW

ΔETS

ΔErxn

−1.93 −2.03 −1.89 −2.02 −1.83 −1.83 −1.85 −1.88 −1.98

6.91 6.94 7.04 7.05 7.65 7.55 7.73 7.89 7.57

−64.70 −64.50 −65.29 −65.35 −63.01 −63.44 −63.17 −64.04 −64.59

−2.11 −2.20 −2.03 −2.04 −2.01 −1.99 −1.84 −1.94 −2.03

2.54 2.54 2.84 2.83 3.01 3.03 3.44 3.51 3.19

−58.23 −58.04 −58.54 −58.70 −56.19 −56.32 −56.53 −57.24 −57.95

a CC3 CCSDT d Including δEcore−valence bDifferent extrapolation schemes were applied. cUsing δECCSDT CCSD(T) instead of δECCSD(T) + δECC3 . Excluding ΔRel and ΔDBOC from that publication. eReplaced the quadruples term from Wheeler et al.16 with other methods. ΔRel and ΔDBOC were subtracted from the best estimates in the reference.

distance is too short for a mere vdW complex. For both TMA and Br−, we included the mono-oxides TMA-O and BrO− as well as 1O2 in the calculations, allowing us to report the dissociation energies from the initial product structures. 3.2. Coupled Cluster Energy Calculations. 3.2.1. Types of Energies Reported. Our aim was to provide accurate electronic energies for the stationary structures shown in Figure 1, which can then be used to evaluate the performances of various DFT methods. As such, the reported energies in the present study do not contain relativistic effects or Born− Oppenheimer corrections, since these contributions would not be included in DFT calculations either. Furthermore, we do not consider vibrational data (zero-point vibrational energies or temperature effects). We report four types of energies, all of which are defined relative to separated, optimized reactants: the binding energy of the vdW/prereactive complex (ΔEvdW), the energy of the transition structure (ΔETS), the energy of the (initial) product (ΔErxn), and the energy of separated product molecules of the monoxide and singlet oxygen in the cases of TMA and Br− (ΔE1O2). For our purposes, there is little need to report the electronic barrier height ΔE⧧ (the difference between ΔETS and ΔEvdW) explicitly. For all systems except Br−, ΔEvdW is accurately described by almost all electronic structure methods considered here. Hence, the accuracy of a method for ΔETS is closely related to that found for ΔE⧧. 3.2.2. Convergence of the Electron Correlation Energies of O3 and 1O2. Before discussing the ΔE results of the reactions of interest, we aimed to identify and quantify the sources of the errors in computed electronic energies of the problematic species. Most past theoretical studies on O3 chemistry employed CC calculations based on a restricted Hartree−Fock (RHF) wave function. This is a reliable approach to compute accurate electronic correlation energies of compounds not affected by nondynamical correlation, but may be inappropriate for systems with a high amount of nondynamical correlation, such as O3. In the case of a species with diradical character such as O3, the other straightforward approach is to use an unrestricted Hartree−Fock (UHF) reference wave function, which, however, will be heavily spincontaminated. In the CC calculations, this choice of reference wave function should become less important when including increasing levels of excitations, finally converging to the exact Full-CI solution in both cases.

gen transfer: ammonium, hydrochloric acid, and (3) monodentate addition/oxygen atom transfer: bromide, trimethylamine (TMA). These cover some of the functional groups (aliphatic amines, inorganic anions, and olefins) and reaction mechanisms that are widespread in aqueous ozone chemistry. This diverse set should enable a general overview of the energetics of many important ozone reactions. Because of their increased size, aromatic systems were not included in the present study. The considered reaction channels were found on the CCSD potential energy surface (PES); nuclear coordinates of all stationary structures are reported in Section S1 in the SI. These are not necessarily the only possible reaction channels for the respective systems: for example, a multitude of reaction channels involving NH3 and O3 are possible.26 For Br−, the employed CCSD geometry optimization located two van der Waals (vdW) complexes: one Cs-symmetric complex was found slightly more stable on the CCSD PES; the other (nonsymmetric) conformation, directly preceding the transition structure, was found slightly more stable on the CCSD(T) PES. Therefore, we computed energies for both conformations. The stationary structures for the C2H2 and C2H4 reactions correspond to those previously reported16 for the respective Criegee mechanisms, featuring a Cs-symmetric vdW complex and transition structure. The vdW complex of HCN is Cssymmetric, with the nitrogen pointing toward the central oxygen of O3. The transition structure of HCN is naturally not symmetric, and the O−N distance is 0.2 Å longer than the O− C distance. The transition structures of NH3 and HCl were initially found by performing a “relaxed scan”, that is, by systematically decreasing the O−N and O−Cl distances while reoptimizing the rest of the system. This resulted in a concerted hydrogen transfer reaction in both cases. Such a hydrogen transfer is not possible for Br− or TMA. In these cases, O3 simply adds to the Br/N atoms. The resulting products appear to be a precursor to the release of 1O2. In both cases, the corresponding O−O distance is elongated (2.11 Å for Br−, 2.87 Å for TMA). The elongated O−O distance for TMA corresponds to the dissociation into TMA−N-oxide (TMA-O) and 1O2, corroborated by the thermodynamic data reported below (Table 1), and there was no formation of an intermediate TMA−O−O−O−. In the case of Br−, the O−O H


Article


Table 3. Errors Relative to Best Estimates, Using CCSD Geometriesa Error [kcal mol−1] van der Waals Binding Energies (vdW) method CCSD/APVTZ + δECCSD(T) CCSD /PVDZ CCSD(T)/APVTZ CCSD(T)/APVTZ + δECC3 CCSD(T)/PVDZ

Transition Structures (TS)

Reaction Energies (rxn)

MUE

MAX

MUE

MAX

MUE

MAX

0.27

0.96

0.67

1.93

2.24

4.06

0.30 0.31

0.80 0.81

0.98 1.01

1.74 2.39

2.37 1.65

3.83 4.00

a

MUE is the mean unsigned error, whereas MAX is the maximum deviation. Errors are reported for vdW binding energies, transition structures (TS), and reaction energies (rxn), which include the 1O2cleaved structures.

energy, with respect to excitation level for the given basis set. This approach is additionally justified by the fact that postCCSD(T) contributions to the atomization energy, even of multireference systems, are basis-set converged to within ∼1 kcal mol−1 with the PVDZ basis.78,111 In Figure 2a, we compared the absolute energy of the O3 molecule computed with increasing CC levels up to CCSDTQ(P) with the PVDZ basis. From the CCSDT level onward, there is little difference between RHF-based and UHF-based energies, with the largest discrepancy of 0.7 kcal/mol observed at the CCSDT(Q) level. However, the RHF-based absolute energies at the CC3 and CCSDT(Q) level are closer to the (deemed most accurate) RCC5 result. Karton et al.78 reported the δECCSDTQ CCSDT(Q)/PVDZ and δECCSDTQ(P) contributions to the atomization energy of O3, CCSDTQ only with a UHF reference, to be −0.8 and 0.5 kcal mol−1, respectively. This agrees well with the differences in electronic energy reported in Figure 2, and we stress that the RCCSDT(Q) result is very close to R/U-CCSDTQ(P) and RCC5 already. Karton et al. also explored the use of RHF vs UHF reference wave functions for total atomization energies at the CCSD(T) and CCSDT levels with the larger pCVTZ basis,78 and found that for O3, this yielded differences of 0.2 kcal mol−1 to 0.3 kcal mol−1 in both cases. However, those data are not directly comparable to the R-/U-CC electronic energies in Figure 2, since calculated atomization energies in this previous study also differed by the use of UHF versus ROHF references for atomic oxygen. In two of the computed reaction channels of the present study, 1O2 is a product, which is, in turn, a multireference system (like O3) that is problematic for single-reference WF methods.112,113 Figure 2b shows that the convergence of the correlation energy of 1O2 within the CC series is not only slow and irregular, but also rather different from that of O3. For O3 and 1O2, CCSDT(Q) is necessary to come within 1 kcal mol−1 of the converged result. For comparison, in the well-behaved systems C2H2 (very minor multireference case) and O(3P), the convergence of the correlation energy is fast: absolute energies are within 1 kcal mol−1 of the converged result already at the CCSD(T) level (Figure 2). The CC3/CC4 methods are usually not used in thermochemical calculations, but we found them of interest for studying the convergence of the correlation energy in the CC hierarchy. For well-behaved systems (C2H2 and O(3P)), the resulting correlation energy is almost identical to that of the perturbative inclusion of excitations of the same order

Figure 4. Performance of approximate CC methods for the different energies (ΔEvdW, ΔETS, ΔErxn), relative to the best estimates (ΔΔE = ΔE(method) − ΔE(best)). Frozen CCSD geometries were used.

From a practical point of view, restricted calculations are computationally less expensive and are, thus, preferred for larger systems. We studied the difference between RHF-CC and UHF-CC calculations with the small PVDZ basis, omitting basis-set convergence. This yields inaccurate absolute energies, but allows for an evaluation of convergence of the correlation I


Article


Figure 5. Performance of approximate methods for the different energies (ΔEvdW, ΔETS, ΔErxn), relative to the best estimates. Frozen CCSD geometries were used. Local functionals do not include orbital exchange, and hybrid functionals are ordered by increasing fraction of orbital exchange. SRC functionals employ orbital exchange in the short range, and LRC functionals employ orbital exchange in the long range.

treatment. Consequently, chemical accuracy (≤1 kcal mol−1 error) of reaction energies can often be reached with CCSD(T).73 However, in the case of O3 and 1O2, the discrepancy between the correlation energy and the converged result at the lower CC levels (

Gas-Phase Ozone Reactions with a Structurally Diverse Set of

Recommend Documents