Investigation of the ozone formation reaction pathway: Comparisons of full configuration interaction quantum Monte Carlo and fixed-node diffusion Monte Carlo with contracted and uncontracted MRCI

Ozone plays several crucial roles in the atmosphere,¹ including protecting life from harmful UV radiation, as well as participating in reactions with a number of trace gases.^2–4 Signatures of ozone’s usual isotope dynamics are imparted into other species and can provide insight into modern hydrological cycles and the dynamics of stratosphere-troposphere exchange, and through biogeochemical pathways they can also tell us about Earth’s environment millions of years ago.⁵ 

The measured thermal rate coefficients of ozone isotope exchange reactions

where x, y, and z represent the different ¹⁶O, ¹⁷O, and ¹⁸O masses, respectively, have steep negative temperature dependencies indicative of a barrierless mechanism to form the O₃^* complex.^6–9 Strong kinetic isotope effects (KIEs) have been observed in these nearly thermoneutral exchange processes, leading to mass-independent fractionation (MIF) in the stratosphere.^10–13 The dynamics are non-statistical and the relative roles of differences in zero-point energy, symmetry and nuclear spin-statistics, unbound resonances, competing stabilization processes, as well as nonadiabatic effects such as spin-orbit and derivative coupling, and geometric phase effects are still under investigation. Spectroscopic evidence also supports a barrierless topography of the potential energy surface (PES). Recent measurements of progressions of vibrational levels approaching the highest-lying bound states, combined with theoretical analysis by Tyuterev et al., are much more consistent with a barrierless PES.¹⁴ 

Efforts to construct a PES for ozone useful to investigate some of the above-mentioned dynamical processes are constrained and guided by the large number of recorded vibrational levels for various isotopologues as well as by the dissociation energy which is specified to remarkable precision (D₀ = 8563 ± 3.5 cm⁻¹) by the Argonne active thermochemical tables approach (ATcT),¹⁵ which also infers D_e ∼ 9219 ± 10 cm⁻¹. The long range interaction between well-separated O₂ and O-atom fragments has been characterized in an electrostatic approach by Lepers et al.¹⁶ The region that is least well-defined is the transition region just below dissociation. No vibrational levels have so far been recorded within 600 cm⁻¹ of dissociation, yet the topography in the uppermost part of the wells crucially determines the dynamics and kinetics, including the observed negative temperature dependence of the exchange reactions.

An accurate global PES for the ground electronic state of ozone was published by Dawes et al. in 2013.¹⁷ The PES lacks the spurious reef feature found in several previous studies^18–20 and was used in wavepacket based quantum scattering calculations to successfully reproduce the negative temperature dependence of exchange rate coefficients as well as the large KIEs.^21–24 More recently, the PES was used in studies of the total number of bound vibrational states, their symmetry, and their density as a function of energy.²⁵ Prior to constructing the PES, an understanding of the origin of the reef feature and its sensitivity to calculational parameters was sought.

Preliminary insight into the origin of the spurious reef was reported in a 2011 study that showed that the height and even the existence of the barrier depend on the details of the ab initio calculations.²⁶ The reef was attributed to a widely avoided (large energy gap) crossing in the transition region between the lowest lying excited ¹A′ state and the ground molecular state connected to the $O2(Σg−3)+O(P3)$ asymptote. If adiabatic dissociation of ozone is considered,

three-fold degeneracy is reached asymptotically, represented as (2A′ + A′′) in the C_s symmetry group. Diabatically, the ground state of ozone connects to excited states of both atomic and molecular oxygen

