THE JOURNAL OF CHEMICAL PHYSICS 139, 201103 (2013)
Communication: An accurate global potential energy surface for the ground electronic state of ozone Richard Dawes,1,a) Phalgun Lolur,1 Anyang Li,2 Bin Jiang,2 and Hua Guo2,a) 1
Department of Chemistry, Missouri University of Science and Technology, Rolla, Missouri 65409, USA Department of Chemistry and Chemical Biology, University of New Mexico, Albuquerque, New Mexico 87131, USA 2
(Received 22 October 2013; accepted 17 November 2013; published online 27 November 2013) We report a new full-dimensional and global potential energy surface (PES) for the O + O2 → O3 ozone forming reaction based on explicitly correlated multireference configuration interaction (MRCI-F12) data. It extends our previous [R. Dawes, P. Lolur, J. Ma, and H. Guo, J. Chem. Phys. 135, 081102 (2011)] dynamically weighted multistate MRCI calculations of the asymptotic region which showed the widely found submerged reef along the minimum energy path to be the spurious result of an avoided crossing with an excited state. A spin-orbit correction was added and the PES tends asymptotically to the recently developed long-range electrostatic model of Lepers et al. [J. Chem. Phys. 137, 234305 (2012)]. This PES features: (1) excellent equilibrium structural parameters, (2) good agreement with experimental vibrational levels, (3) accurate dissociation energy, and (4) most-notably, a transition region without a spurious reef. The new PES is expected to allow insight into the still unresolved issues surrounding the kinetics, dynamics, and isotope signature of ozone. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4837175] Ozone has been extensively studied for many years both experimentally and theoretically due to its central importance to atmospheric chemistry. Understanding ozone in terms of its spectroscopy, kinetics, and dynamics has been very challenging for theory despite being only a triatomic. There is still a substantial discrepancy between predicted and observed low-temperature rates of O + O2 isotope exchange.1–6 Experimental rates are 3-5 times larger than predictions at low-temperature and exhibit a negative temperature dependence that has proven difficult to reproduce theoretically. Another outstanding issue of considerable ongoing interest is the anomalous isotope enrichment in atmospheric ozone, the so called mass-independent fractionation (MIF) effect.7–10 Perhaps the biggest limitation for achieving agreement between theory and experiment has been the lack of an entirely satisfactory potential energy surface or surfaces (PESs).11–16 Here, we focus on the development of an accurate PES for the ground state based on explicitly correlated multireference configuration interaction (MRCI-F12)17, 18 data computed using MOLPRO,19 which can be used to make accurate predictions of the kinetics and dynamics. The electronic structure of ozone is very challenging to describe accurately with ab initio calculations. It has been noted as a problematic case for highly accurate singlereference-based schemes to calculate atomization energies.20 The best current experimental value for the dissociation energy from the latest active thermochemical tables (ATcT) is D0 = 8563 ± 3.5 cm−1 , which corresponds to De ∼ 9219 ± 10.0 cm−1 (we assign larger uncertainty to the zero-point energy, or ZPE of O3 ).21 As we have reported previously, De for standard MRCI varies with basis set from 7927 cm−1 a) Electronic addresses: [email protected] and [email protected]
0021-9606/2013/139(20)/201103/4/$30.00
at the AVTZ level to 8768 cm−1 at the AV6Z level, an increase of more than 800 cm−1 .15 This makes complete basis set (CBS) extrapolation uncertain and we previously noted a variation of ∼150 cm−1 between various CBS schemes.15 The large effect of the Davidson correction (>300 cm−1 ) on De is evidence of the importance of high-order correlation.22 Dynamic correlation is very important in general, since De at the full-valence complete active space self-consistent field (CASSCF) level is only 3003 and 3104 cm−1 for the VTZ-F12 and VQZ-F12 basis sets,23 respectively, and the shape of the minimum energy path (MEP) towards dissociation changes qualitatively with the inclusion of dynamic correlation. Watts and Bartlett24 found that for single reference coupled cluster, the effect of adding triple excitations, contributed changes in the two harmonic stretching frequencies ω1 and ω3 of more than 100 cm−1 , a rather dramatic effect of high-order correlation. For ozone, we find that at geometries away from the global minimum the system becomes increasingly multiconfigurational. This puts single-reference methods out of the question for a global PES making a multireference method such as MRCI the appropriate choice. There are also reasons to consider multiple electronic states. Adiabatic dissociation of ozone 1 3 O3 A1 → O 3 g− +O Pg , leads to a threefold degeneracy of singlet states (2 1 A + 1 A in the Cs group). Diabatically, the ground state of ozone correlates to excited states of both the atomic and molecular oxygen products 1 O3 (1 A1 ) → O2 1 g + O Dg ,
1
yielding another block of 10 degenerate singlet states (5 1 A + 5 1 A in the Cs group). 139, 201103-1
© 2013 AIP Publishing LLC
201103-2
Dawes et al.
A critical issue for the PES of ozone in the bottleneck region is an avoided crossing with an excited state that can result in disruptions in the ground state surface seen as a submerged reef and an associated van der Waals (vdW) well in the exit channel. The electronic structure in the bottleneck region was the main subject of our previous study with many details provided in the supplementary material of that paper.15 To briefly summarize, we found that the height of the submerged reef and depth of the vdW well are substantially reduced by both basis set and active space completeness. We followed the switch in character between the ground state and an excited state correlated with O2 (1 g ) + O(1 Dg ) products and noted that if all 13 singlet states were included using a dynamically weighted state-averaged method25, 26 (DW-SACASSCF), then the transition was better behaved and the reef all but disappeared. Finally, by including more than one reference state in the internally contracted icMRCI (thus increasing the number of contracted configurations and reducing the effect of internal contraction), we obtained a MEP that is a monotonic function of the distance to the dissociating O-atom (no reef feature). Given the sensitivity of the kinetics (noted in our previous study) to the shape of the PES in the bottleneck region, our goal has been to construct a full-dimensional global PES useful to explore the unresolved issues surrounding the kinetics and dynamics of this important system. We identified four features desired in a PES suitable for dynamics calculations: (1) excellent equilibrium structural parameters, (2) good agreement with experimental vibrational levels, (3) accurate dissociation energy, and (4) most-notably, a transition region with accurate topography (especially without a spurious reef). Until recently, the best MRCI-based PES for ozone was that initially developed by Siebert, Schinke, and Bittererova (SSB) in 2002.11, 12 The SSB PES is nominally based on 5000 spline-fit points at the MRCI/VQZ level starting from a one-state CASSCF calculation with the 2s orbitals closed in the reference space (12,9). As discussed above, direct use of MRCI at that basis level substantially underestimates the dissociation energy, so an adjusted hybrid PES was constructed.13 Schinke et al.4 were also aware of the sensitivity of the kinetics to the height of the pronounced reef found in their ab initio data and some numerical experiments were reported using versions of the PES in which the reef was removed. The low-lying vibrational levels on the SSB PES are very realistic matching experiment to a root-mean-squared error (RMSE) of ∼ 4 cm−1 for levels up to 4000 cm−1 above the ZPE.12 Straightforward application of the MRCI method at this basis set level does not perform so well and the recent (2013) PES of Ayouz and Babikov27 employing the same (12,9) reference space, but extrapolating the MRCI data to the CBS limit (MRCI/CBS) overshoots the fundamentals ν 1 and ν 3 by roughly 20 and 30 cm−1 , respectively, with larger errors recorded for overtones and combinations. Those results closely match our unpublished tests in which we have produced accurate fits at the MRCI/CBS level. The new Babikov PES includes the full symmetry and has an accurate dissociation energy of 9177 cm−1 (close to the experimental value of 9219 ± 10.0 cm−1 discussed above). However, it lacks a spin-orbit (SO) correction which would
J. Chem. Phys. 139, 201103 (2013)
lower the asymptote by ∼80 cm−1 (pushing De further from experiment) and influence the shape of the long-range PES. It also has a pronounced reef feature which means that overall it meets two of our four above-stated criteria for our new PES. Very recently, Tyuterev et al.28 published a spectroscopic PES fit to an analytic form describing one of the three global minimum isomer wells. This PES uses one-state Davidson corrected MRCI data starting with a full-valence (18,12) CASSCF reference. The PES combines data generated with the AV5Z basis (in the well), with other data extrapolated to CBS(5,6) (along the dissociation coordinate) thus achieving a dissociation energy of De = 9164 cm−1 (similar to that obtained by Babikov27 ). The authors report the submerged reef feature in their data, but produced two versions of their PES including one for which a correction was introduced based on our previous report (thus removing the reef). Interestingly, the presence or absence of the reef was found to influence the high-lying vibrational levels. Ultimately, they arrived at a PES without a reef feature that matches the experimental vibrational levels including high-lying levels up to ∼7900 cm−1 above the ZPE to within about 1 cm−1 RMSE.28 The highest lying experimental levels around 7900 cm−1 are only about 650 cm−1 away from the experimental value for D0 of 8563 ± 3.5 cm−1 .21 The Tyuterev PES undershoots D0 by only ∼55 cm−1 , but adding a correction for spin-orbit (not yet included) would affect the shape of the bottleneck and asymptotic region and bring D0 on the PES to ∼135 cm−1 below experiment. Thus, if higher-lying experimental levels were known, the error in calculated levels would have to increase rapidly from sub-wavenumber up to 7900 cm−1 above ZPE, to more than 100 cm−1 over the next ∼500 cm−1 toward dissociation. The Tyuterev PES is already proving very useful for understanding the spectroscopy of the known levels including determining the states involved in particular resonance polyads and helping to determine the terms and parameters in effective spectroscopic Hamiltonians.29 All of which aid in the development of line-lists for atmospheric databases such as HITRAN.30 Thus, the Tyuterev PES meets at least 3 of our 4 criteria desired for our planned studies of kinetics and dynamics, and will likely be useful for future comparisons. In view of our previous study and the observed sensitivity of the shape of the MEP and dissociation energy to basis set completeness and correlation treatment, we tested the recently developed MRCI-F12 method along the MEP. Explicitly correlated F12 methods have proven to greatly improve convergence with basis set, and been shown to produce near CBS quality results with relatively small basis sets.31 Figure S1 in the supplementary material32 plots a 1D cut using MRCI-F12 showing that the pronounced reef seen at the triple-zeta level for standard MRCI is nearly absent in the MRCI-F12 data, and at the quadruple-zeta level a monotonic cut is obtained (a result only achieved at the CBS level for standard MRCI). The dissociation energy De increases from 9162 cm−1 at the VTZF12 level to 9354 cm−1 with VQZ-F12 (without spin-orbit or other corrections). This is a much smaller step than found in standard MRCI for which De increases by nearly 700 cm−1 between AVTZ and AVQZ. The value for De of 9354 cm−1 obtained with VQZ-F12 already slightly exceeds experiment. Adding a correction for spin-orbit (which we do) lowers De by
201103-3
Dawes et al.
FIG. 1. Contour plot of the asymptotic region with rOO fixed at 2.282 bohrs, with R = [3.5, 12] bohrs and θ = [0◦ , 90◦ ]. The contours are in 10 cm−1 steps, the one labeled 9270 cm−1 would be the last to appear at any distance since De = 9275.12 cm−1 on the fitted PES with SO-correction.
∼80 cm−1 and adding a further scalar relativistic correction (which we neglect), lowers De by another 21 cm−1 ,15 reaching De = 9253 cm−1 (just above the estimated experimental value of 9219 ± 10.0 cm−1 discussed above). Feller and Peterson33 have reported systematic overestimation of atomization energies in some systems by F12 methods when extrapolated to the CBS limit consistent with our results. They also conclude however, that unextrapolated F12 data at the VQZF12 basis level “provides a very cost effective alternative to achieving results comparable to much larger standard method calculations.” Thus, we fit a global PES using MRCI-F12/VQZ-F12 data without CBS extrapolation (some additional details are given in the supplementary material).32 Our automated PES generation method described previously34–36 was used to fit the surface out to R(O–OO) = 10 bohrs in Jacobi coordinates including 2663 ab initio data (some obtained by permutation). No energy-dependent bias was used for the point selection, so the four wells (including the high-energy ring-minimum) and the asymptotic region were treated equally. To include the long-range part of the PES, we switch after 8 bohrs to the 2D (R, θ ) long-range interaction model of O2 + O by Lepers et al.37 To describe the r-dependence of the O2 fragment asymptotically, we switch simultaneously to an accurate potential energy curve for O2 by Ruedenberg et al.38 We also added a correction for SO coupling (some additional details are given in the supplementary material).32 The long-range model of Lepers et al.37 has a similar optional SOcorrection, so conveniently a consistent switch could be made to the long-range, with or without the SO-correction. No reef or vdW well is observed with or without the SO-correction. Figure 1 shows a contour plot of the PES with SO-correction in the asymptotic region. The C2V global minimum on the fitted PES has geometric parameters of re = 2.4031 a.u. [r0 (expt) = 2.4052 a.u.], and θ e = 116.84◦ [θ 0 (expt) = 116.75◦ ]39, 40 while the D3h ring minimum is found at an optimized structure of re = 2.715 a.u. (a value very similar to previous reports),12 with an energy of 10 755.1 cm−1 (1580 cm−1 above De , see the supplementary material32 ). We calculated the low-lying vibrational levels (up to 4000 cm−1 above ZPE) for the 16 O3 and 18 O3 mass combi-
J. Chem. Phys. 139, 201103 (2013)
FIG. 2. Cut along the bond stretching coordinate in valence coordinates comparing our new MRCI-F12 PES with the Babikov 2013 PES.27 The other bond is fixed at r1 = 2.282 a.u. and the angle is fixed at θ = 116.75◦ .
nations localized in one well, using the parallel ro-vibrational DVR code RTR3D by Wang and Carrington.41 Reasonably accurate vibrational levels are produced (RMSE for 16 O3 and 18 O3 are 14.4 and 7.8 cm−1 , respectively) that do not degrade substantially toward higher energy. They are not as accurate as the spectroscopic PES of Tyuterev,28 but indicate good representation of the wells in a PES intended for studies of dynamics. Tyuterev et al.28 report some influence of the reef (or its removal) even on levels 600-700 cm−1 below dissociation. Grebenshchikov et al.42 previously studied the highest lying states on a PES with the reef and reported many states localized in the vdW wells, and many others exhibiting mixed or delocalized character. In the supplementary material,32 we compare all of the bound states for the 16 O17 O18 O mass combination using the Babikov 2013 PES27 (with reef) and ours (without reef). Figure 2 shows a simple unrelaxed cut along the bond stretching coordinate in valence coordinates comparing our new MRCI-F12 PES with the Babikov 2013 PES. The pronounced reef feature in the Babikov PES changes the high-lying vibrational levels qualitatively (see text in the supplementary material),32 while for our new PES, in addition to supporting more bound states, even the high-lying states are mostly recognizable from one-well calculations. We applied a three-dimensional (3D) quantum statistic model (QSM) to predict the temperature dependent rate of the O + O2 exchange reaction, extending our previous twodimensional (2D) work. Figure 3 compares 3D results using the new PES to our previous 2D results and experiment. As seen in Figure 3, the rate predicted by the new 3D calculation is slightly lower than the 2D result. It underestimates the experiment at low temperatures. However, the negative temperature dependence qualitatively reproduced in the new 3D results is a direct result of the reef-less MEP. It is known that the QSM is a crude approximation to the dynamics and not strictly applicable in this case because of the known nonstatistical behaviors of this reaction.6, 43 Nevertheless, we take these results as an indication that more sophisticated stateto-state scattering calculations might be expected to display the desired negative temperature dependence on this PES, and hopefully produce much better agreement overall. In addition, a time-dependent wavepacket calculation for J = 0 indicates a rich resonance structure shown in Figure S6 of the
201103-4
Dawes et al.
J. Chem. Phys. 139, 201103 (2013) 5 S.
Y. Lin and H. Guo, J. Phys. Chem. A 110, 5305 (2006). Sun, L. Liu, S. Y. Lin, R. Schinke, H. Guo, and D. H. Zhang, Proc. Natl. Acad. Sci. U.S.A. 107, 555 (2010). 7 M. H. Thiemens, Science 293, 226 (2001). 8 M. H. Thiemens, Annu. Rev. Earth Planet. Sci. 34, 217 (2006). 9 K. Mauersberger, D. Krankowsky, C. Janssen, and R. Schinke, Adv. At., Mol., Opt. Phys. 50, 1 (2005). 10 Y. Q. Gao and R. A. Marcus, Science 293, 259 (2001). 11 R. Siebert, R. Schinke, and M. Bittererova, Phys. Chem. Chem. Phys. 3, 1795 (2001). 12 R. Siebert, P. Fleurat-Lessard, R. Schinke, M. Bittererova, and S. C. J. Farantos, J. Chem. Phys. 116, 9749 (2002). 13 D. Babikov, B. K. Kendrick, R. B. Walker, R. T. Pack, P. Fleurat-Lessard, and R. J. Schinke, J. Chem. Phys. 118, 6298 (2003). 14 F. Holka, P. G. Szalay, T. Müller, and Vl. G. Tyuterev, J. Phys. Chem. A 114, 9927 (2010). 15 R. Dawes, P. Lolur, J. Ma, and H. Guo, J. Chem. Phys. 135, 081102 (2011). 16 D. Xie, H. Guo, and K. A. Peterson, J. Chem. Phys. 112, 8378 (2000). 17 T. Shiozaki, G. Knizia, and H.-J. Werner, J. Chem. Phys. 134, 034113 (2011). 18 T. Shiozaki and H.-J. Werner, Mol. Phys. 111, 607 (2013). 19 H.-J. Werner, P. J. Knowles, G. Knizia et al., MOLPRO , version 2012.1, a package of ab initio programs, 2012, see http://www.molpro.net. 20 D. Feller, K. A. Peterson, and D. A. Dixon, J. Chem. Phys. 129, 204105 (2008). 21 B. Ruscic, private communication of ATcT result based on version 1.110 of the Core (Argonne) Thermochemical Network (17 May, 2010). 22 Y. G. Khait, W. Jiang, and M. R. Hoffmann, Chem. Phys. Lett. 493, 1 (2010). 23 K. A. Peterson, T. B. Adler, and H.-J. Werner, J. Chem. Phys. 128, 084102 (2008). 24 J. D. Watts and R. J. Bartlett, J. Chem. Phys. 108, 2511 (1998). 25 R. Dawes, A. W. Jasper, C. Tao, C. Richmond, C. Mukarakate, S. H. Kable, and S. A. Reid, J. Phys. Chem. Lett. 1, 641 (2010). 26 M. P. Deskevich, D. J. Nesbitt, and H.-J. Werner, J. Chem. Phys. 120, 7281 (2004). 27 M. Ayouz and D. Babikov, J. Chem. Phys. 138, 164311 (2013). 28 V. G. Tyuterev, R. V. Kochanov, S. A. Tashkun, F. Holka, and P. G. Szalay, J. Chem. Phys. 139, 134307 (2013). 29 A. Barbe, S. Mikhailenko, E. Starikova, M.-R. DeBacker, V. G. Tyuterev, D. Mondelain, S. Kassi, A. Campargue, C. Janssen, S. Tashkun, R. Kochanov, R. Gamache, and J. Orphal, J. Quant. Spectrosc. Radiat. Transf. 130, 172 (2013). 30 L. S. Rothman, I. E. Gordon, Y. Babikov, A. Barbe, D. C. Benner, P. F. Bernath, M. Birk, L. Bizzocchi, V. Boudon, L. R. Brown et al., J. Quant. Spectrosc. Radiat. Transf. 130, 4 (2013). 31 H.-J. Werner, T. B. Adler, G. Knizia, and F. R. Manby, in Recent Progress in Coupled-Cluster Methods, edited by P. Cársky, J. Paldus, and J. Pittner (Springer, 2010). 32 See supplementary material at http://dx.doi.org/10.1063/1.4837175 for more details of electronic structure calculations, vibrational levels and QSM and scattering calculations. 33 D. Feller and K. A. Peterson, J. Chem. Phys. 139, 084110 (2013). 34 R. Dawes, X.-G. Wang, A. W. Jasper, and T. Carrington, Jr., J. Chem. Phys. 133, 134304 (2010). 35 R. Dawes, D. L. Thompson, A. F. Wagner, and M. Minkoff, J. Phys. Chem. A 113, 4709 (2009). 36 R. Dawes, D. L. Thompson, A. F. Wagner, and M. Minkoff, J. Chem. Phys. 128, 084107 (2008). 37 M. Lepers, B. Bussery-Honvault, and O. Dulieu, J. Chem. Phys. 137, 234305 (2012). 38 L. Bytautas, N. Matsunaga, and K. Ruedenberg, J. Chem. Phys. 132, 074307 (2010). 39 Vl. G. Tyuterev, S. Tashkun, P. Jensen, and A. Barbe, J. Mol. Spectrosc. 198, 57 (1999). 40 Vl. G. Tyuterev, S. Tashkun, D. W. Schwenke, P. Jensen, T. Cours, A. Barbe, and M. Jacon, Chem. Phys. Lett. 316, 271 (2000). 41 RTR is a package of programs to compute Rovibrational levels and wavefunctions of TRiatomic molecules, X.-G. Wang and T. Carrington, Jr. 42 S. Y. Grebenshchikov, R. Schinke, P. Fleurat-Lessard, and M. Joyeux, J. Chem. Phys. 119, 6512 (2003). 43 A. L. Van Wyngarden, K. A. Mar, K. A. Boering, J. J. Lin, Y. T. Lee, S.-Y. Lin, H. Guo, and G. Lendvay, J. Am. Chem. Soc. 129, 2866 (2007). 6 Z.
FIG. 3. Comparison of 3D QSM calculation for rate of exchange reaction (O + O2 ), shown in black solid line, with previous 2D result (red curve) and experiment (dashed curve).
supplementary material32 which compares results for three PESs including our new one. To summarize, we report a new PES for the ground state of ozone based on explicitly correlated MRCI-F12 data with four key features: (1) excellent equilibrium structural parameters, (2) good agreement with experimental vibrational levels, (3) accurate dissociation energy, and (4) most-notably, a transition region with accurate topography (especially without a spurious reef). D0 = 8610.9 cm−1 on the new PES exceed the best experimental value of 8563.6 cm−1 by just 47.3 cm−1 , which could be improved another 21 cm−1 by including a so far neglected scalar relativistic correction. The PES includes a spin-orbit correction and switches to a longrange electrostatic model defined out to infinity. The absence of the submerged barrier has important implications for kinetics. As we have demonstrated with an approximate quantum scattering model, the new PES brings the calculated rate constant for the O + O2 exchange reaction into much better agreement with experiment and produces the negative temperature dependence observed experimentally. The new 3D PES will be useful for quantitative characterization of the O + O2 exchange reaction and ozone formation in full dimensionality. This work was supported by (U.S.) Department of Energy (DOE) (DE-FG02-05ER15694 to H.G.) and National Science Foundation (NSF) (CHE-1300945 to R.D.). R.D. thanks Maxence Lepers, Béatrice Bussery-Honvault, and Olivier Dulieu for a Fortran coding of their long range electrostatic model and Vladimir Tyuterev and Dmitri Babikov for useful discussions. A Fortran coding of the PES is available upon request from R. Dawes ([email protected]). 1 S.
M. Anderson, F. S. Klein, and F. Kaufman, J. Chem. Phys. 83, 1648 (1985). 2 M. R. Wiegell, N. W. Larsen, T. Pedersen, and H. Egsgaard, Int. J. Chem. Kinet. 29, 745 (1997). 3 P. Fleurat-Lessard, S. Y. Grebenshchikov, R. Schinke, C. Janssen, and D. Krankowsky, J. Chem. Phys. 119, 4700 (2003). 4 P. Fleurat-Lessard, S. Y. Grebenshchikov, R. Siebert, R. Schinke, and N. Halberstadt, J. Chem. Phys. 118, 610 (2003).
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
Communication: An accurate global potential energy surface for the ground electronic state of ozone Richard Dawes,1,a) Phalgun Lolur,1 Anyang Li,2 Bin Jiang,2 and Hua Guo2,a) 1
Department of Chemistry, Missouri University of Science and Technology, Rolla, Missouri 65409, USA Department of Chemistry and Chemical Biology, University of New Mexico, Albuquerque, New Mexico 87131, USA 2
(Received 22 October 2013; accepted 17 November 2013; published online 27 November 2013) We report a new full-dimensional and global potential energy surface (PES) for the O + O2 → O3 ozone forming reaction based on explicitly correlated multireference configuration interaction (MRCI-F12) data. It extends our previous [R. Dawes, P. Lolur, J. Ma, and H. Guo, J. Chem. Phys. 135, 081102 (2011)] dynamically weighted multistate MRCI calculations of the asymptotic region which showed the widely found submerged reef along the minimum energy path to be the spurious result of an avoided crossing with an excited state. A spin-orbit correction was added and the PES tends asymptotically to the recently developed long-range electrostatic model of Lepers et al. [J. Chem. Phys. 137, 234305 (2012)]. This PES features: (1) excellent equilibrium structural parameters, (2) good agreement with experimental vibrational levels, (3) accurate dissociation energy, and (4) most-notably, a transition region without a spurious reef. The new PES is expected to allow insight into the still unresolved issues surrounding the kinetics, dynamics, and isotope signature of ozone. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4837175] Ozone has been extensively studied for many years both experimentally and theoretically due to its central importance to atmospheric chemistry. Understanding ozone in terms of its spectroscopy, kinetics, and dynamics has been very challenging for theory despite being only a triatomic. There is still a substantial discrepancy between predicted and observed low-temperature rates of O + O2 isotope exchange.1–6 Experimental rates are 3-5 times larger than predictions at low-temperature and exhibit a negative temperature dependence that has proven difficult to reproduce theoretically. Another outstanding issue of considerable ongoing interest is the anomalous isotope enrichment in atmospheric ozone, the so called mass-independent fractionation (MIF) effect.7–10 Perhaps the biggest limitation for achieving agreement between theory and experiment has been the lack of an entirely satisfactory potential energy surface or surfaces (PESs).11–16 Here, we focus on the development of an accurate PES for the ground state based on explicitly correlated multireference configuration interaction (MRCI-F12)17, 18 data computed using MOLPRO,19 which can be used to make accurate predictions of the kinetics and dynamics. The electronic structure of ozone is very challenging to describe accurately with ab initio calculations. It has been noted as a problematic case for highly accurate singlereference-based schemes to calculate atomization energies.20 The best current experimental value for the dissociation energy from the latest active thermochemical tables (ATcT) is D0 = 8563 ± 3.5 cm−1 , which corresponds to De ∼ 9219 ± 10.0 cm−1 (we assign larger uncertainty to the zero-point energy, or ZPE of O3 ).21 As we have reported previously, De for standard MRCI varies with basis set from 7927 cm−1 a) Electronic addresses: [email protected] and [email protected]
0021-9606/2013/139(20)/201103/4/$30.00
at the AVTZ level to 8768 cm−1 at the AV6Z level, an increase of more than 800 cm−1 .15 This makes complete basis set (CBS) extrapolation uncertain and we previously noted a variation of ∼150 cm−1 between various CBS schemes.15 The large effect of the Davidson correction (>300 cm−1 ) on De is evidence of the importance of high-order correlation.22 Dynamic correlation is very important in general, since De at the full-valence complete active space self-consistent field (CASSCF) level is only 3003 and 3104 cm−1 for the VTZ-F12 and VQZ-F12 basis sets,23 respectively, and the shape of the minimum energy path (MEP) towards dissociation changes qualitatively with the inclusion of dynamic correlation. Watts and Bartlett24 found that for single reference coupled cluster, the effect of adding triple excitations, contributed changes in the two harmonic stretching frequencies ω1 and ω3 of more than 100 cm−1 , a rather dramatic effect of high-order correlation. For ozone, we find that at geometries away from the global minimum the system becomes increasingly multiconfigurational. This puts single-reference methods out of the question for a global PES making a multireference method such as MRCI the appropriate choice. There are also reasons to consider multiple electronic states. Adiabatic dissociation of ozone 1 3 O3 A1 → O 3 g− +O Pg , leads to a threefold degeneracy of singlet states (2 1 A + 1 A in the Cs group). Diabatically, the ground state of ozone correlates to excited states of both the atomic and molecular oxygen products 1 O3 (1 A1 ) → O2 1 g + O Dg ,
1
yielding another block of 10 degenerate singlet states (5 1 A + 5 1 A in the Cs group). 139, 201103-1
© 2013 AIP Publishing LLC
201103-2
Dawes et al.
A critical issue for the PES of ozone in the bottleneck region is an avoided crossing with an excited state that can result in disruptions in the ground state surface seen as a submerged reef and an associated van der Waals (vdW) well in the exit channel. The electronic structure in the bottleneck region was the main subject of our previous study with many details provided in the supplementary material of that paper.15 To briefly summarize, we found that the height of the submerged reef and depth of the vdW well are substantially reduced by both basis set and active space completeness. We followed the switch in character between the ground state and an excited state correlated with O2 (1 g ) + O(1 Dg ) products and noted that if all 13 singlet states were included using a dynamically weighted state-averaged method25, 26 (DW-SACASSCF), then the transition was better behaved and the reef all but disappeared. Finally, by including more than one reference state in the internally contracted icMRCI (thus increasing the number of contracted configurations and reducing the effect of internal contraction), we obtained a MEP that is a monotonic function of the distance to the dissociating O-atom (no reef feature). Given the sensitivity of the kinetics (noted in our previous study) to the shape of the PES in the bottleneck region, our goal has been to construct a full-dimensional global PES useful to explore the unresolved issues surrounding the kinetics and dynamics of this important system. We identified four features desired in a PES suitable for dynamics calculations: (1) excellent equilibrium structural parameters, (2) good agreement with experimental vibrational levels, (3) accurate dissociation energy, and (4) most-notably, a transition region with accurate topography (especially without a spurious reef). Until recently, the best MRCI-based PES for ozone was that initially developed by Siebert, Schinke, and Bittererova (SSB) in 2002.11, 12 The SSB PES is nominally based on 5000 spline-fit points at the MRCI/VQZ level starting from a one-state CASSCF calculation with the 2s orbitals closed in the reference space (12,9). As discussed above, direct use of MRCI at that basis level substantially underestimates the dissociation energy, so an adjusted hybrid PES was constructed.13 Schinke et al.4 were also aware of the sensitivity of the kinetics to the height of the pronounced reef found in their ab initio data and some numerical experiments were reported using versions of the PES in which the reef was removed. The low-lying vibrational levels on the SSB PES are very realistic matching experiment to a root-mean-squared error (RMSE) of ∼ 4 cm−1 for levels up to 4000 cm−1 above the ZPE.12 Straightforward application of the MRCI method at this basis set level does not perform so well and the recent (2013) PES of Ayouz and Babikov27 employing the same (12,9) reference space, but extrapolating the MRCI data to the CBS limit (MRCI/CBS) overshoots the fundamentals ν 1 and ν 3 by roughly 20 and 30 cm−1 , respectively, with larger errors recorded for overtones and combinations. Those results closely match our unpublished tests in which we have produced accurate fits at the MRCI/CBS level. The new Babikov PES includes the full symmetry and has an accurate dissociation energy of 9177 cm−1 (close to the experimental value of 9219 ± 10.0 cm−1 discussed above). However, it lacks a spin-orbit (SO) correction which would
J. Chem. Phys. 139, 201103 (2013)
lower the asymptote by ∼80 cm−1 (pushing De further from experiment) and influence the shape of the long-range PES. It also has a pronounced reef feature which means that overall it meets two of our four above-stated criteria for our new PES. Very recently, Tyuterev et al.28 published a spectroscopic PES fit to an analytic form describing one of the three global minimum isomer wells. This PES uses one-state Davidson corrected MRCI data starting with a full-valence (18,12) CASSCF reference. The PES combines data generated with the AV5Z basis (in the well), with other data extrapolated to CBS(5,6) (along the dissociation coordinate) thus achieving a dissociation energy of De = 9164 cm−1 (similar to that obtained by Babikov27 ). The authors report the submerged reef feature in their data, but produced two versions of their PES including one for which a correction was introduced based on our previous report (thus removing the reef). Interestingly, the presence or absence of the reef was found to influence the high-lying vibrational levels. Ultimately, they arrived at a PES without a reef feature that matches the experimental vibrational levels including high-lying levels up to ∼7900 cm−1 above the ZPE to within about 1 cm−1 RMSE.28 The highest lying experimental levels around 7900 cm−1 are only about 650 cm−1 away from the experimental value for D0 of 8563 ± 3.5 cm−1 .21 The Tyuterev PES undershoots D0 by only ∼55 cm−1 , but adding a correction for spin-orbit (not yet included) would affect the shape of the bottleneck and asymptotic region and bring D0 on the PES to ∼135 cm−1 below experiment. Thus, if higher-lying experimental levels were known, the error in calculated levels would have to increase rapidly from sub-wavenumber up to 7900 cm−1 above ZPE, to more than 100 cm−1 over the next ∼500 cm−1 toward dissociation. The Tyuterev PES is already proving very useful for understanding the spectroscopy of the known levels including determining the states involved in particular resonance polyads and helping to determine the terms and parameters in effective spectroscopic Hamiltonians.29 All of which aid in the development of line-lists for atmospheric databases such as HITRAN.30 Thus, the Tyuterev PES meets at least 3 of our 4 criteria desired for our planned studies of kinetics and dynamics, and will likely be useful for future comparisons. In view of our previous study and the observed sensitivity of the shape of the MEP and dissociation energy to basis set completeness and correlation treatment, we tested the recently developed MRCI-F12 method along the MEP. Explicitly correlated F12 methods have proven to greatly improve convergence with basis set, and been shown to produce near CBS quality results with relatively small basis sets.31 Figure S1 in the supplementary material32 plots a 1D cut using MRCI-F12 showing that the pronounced reef seen at the triple-zeta level for standard MRCI is nearly absent in the MRCI-F12 data, and at the quadruple-zeta level a monotonic cut is obtained (a result only achieved at the CBS level for standard MRCI). The dissociation energy De increases from 9162 cm−1 at the VTZF12 level to 9354 cm−1 with VQZ-F12 (without spin-orbit or other corrections). This is a much smaller step than found in standard MRCI for which De increases by nearly 700 cm−1 between AVTZ and AVQZ. The value for De of 9354 cm−1 obtained with VQZ-F12 already slightly exceeds experiment. Adding a correction for spin-orbit (which we do) lowers De by
201103-3
Dawes et al.
FIG. 1. Contour plot of the asymptotic region with rOO fixed at 2.282 bohrs, with R = [3.5, 12] bohrs and θ = [0◦ , 90◦ ]. The contours are in 10 cm−1 steps, the one labeled 9270 cm−1 would be the last to appear at any distance since De = 9275.12 cm−1 on the fitted PES with SO-correction.
∼80 cm−1 and adding a further scalar relativistic correction (which we neglect), lowers De by another 21 cm−1 ,15 reaching De = 9253 cm−1 (just above the estimated experimental value of 9219 ± 10.0 cm−1 discussed above). Feller and Peterson33 have reported systematic overestimation of atomization energies in some systems by F12 methods when extrapolated to the CBS limit consistent with our results. They also conclude however, that unextrapolated F12 data at the VQZF12 basis level “provides a very cost effective alternative to achieving results comparable to much larger standard method calculations.” Thus, we fit a global PES using MRCI-F12/VQZ-F12 data without CBS extrapolation (some additional details are given in the supplementary material).32 Our automated PES generation method described previously34–36 was used to fit the surface out to R(O–OO) = 10 bohrs in Jacobi coordinates including 2663 ab initio data (some obtained by permutation). No energy-dependent bias was used for the point selection, so the four wells (including the high-energy ring-minimum) and the asymptotic region were treated equally. To include the long-range part of the PES, we switch after 8 bohrs to the 2D (R, θ ) long-range interaction model of O2 + O by Lepers et al.37 To describe the r-dependence of the O2 fragment asymptotically, we switch simultaneously to an accurate potential energy curve for O2 by Ruedenberg et al.38 We also added a correction for SO coupling (some additional details are given in the supplementary material).32 The long-range model of Lepers et al.37 has a similar optional SOcorrection, so conveniently a consistent switch could be made to the long-range, with or without the SO-correction. No reef or vdW well is observed with or without the SO-correction. Figure 1 shows a contour plot of the PES with SO-correction in the asymptotic region. The C2V global minimum on the fitted PES has geometric parameters of re = 2.4031 a.u. [r0 (expt) = 2.4052 a.u.], and θ e = 116.84◦ [θ 0 (expt) = 116.75◦ ]39, 40 while the D3h ring minimum is found at an optimized structure of re = 2.715 a.u. (a value very similar to previous reports),12 with an energy of 10 755.1 cm−1 (1580 cm−1 above De , see the supplementary material32 ). We calculated the low-lying vibrational levels (up to 4000 cm−1 above ZPE) for the 16 O3 and 18 O3 mass combi-
J. Chem. Phys. 139, 201103 (2013)
FIG. 2. Cut along the bond stretching coordinate in valence coordinates comparing our new MRCI-F12 PES with the Babikov 2013 PES.27 The other bond is fixed at r1 = 2.282 a.u. and the angle is fixed at θ = 116.75◦ .
nations localized in one well, using the parallel ro-vibrational DVR code RTR3D by Wang and Carrington.41 Reasonably accurate vibrational levels are produced (RMSE for 16 O3 and 18 O3 are 14.4 and 7.8 cm−1 , respectively) that do not degrade substantially toward higher energy. They are not as accurate as the spectroscopic PES of Tyuterev,28 but indicate good representation of the wells in a PES intended for studies of dynamics. Tyuterev et al.28 report some influence of the reef (or its removal) even on levels 600-700 cm−1 below dissociation. Grebenshchikov et al.42 previously studied the highest lying states on a PES with the reef and reported many states localized in the vdW wells, and many others exhibiting mixed or delocalized character. In the supplementary material,32 we compare all of the bound states for the 16 O17 O18 O mass combination using the Babikov 2013 PES27 (with reef) and ours (without reef). Figure 2 shows a simple unrelaxed cut along the bond stretching coordinate in valence coordinates comparing our new MRCI-F12 PES with the Babikov 2013 PES. The pronounced reef feature in the Babikov PES changes the high-lying vibrational levels qualitatively (see text in the supplementary material),32 while for our new PES, in addition to supporting more bound states, even the high-lying states are mostly recognizable from one-well calculations. We applied a three-dimensional (3D) quantum statistic model (QSM) to predict the temperature dependent rate of the O + O2 exchange reaction, extending our previous twodimensional (2D) work. Figure 3 compares 3D results using the new PES to our previous 2D results and experiment. As seen in Figure 3, the rate predicted by the new 3D calculation is slightly lower than the 2D result. It underestimates the experiment at low temperatures. However, the negative temperature dependence qualitatively reproduced in the new 3D results is a direct result of the reef-less MEP. It is known that the QSM is a crude approximation to the dynamics and not strictly applicable in this case because of the known nonstatistical behaviors of this reaction.6, 43 Nevertheless, we take these results as an indication that more sophisticated stateto-state scattering calculations might be expected to display the desired negative temperature dependence on this PES, and hopefully produce much better agreement overall. In addition, a time-dependent wavepacket calculation for J = 0 indicates a rich resonance structure shown in Figure S6 of the
201103-4
Dawes et al.
J. Chem. Phys. 139, 201103 (2013) 5 S.
Y. Lin and H. Guo, J. Phys. Chem. A 110, 5305 (2006). Sun, L. Liu, S. Y. Lin, R. Schinke, H. Guo, and D. H. Zhang, Proc. Natl. Acad. Sci. U.S.A. 107, 555 (2010). 7 M. H. Thiemens, Science 293, 226 (2001). 8 M. H. Thiemens, Annu. Rev. Earth Planet. Sci. 34, 217 (2006). 9 K. Mauersberger, D. Krankowsky, C. Janssen, and R. Schinke, Adv. At., Mol., Opt. Phys. 50, 1 (2005). 10 Y. Q. Gao and R. A. Marcus, Science 293, 259 (2001). 11 R. Siebert, R. Schinke, and M. Bittererova, Phys. Chem. Chem. Phys. 3, 1795 (2001). 12 R. Siebert, P. Fleurat-Lessard, R. Schinke, M. Bittererova, and S. C. J. Farantos, J. Chem. Phys. 116, 9749 (2002). 13 D. Babikov, B. K. Kendrick, R. B. Walker, R. T. Pack, P. Fleurat-Lessard, and R. J. Schinke, J. Chem. Phys. 118, 6298 (2003). 14 F. Holka, P. G. Szalay, T. Müller, and Vl. G. Tyuterev, J. Phys. Chem. A 114, 9927 (2010). 15 R. Dawes, P. Lolur, J. Ma, and H. Guo, J. Chem. Phys. 135, 081102 (2011). 16 D. Xie, H. Guo, and K. A. Peterson, J. Chem. Phys. 112, 8378 (2000). 17 T. Shiozaki, G. Knizia, and H.-J. Werner, J. Chem. Phys. 134, 034113 (2011). 18 T. Shiozaki and H.-J. Werner, Mol. Phys. 111, 607 (2013). 19 H.-J. Werner, P. J. Knowles, G. Knizia et al., MOLPRO , version 2012.1, a package of ab initio programs, 2012, see http://www.molpro.net. 20 D. Feller, K. A. Peterson, and D. A. Dixon, J. Chem. Phys. 129, 204105 (2008). 21 B. Ruscic, private communication of ATcT result based on version 1.110 of the Core (Argonne) Thermochemical Network (17 May, 2010). 22 Y. G. Khait, W. Jiang, and M. R. Hoffmann, Chem. Phys. Lett. 493, 1 (2010). 23 K. A. Peterson, T. B. Adler, and H.-J. Werner, J. Chem. Phys. 128, 084102 (2008). 24 J. D. Watts and R. J. Bartlett, J. Chem. Phys. 108, 2511 (1998). 25 R. Dawes, A. W. Jasper, C. Tao, C. Richmond, C. Mukarakate, S. H. Kable, and S. A. Reid, J. Phys. Chem. Lett. 1, 641 (2010). 26 M. P. Deskevich, D. J. Nesbitt, and H.-J. Werner, J. Chem. Phys. 120, 7281 (2004). 27 M. Ayouz and D. Babikov, J. Chem. Phys. 138, 164311 (2013). 28 V. G. Tyuterev, R. V. Kochanov, S. A. Tashkun, F. Holka, and P. G. Szalay, J. Chem. Phys. 139, 134307 (2013). 29 A. Barbe, S. Mikhailenko, E. Starikova, M.-R. DeBacker, V. G. Tyuterev, D. Mondelain, S. Kassi, A. Campargue, C. Janssen, S. Tashkun, R. Kochanov, R. Gamache, and J. Orphal, J. Quant. Spectrosc. Radiat. Transf. 130, 172 (2013). 30 L. S. Rothman, I. E. Gordon, Y. Babikov, A. Barbe, D. C. Benner, P. F. Bernath, M. Birk, L. Bizzocchi, V. Boudon, L. R. Brown et al., J. Quant. Spectrosc. Radiat. Transf. 130, 4 (2013). 31 H.-J. Werner, T. B. Adler, G. Knizia, and F. R. Manby, in Recent Progress in Coupled-Cluster Methods, edited by P. Cársky, J. Paldus, and J. Pittner (Springer, 2010). 32 See supplementary material at http://dx.doi.org/10.1063/1.4837175 for more details of electronic structure calculations, vibrational levels and QSM and scattering calculations. 33 D. Feller and K. A. Peterson, J. Chem. Phys. 139, 084110 (2013). 34 R. Dawes, X.-G. Wang, A. W. Jasper, and T. Carrington, Jr., J. Chem. Phys. 133, 134304 (2010). 35 R. Dawes, D. L. Thompson, A. F. Wagner, and M. Minkoff, J. Phys. Chem. A 113, 4709 (2009). 36 R. Dawes, D. L. Thompson, A. F. Wagner, and M. Minkoff, J. Chem. Phys. 128, 084107 (2008). 37 M. Lepers, B. Bussery-Honvault, and O. Dulieu, J. Chem. Phys. 137, 234305 (2012). 38 L. Bytautas, N. Matsunaga, and K. Ruedenberg, J. Chem. Phys. 132, 074307 (2010). 39 Vl. G. Tyuterev, S. Tashkun, P. Jensen, and A. Barbe, J. Mol. Spectrosc. 198, 57 (1999). 40 Vl. G. Tyuterev, S. Tashkun, D. W. Schwenke, P. Jensen, T. Cours, A. Barbe, and M. Jacon, Chem. Phys. Lett. 316, 271 (2000). 41 RTR is a package of programs to compute Rovibrational levels and wavefunctions of TRiatomic molecules, X.-G. Wang and T. Carrington, Jr. 42 S. Y. Grebenshchikov, R. Schinke, P. Fleurat-Lessard, and M. Joyeux, J. Chem. Phys. 119, 6512 (2003). 43 A. L. Van Wyngarden, K. A. Mar, K. A. Boering, J. J. Lin, Y. T. Lee, S.-Y. Lin, H. Guo, and G. Lendvay, J. Am. Chem. Soc. 129, 2866 (2007). 6 Z.
FIG. 3. Comparison of 3D QSM calculation for rate of exchange reaction (O + O2 ), shown in black solid line, with previous 2D result (red curve) and experiment (dashed curve).
supplementary material32 which compares results for three PESs including our new one. To summarize, we report a new PES for the ground state of ozone based on explicitly correlated MRCI-F12 data with four key features: (1) excellent equilibrium structural parameters, (2) good agreement with experimental vibrational levels, (3) accurate dissociation energy, and (4) most-notably, a transition region with accurate topography (especially without a spurious reef). D0 = 8610.9 cm−1 on the new PES exceed the best experimental value of 8563.6 cm−1 by just 47.3 cm−1 , which could be improved another 21 cm−1 by including a so far neglected scalar relativistic correction. The PES includes a spin-orbit correction and switches to a longrange electrostatic model defined out to infinity. The absence of the submerged barrier has important implications for kinetics. As we have demonstrated with an approximate quantum scattering model, the new PES brings the calculated rate constant for the O + O2 exchange reaction into much better agreement with experiment and produces the negative temperature dependence observed experimentally. The new 3D PES will be useful for quantitative characterization of the O + O2 exchange reaction and ozone formation in full dimensionality. This work was supported by (U.S.) Department of Energy (DOE) (DE-FG02-05ER15694 to H.G.) and National Science Foundation (NSF) (CHE-1300945 to R.D.). R.D. thanks Maxence Lepers, Béatrice Bussery-Honvault, and Olivier Dulieu for a Fortran coding of their long range electrostatic model and Vladimir Tyuterev and Dmitri Babikov for useful discussions. A Fortran coding of the PES is available upon request from R. Dawes ([email protected]). 1 S.
M. Anderson, F. S. Klein, and F. Kaufman, J. Chem. Phys. 83, 1648 (1985). 2 M. R. Wiegell, N. W. Larsen, T. Pedersen, and H. Egsgaard, Int. J. Chem. Kinet. 29, 745 (1997). 3 P. Fleurat-Lessard, S. Y. Grebenshchikov, R. Schinke, C. Janssen, and D. Krankowsky, J. Chem. Phys. 119, 4700 (2003). 4 P. Fleurat-Lessard, S. Y. Grebenshchikov, R. Siebert, R. Schinke, and N. Halberstadt, J. Chem. Phys. 118, 610 (2003).
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