which combine asymptotically to give (5A′ + 5A″). These two lowest singlet blocks combine to give a total of 13 singlet states (7A′, 6A″). In that study, to facilitate the switch in state character and to represent the asymptotically degenerate states on an equal footing, the 13 lowest singlet states were included in state-averaged complete active space self-consistent field (SA-CASSCF) calculations with dynamic weighting (DW).²⁷ (In this paper, we will refer to single-state and state-averaged CASSCF as SS-CASSCF and SA-CASSCF, respectively.) Subsequent internally contracted multireference configuration interaction (icMRCI) calculations were found to be sensitive to several additional factors: (1) active space, (2) basis set completeness, (3) Davidson correction, (4) and perhaps the internal contraction error. The height of the reef feature was found to be lower using a full-valence active space (18e, 12o) than for a reduced space (12e, 9o) in which the 2s orbitals are held doubly occupied. The height of the reef is progressively lower for each basis set in the correlation consistent series aug-cc-pVnZ (n = 3-6) approaching the CBS limit. The Davidson correction yields a significantly more attractive PES with a lower barrier and a dissociation energy in much better agreement with the experiment, confirming the importance of high-order dynamic correlation. However, despite all of that, results at the icMRCI(Q)/CBS level (using a full-valence active space) still exhibit a very slight reef of about 10 cm⁻¹ in height. A 2013 PES reported by Ayouz and Babikov²⁸ fit to data at the icMRCI(Q)/CBS level, used a SS-CASSCF reference state with the reduced (12e, 9o) active space, and has a reef feature with a height of more than 100 cm⁻¹. The last important factor appeared to be the internal contraction error in the MRCI calculations. Internal contraction has since been shown to introduce small but kinetically significant errors in the bottleneck region of the PES for other systems.^29,30 In the 2011 study by Dawes,²⁶ it was found that by using two reference states in the icMRCI calculation, the transition was made slightly more attractive, just enough to be monotonically attractive (barrierless). The number of employed reference states affects but does not systematically control the internal contraction error. Uncontracted MRCI calculations are much more computationally expensive and so no explicit tests were performed at that time. Nevertheless it was concluded that by considering the four factors listed above, a more realistic barrierless PES could be obtained.

In 2013, Tyuterev et al.³¹ published a spectroscopic PES fit to an analytic form describing one of the three global minimum isomer wells. This PES used a full valence (18e, 12o) SS-CASSCF reference to calculate the one-state Davidson-corrected MRCI data. The PES combines data from the aug-cc-pV5Z (AV5V) basis with other data extrapolated to the CBS(5,6) limit. The authors reported a submerged reef in their data but produced two versions of the PES, one of which included a Dawes correction to remove the reef.^17,31 They found that the reef feature or its absence affected the highest-lying vibrational levels. The PES without a reef produced much better agreement with experimental level progressions.³² 

The 2013 PES by Dawes et al.¹⁷ (denoted DLLJG) was constructed using the insights reported in 2011. To promote convergence with respect to basis set completeness, the newly available explicitly correlated multireference configuration interaction (icMRCI-F12)³³ method was used. Explicitly correlated F12 methods have proven to greatly improve convergence with respect to the basis set and have been shown to provide near CBS quality with relatively small basis sets.³⁴ Using the VTZ-F12 basis³⁵ directly produced a tiny barrier (∼3 cm⁻¹), while bases at or beyond VQZ-F12 yielded barrierless results. The effect of the avoided crossing is reduced to a shoulder along the monotonically attractive cut.¹⁷ For the DLLJG PES, the VQZ-F12 basis³⁵ was used without further extrapolation as this best matched the dissociation energy. The full-valence (18e, 12o) active space was used, and to promote orbital stability for some stretched geometries, 20 singlet states (rather than 13) were included in the dynamically weighted DW-SA-CASSCF procedure. Further tests were conducted with respect to the number of reference states used in the icMRCI-F12 calculations (which affects the internal contraction error as mentioned above). Ultimately, 7 reference states were used (corresponding to all of the ¹A′ states from the two lowest blocks), making the PES still more attractive in the transition region. The DLLJG PES was used in time-dependent wave packet scattering calculations of the thermal rate constants for the O + O₂ isotope exchange reactions, which agree well with the experiment, including their negative temperature dependence.^21–24 We use the PES as a reference in comparisons of various computational methods, while recognizing that cancellation of errors clearly contributes to the accuracy of the PES.

Here, to more systematically determine the effect of internal contraction on the minimum energy reaction path (MEP), large basis set uncontracted MRCI (ucMRCI) and ucMR averaged quadratic coupled cluster (ucMR-AQCC)³⁶ calculations were performed with the COLUMBUS³⁷ program and compared with internally contracted MRCI (icMRCI) calculations performed with the MOLPRO³⁸ program. In the case of the uc calculations, the interacting space restriction³⁹ was applied except where noted differently. In an effort to assess the high-order correlation contribution and role of size-extensivity corrections such as the Davidson^40,41 and Pople⁴² corrections, denoted (Q) and (P), respectively, initiator full configuration interaction quantum Monte Carlo (i-FCIQMC) calculations were performed with the NECI⁴³ code, and density matrix renormalization group (DMRG) calculations were performed using the BLOCK code⁴⁴ interfaced to MOLPRO. Finally, fixed-node diffusion Monte Carlo (DMC) calculations were performed with the quantum Monte Carlo (QMC) package QMCPACK⁴⁵ to test the routine feasibility of benchmarking challenging reaction paths using that method. Section II briefly describes the Monte Carlo based methods, while Sec. III provides results and discussion, followed by a conclusion in Sec. IV.

Variational Monte Carlo (VMC) uses an approximate trial wave function $ΨT$⁠, an initial reference for QMC, to calculate the expectation value of the Hamiltonian, the integration of which is performed by a Monte Carlo method.⁴⁶ In this study, VMC is primarily used to optimize parameters of the wave function for subsequent use in the more accurate diffusion Monte Carlo (DMC) method. DMC uses the importance-sampled imaginary-time Schrödinger equation to evolve an ensemble of electronic configurations toward the ground state.⁴⁷ Due to the fermion sign problem, DMC adopts the fixed-node approximation.⁴⁸ Since both VMC and DMC are Monte Carlo methods, statistical uncertainty decreases as $1/N$⁠, where $N$ is the number of samples.

For all of the QMC calculations, multi-determinant Slater-Jastrow (MD-SJ) trial wave functions were used, which can be written mathematically as

where $ΨT$ is the trial wave function, $eJ$ is the Jastrow factor, $ck$ are the determinant coefficients for the multi-determinant expansion describing static correlation, and $Dk↑$ and $Dk↓$ are the spin-up and spin-down Slater determinants, respectively. The multi-configurational trial wave functions were generated with the multi-configurational self-consistent field (MCSCF) and configuration interaction (CI) methods from GAMESS.⁴⁹ The aug-cc-pVQZ (AVQZ) basis set was used with the full valence (18e, 12o) active space in a SS-CASSCF calculation. Each trial wave function was combined with a three-body Jastrow factor containing electron-electron, electron-nucleus, and electron-electron-nucleus terms. An expansion order of 10 was used for the electron-electron and electron-nucleus terms, and an expansion order of 3 was used for the electron-electron-nucleus terms, resulting in a total of 82 optimized parameters. VMC calculations with energy minimization were used to simultaneously optimize the Jastrow factor parameters and the coefficients of the configuration state functions (CSFs). For the initial optimization, the default cutoff lengths were used with the initial Jastrow parameters set to zero. The Jastrow parameters and CSFs were optimized simultaneously for 10 cycles. The scheme of Ma et al.⁵⁰ was used to correct for the electron-nuclear cusps. Selecting the number of CSFs used in the trial wave functions involves balancing cost with accuracy. In this study, to limit the QMC cost, each trial wave function employed a fixed number of 750 CSFs with the largest squared coefficients. The DMC calculations were performed with a target population of 2880 and a time step of 0.0005 a.u. for all geometries. All electrons were correlated in the DMC calculations.

FCIQMC is a quantum Monte Carlo method⁵¹ designed to converge to the full configuration-interaction (FCI) energy, i.e., the exact solution to the Schrödinger equation for a given basis set. Thus in contrast to DMC, which was described in Subsection II A, FCIQMC results are directly comparable to the basis set dependent results obtained by standard electronic structure methods. The method simulates the imaginary-time Schrödinger equation of the interacting Hamiltonian based on the stochastic population dynamics of an evolving set of walkers which live and propagate in Slater determinant space. The method is able to converge to the FCI energy of a system once a sufficient system-dependent number of walkers are reached. Its initiator extension (i-FCIQMC) is designed to accelerate convergence to the FCI energy.⁵² Both FCIQMC and i-FCIQMC methods have been used to compute FCI energies in several benchmark studies.^53–55 

For this study, the FCIQMC and i-FCIQMC calculations were performed with the NECI code using the aug-cc-pVDZ (AVDZ) basis set.

For all of the calculations, the O-O (O₂) fragment bond distance was fixed at its asymptotic equilibrium value of 2.282 a.u. The bond angle varies negligibly along the MEP in valence coordinates and so was held constant at 116.8°. Calculations were carried out for the dissociating bond at a series of distances between 3.30 a.u. and 4.95 a.u. which cover the transition region. An optional correction for spin-orbit (SO) coupling reduces the dissociation energy of the reference DLLJG PES by ∼80 cm⁻¹ at the asymptote. Since the SO correction was not applied to the energies compared here, the uncorrected PES with D_e = 9355 cm⁻¹ was used in all comparisons. The long range region of the DLLJG PES is consistent with other PESs and also with electrostatic treatments,^17,16 so to focus on the transition region, the zero of energy in the comparisons presented here was set at 4.95 a.u. along the dissociation coordinate. (This is ∼243 cm⁻¹ below the asymptote on the non-SO-corrected DLLJG PES.) At each point, icMRCI energies were computed with MOLPRO, and ucMRCI and ucMR-AQCC energies were computed using the COLUMBUS code. For both the icMRCI and the ucMRCI methods, a one-state full valence (18e, 12o) SS-CASSCF reference wave function was used followed by a one-state MRCI calculation with the standard relaxed Davidson correction (Q). Additionally, SA-CASSCF using three states (two ¹A′ and one ¹A″) was used in ucMRCI and ucMR-AQCC calculations. In all cases, precise agreement was obtained between the two codes for the SS-CASSCF reference energy. Three basis sets such as aug-cc-pVDZ (AVDZ), aug-cc-pVTZ (AVTZ), and aug-cc-pVQZ (AVQZ) were used,^56,57 and both the icMRCI and the ucMRCI energies (with and without the Davidson correction) were extrapolated to the CBS limit (l⁻³ formula)⁵⁸ using the AVTZ and AVQZ bases.

For i-FCIQMC, the restricted Hartree-Fock (RHF) method at the AVDZ basis level was used as the reference wave function. The simulation was initialized with a single walker on the reference determinant D₀, and then the calculation proceeded with a shift value of zero, allowing for an initial exponential growth of walkers. The initial time step was chosen to be 0.000 14 a.u., which was allowed to be dynamically updated so that multiple walkers spawning from the same attempt would be rare. For the initiator method, a parameter n_a governing which determinants to include in the initiator space was selected as n_a = 3 (walkers). If n_a = 0, then all determinants would have been included in the initiator space. To reduce the memory requirements and CPU time due to a large number of additional spawned walkers, a cutoff κ of 0.01 was applied, meaning that the walkers with weights greater than κ would be left unchanged. For each point, the total number of walkers was initially grown to 8.0 × 10⁶, after which point, a series of calculations using semi-stochastic propagation⁵⁹ were performed, which contained 450 000 of the largest weighted determinants in the core space. The number of walkers was then doubled and subsequently followed by semi-stochastic calculations until the total number of walkers reached 2.56 × 10⁸. When the walkers reached this value, in order to reduce the CPU cost, the core space was reduced to 200 000 of the most dominant determinants. After equilibrating semi-stochastically, the number of walkers was then doubled again to 5.12 × 10⁸ for all geometries. To determine how further increases in the total number of walkers could affect the FCIQMC energy, for the geometry of 3.75 a.u., the number of walkers was expanded to 1.0 × 10⁹, then to 8.0 × 10⁹ and even a final test at 16.0 × 10⁹.

Complete results including tables of energies and numerous comparison plots are given in the supplementary material. Here we provide a few relevant comparisons. As plotted in Fig. S1 of the supplementary material, energies computed with small basis, one-state icMRCI/AVDZ, exhibit a pronounced barrier with a height of ∼556 cm⁻¹. [In this discussion, we refer to the barrier height (if any) relative to what would be a van der Waals minimum if the barrier was not spurious. Depending on the basis set and method, the barrier under discussion might be submerged with respect to the asymptote.] With the Davidson correction, the barrier height is reduced to ∼272 cm⁻¹. Coincidentally, uncontracted MRCI/AVDZ without Davidson correction has a similar barrier as icMRCI(Q) (with Davidson correction) with a barrier height of ∼250 cm⁻¹. However, Davidson corrected ucMRCI(Q) is already barrierless even at the AVDZ basis level. For the AVDZ basis only, FCIQMC calculations with the initiator extension were used in an effort to benchmark the FCI limit and hence provide insight into the accuracy of the various calculations and the Davidson correction. Due to their enormous cost (10s of thousands of CPU hours), all of the geometries except 3.75 a.u. were converged at 5.12 × 10⁸ walkers with nominal statistical uncertainties on the order of 40 cm⁻¹. The i-FCIQMC energies at those points are between the icMRCI(Q) and ucMRCI(Q) results, which might lead one to conclude that the Davidson correction slightly overcorrects in this case. However, it is noteworthy that the uncertainties given for the i-FCIQMC method in Table SI (supplementary material) should be interpreted cautiously. The derived uncertainties at intermediate stages of the calculations (smaller numbers of walkers) did not accurately reflect the range of where the energy might end up. In fact, in these results, the point at 3.75 a.u. is much better converged than the others. For that point, no further lowering of the energy was obtained beyond 8 × 10⁹ walkers, but drops much larger than the nominal uncertainty at earlier stages were noted upon each doubling of the walkers. For the point at 3.75 a.u., the i-FCIQMC energy is below that of all the other methods including the Davidson corrected uncontracted MRCI(Q) result. Thus, the best estimate of the FCI/AVDZ energy is below ucMRCI(Q), but unfortunately it was deemed too computationally expensive to converge all of the points to that same degree. The total CPU time to converge to 5.12 × 10⁸ walkers for the points at 3.60, 4.05, and 4.50 a.u. was ∼83 000 CPU-hrs, ∼128 000 CPU-hrs, and ∼103 000 CPU-hrs, respectively. In order to reach 8.0 × 10⁹ walkers at 3.75 a.u., the computational cost was >175 000 CPU-hrs. Note that these results are for the AVDZ basis set and of course capturing a large fraction of the dynamic correlation energy requires much larger basis sets.

The i-FCIQMC method was able to provide insight into the multireference character of the electronic structure of ozone along the studied pathway. Note that the details of specific configurations important to the bonding description of ozone have been discussed before by Ruedenberg and co-workers.^60,61 Here we confine our discussion to the excitation levels and the implications for various computational approaches. At 3.75 a.u. with 1.0 × 10⁹ walkers (Table I), the two most heavily weighted determinants were the HF reference and a configuration related by a quadruple excitation, both of which had nearly equal populations of ∼40 000 walkers. As given in Table I, the next seven leading determinants, in terms of decreasing weights, were related by a single excitation, a triple excitation, and several double excitations. Tables II and III list the excitation levels for the most relevant configurations for various geometries along the path. At r_OO = 3.75 a.u., as shown in Table III, considering the first 10 000 determinants, ∼1800 determinants were found to be quadruple excitations with respect to the HF reference, while ∼1700 were pentuple excitations and ∼900 were hextuple excitations. This highlights the challenge that ozone presents for single-reference based approaches. Indeed, it has been noted previously that the triples contribution shifts the harmonic frequencies by more than 100 cm⁻¹.⁶² At the geometries of 3.60, 4.05, and 4.95 a.u. with 5.12 × 10⁸ walkers, many important contributions from quadruply excited configurations were also found. The second most heavily weighted determinant at 4.95 a.u. was a quadruple excitation, while at 3.60 and 4.05 a.u., it was found to be the third most weighted determinant. Interestingly at 4.95 a.u., of the top nine most heavily weighted determinants, four were quadruple excitations, and in the first 1000 determinants, the number of pentuple excitations (∼160) was found to be greater than the number of quadruple excitations (∼130). These examples of important contributions from high excitation levels indicate why even triple-excitation single reference based methods will not suffice for ozone.

There has been some speculation about a significant strong/static correlation contribution from the 3s and 3p orbitals in this system. The cost to perform CASSCF calculations with larger than full-valence active spaces has been prohibitive for ozone even for small basis sets. Here to explore this issue, the DMRG method (as implemented in the BLOCK program)⁴⁴ was used to perform SS-CASSCF calculations also opening the 3s and 3p orbitals to construct an (18e, 24o) active space. [We confirmed that for the (18e, 12o) active space numerically identical energies are obtained with MOLPRO and BLOCK.] The calculation was performed for a dissociating bond distance of 3.90 a.u. using the AVDZ basis set. The progression of energies is rather interesting. The RHF/AVDZ energy of −224.2362 a.u. is lowered to −224.5089 a.u. for SS-CASSCF with the usual full-valence active space (18e, 12o). This large drop is due to the fact that no particular configuration is very dominant and many configurations contribute significantly. The leading squared coefficient is only 0.34, with the next few values being 0.22, 0.11, 0.10, 0.05, etc. However, the SS-CASSCF energy for the (18e, 24o) active space including the 3s and 3p orbitals only lowers the energy to −224.5357. Thus, the further drop in energy from opening the active space from (18e, 12o) to (18e, 24o) is less than 10% of that observed going from RHF to (18e, 12o). This is perhaps not as significant a drop as might be expected if the 3s and 3p orbitals contribute significantly to the bonding description. Note that the corresponding icMRCI and icMRCI(Q) values are −224.9050 and −224.9443 a.u., respectively [the ucMRCI and ucMRCI(Q) values are −224.9134 and −224.9535 a.u.]. Thus even for the small AVDZ basis set, more than half of the correlation energy comes from higher virtual orbitals beyond 3s and 3p. Overall this indicates that while strong correlation within the valence space is very important, strong correlation contributions from 3s and 3p orbitals are relatively insignificant. On the other hand, dynamic correlation is also very large, and based on the large size of the Davidson correction, so is the high-order contribution not directly captured in the MRCI(SD) treatment.

The icMRCI method yields a significant barrier at each basis set with a barrier height of ∼240 cm⁻¹ remaining at CBS. The icMRCI(Q) and ucMRCI methods follow similar curves to each other at each basis level from AVDZ towards the CBS limit, likely due to fortuitous error cancellation. After extrapolation toward CBS, icMRCI(Q) still has a slight barrier of ∼15 cm⁻¹, while with ucMRCI, the barrier shrinks from ∼250 cm⁻¹ at AVDZ to a height of ∼25 cm⁻¹ after extrapolation. Comparisons of the various methods are plotted for AVTZ and CBS in Figs. 1 and 3 (many more plots are provided in the supplementary material).

Comparing Figs. 1 and 2 (and also Figs. S5 and S6 in the supplementary material), the reef is seen to be strongly dependent upon the basis set. The icMRCI(Q) has a barrier height of ∼272 cm⁻¹ for AVDZ, while at AVTZ, it is reduced to ∼70 cm⁻¹.

The Davidson corrected ucMRCI(Q) does not produce a barrier at any basis set level. Thus, as noted previously, the use of internally contracted methods can certainly have a dynamically relevant impact. Similar to the reference PES, the ucMRCI(Q) has monotonically decreasing energies. However, the inclusion of the Davidson correction to the one-state ucMRCI(Q) calculation makes the method too attractive in this application when compared to the DLLJG PES. It was determined to explore this further, since the ucMRCI calculations are more flexible (hence more accurate) than internally contracted ones. Figure 3 plots results for ucMRCI and ucMR-AQCC comparing energies computed using a single-state SS-CASSCF reference vs a three-state state-averaged SA-CASSCF reference and the AVTZ basis. Basis set dependence and CBS limit values are given in Fig. S9 of the supplementary material for the ucMR-AQCC method. For each of the two orbital choices (SS and SA), the three methods containing size-extensivity contributions (Davidson, Pople, and MR-AQCC) lead to closely spaced curves. This is a very satisfactory result since the three methods approach the size-extensivity problem in different ways where especially the MR-AQCC method has the advantage of including the respective corrections in an internally consistent way whereas the Davidson corrections are added only a posteriori. It is noted that in all these cases, a qualitative change to a barrierless MEP is observed. The use of state-averaged orbitals has a rather significant impact, making the PES significantly less attractive through the transition region and thus closer to the reference DLLJG PES. However, this dependence on the orbital reference raises concern about the completeness of even these high-level calculations since when approaching the FCI limit, energies should become independent of the chosen orbital set.

Given the overly attractive behavior of the uncontracted MRCI calculations [especially ucMRCI(Q)] through the transition region, it was decided to benchmark the corresponding dissociation energy D_e to see if it is too large as might be anticipated. As discussed above, the best determination of D_e is 9219 ± 10 cm⁻¹.¹⁵ Thus, prior to correcting for spin-orbit effects, a highly accurate calculation should be in the neighborhood of D_e = 9299 cm⁻¹. As mentioned above, the DLLJG PES based on explicitly correlated MRCI (icMRCI-F12/VQZ-F12) slightly overshoots D_e by 56 cm⁻¹. The performance of standard internally contracted MRCI has been discussed previously.²⁶ Benefitting from the cancellation of errors, icMRCI(Q)/CBS produces estimates of D_e, that are quite accurate, ranging from only about 2–140 cm⁻¹ below the benchmark, depending which CBS extrapolation scheme is employed (D_e varies substantially with the basis set resulting in significant uncertainty).²⁶ Here we find that single state ucMRCI(Q)/SS overshoots D_e rather significantly at the CBS limit, following the progression of D_e = 8670 cm⁻¹, 9330 cm⁻¹, and 9810 cm⁻¹ for AVTZ, AVQZ, and CBS, respectively (more than 500 cm⁻¹ too large). However, for single state calculations without Davidson correction (ucMRCI/SS), D_e follows the progression of 8156 cm⁻¹, 8729 cm⁻¹, and 9147 cm⁻¹ for AVTZ, AVQZ, and CBS, respectively (ending up slightly below the 9299 cm⁻¹ benchmark). If a three-state state-averaged calculation is performed at the asymptote (which provides a balanced representation of the asymptotic electronic structure similar to what is obtained using the dynamic weighting procedure employed in the DLLJG PES), then D_e is slightly increased even further and at the CBS limit overshoots the benchmark even without the Davidson correction. The progressions are 8405 cm⁻¹, 9000 cm⁻¹, and 9434 cm⁻¹ for AVTZ, AVQZ, and CBS without Davidson correction and 8748 cm⁻¹, 9433 cm⁻¹, and 9933 cm⁻¹ for AVTZ, AVQZ, and CBS with Davidson correction. Thus it was determined that the uncontracted MRCI(Q) calculations do indeed overshoot D_e as implied by the overly attractive behavior in the transition region.

The all-electron fixed-node DMC results (Table IV) predict a barrier within the uncertainties, in contrast to the reference PES. However, in these tests, only 750 CSFs were included in the trial wave functions, and it is not known at this time how sensitive the topography of the PES is to this restriction. The energies are compared with all-electron (AE)-icMRCI and (AE)-icMRCI(Q) energies computed using the aug-cc-pwCVnZ (n = 3-5) basis sets and extrapolated to CBS using an optimized power-law. The trial wavefunction was constructed using SS-CASSCF/AVQZ. The DMC energies are well below the icMRCI/CBS values, but not as low as those of the icMRCI(Q)/CBS method. (In fact they are quite similar to the icMRCI(Q)/AVQZ results.) Note that core-correlation does not have a large effect on the PES (e.g., vibrational levels, MEP, or D₀).

As an aside, both the restricted Hartree-Fock (RHF) and the configuration interaction with single and double excitations (CISD) methods produce no submerged reef or barrier at any basis set level (in contrast to most of the MRCI methods).

Calculations were performed with internally contracted and uncontracted MRCI, i-FCIQMC, and fixed-node DMC along the association/dissociation MEP and were compared to a spectroscopically and dynamically accurate barrierless PES. Comparing icMRCI and ucMRCI, it was found that the internal contraction error indeed plays a significant role in producing the reef feature. One-state calculations with icMRCI, icMRCI(Q), and ucMRCI all produced a barrier in contrast to the reference DLLJG PES. The different ucMRCI approaches including size-extensivity contributions (Davidson and Pople corrections, MR-AQCC) all produced monotonically decreasing energies, but with respect to the PES, were too attractive, which would be inconsistent with the spectroscopy and dynamics. The use of a state-averaged CASSCF reference had a fairly significant effect on the uncontracted MRCI(Q) calculations, reducing the attractiveness through the transition region. Results obtained using a three-state SA-CASSCF reference were much closer to the reference PES, raising concern, though, about the sensitivity to the orbital reference. CASSCF calculations performed using the DMRG method with an active space expanded beyond full-valence to include the 3s and 3p orbitals, obtained a relatively small strong-correlation contribution from those orbitals.

The i-FCIQMC method was used to benchmark energies at the AVDZ basis set. However, due to its very high computational cost, only the geometry at 3.75 a.u. was fully converged. The best resulting energy at that point was found to be even lower than ucMRCI(Q), but the other points were not well enough converged to draw conclusions about the reef feature. The important configurations determined by the i-FCIQMC method reflect the multireference character of ozone, indicating important determinants of quadruple (and even hextuple) excitation levels from the single reference at geometries along the pathway.

Thus, while strong correlation within the valence space is very important, strong correlation contributions from 3s and 3p orbitals are relatively insignificant. On the other hand, dynamic correlation is also very large and based on the large size of the Davidson correction, so too is the high-order contribution not directly captured in the MRCI(SD) treatment. These results highlight the challenging nature of ozone’s electronic structure. The present paper shows a snapshot of the current situation concerning the accurate calculation of the asymptotic region of the surface and explores some of the factors contributing to a spurious entrance channel energy barrier. The most desirable and definitive level of calculation would be uncontracted MRCI with large basis sets and excitations beyond doubles, which at this time is still prohibitively expensive. It may be valuable to explore further (to the extent that it is affordable) the relative roles of active and virtual orbitals, perhaps through a series of steps in a restricted active space RAS-type approach, and approaching large basis sets.

In the development of new PESs, usually requiring thousands of data points, it is necessary to find an acceptable balance of accuracy and cost. For multireference calculations, internal contraction, size-extensivity corrections, and explicitly correlated methods control cost and facilitate convergence with respect to the basis set. This often means exploiting cancelation of errors, which are worth investigating prior to investing huge numbers of CPU-hrs in a PES.

This sort of insight into ozone and other systems emphasizes the kinetically and dynamically relevant effects that may arise from subtle aspects of the topography in key regions of the PES. Clearly, the notion of particular RMS error measures corresponding to “chemical accuracy” (e.g., 1 kcal/mol) is misleading in the context of even very small spurious barriers, or large relative errors, particularly for reactions at low temperature.

See supplementary material for tables of electronic energies for all of the employed methods and a variety of additional plots. For the i-FCIQMC method, tables of determinant weights and excitation levels are provided.

R.D. acknowledges support from the U.S. Department of Energy (Grant No. DE-SC0010616). F.B.C.M. wishes to thank the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under Process No. 307052/2016-8. F.B.C.M. and H.L. thank the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP)/Texas Tech University SPRINT program (Project No. 2015/50018-9) for travel support. H.L., R.F.K.S., and F.B.C.M. are grateful for computer time at the Vienna Scientific Cluster (VSC), Project No. 70376, and on the computer cluster array of the School of Pharmaceutical Science and Technology, Tianjin University. We would like to thank Ron Shepard, Stephen Klippenstein, and Larry Harding for useful discussions.