A benchmark study of dioxygen complexes based on coupled cluster and density functional theory

A set of five compounds containing peroxo, superoxo or bis-μ-oxo moieties has been studied in the gas phase using CCSD(T)/aug-cc-pVQZ, also in combination with Goodson’s continued fraction approach. The corresponding analytical frequencies corroborate assignments of bands from experiments, and thus provide a consistent set of reference data that can be used for benchmarking a range of density functional approximations. A total of 100 density functionals have been checked for the general bond lengths, the specific peroxo/superoxo bond lengths, angles, and vibrational frequencies. There is not one density functional that performs equally well for all of these properties, not even within one class of density functionals.


Introduction
Dioxygen can coordinate to first-row transition metals (in different oxidation and spin states) with several protonation states, and hapticity, leading to e.g.peroxo, superoxo, hydroperoxo or bis-μ-oxo species (see Scheme 1).Although crystal structures exist for all of these species [1,2], the corresponding complexes are too large for treatment by accurate wavefunction methods like CCSD(T) with large basis sets (aug-cc-pVQZ or better).This prevents the systematic validation [3][4][5][6][7][8] of more efficient yet less accurate quantum-chemistry methods, e.g. based on density functional theory [9], leading to uncertainty about the reliability and appropriateness of different density functionals [10,11] for this metal-dioxygen chemistry.Fortunately, in the early 2000s a series of papers by Schnöckel [12,13] and co-workers reported the detection of peroxo and bis-μ-oxo complexes of aluminium (see Scheme 2), with vibrational spectroscopic data available for characterization.Furthermore, Tremblay and Roy characterized silicon trioxide compounds spectroscopically [14].Since these complexes are small (maximum six atoms), and highly symmetric, it is possible to study them by high-level wavefunction methods with sufficiently large basis sets.Hence, here we used CCSD(T) [15] with the aug-cc-pVQZ (and aug-cc-pV(X+d)Z [16], X=T,Q) basis sets for complexes 1-5 (Scheme 2) to optimize their structure and compute the corresponding vibrational frequencies.Because of the well-known instability of coupled cluster frequencies for close-lying electronic states, in particular when using unrestricted Hartree-Fock (HF) orbitals [17], we have complemented these data with CCSD(T) results based on restricted open-shell Kohn-Sham (ROKS) PBE orbitals, and CASPT2 frequencies.These results match the standard CCSD(T) results very well, when the underlying HF solution exhibits no close-lying electronic states.These data can then serve as reference data for the development of new density functionals.Here, the results were used as reference data for the benchmarking of density functional approximations, to explore how well the latter reproduce the metal-O and O-O distances, the angles, and the vibrational frequencies.

Results and Discussion
Coupled cluster results.The experimentally observed vibrational data [12][13][14] for 1, 3, 5 are reported in Table 1, together with the CCSD(T)/aug-cc-pVQZ values.Also included are two related compounds, 2 and 4, for which only the coupled cluster data are available.For 1 and 2 we also report CASPT2 results based on a CASSCF (14,8) reference.In general, there is a very good agreement between the experimentally observed peaks and the coupled cluster data.For instance, for [FAlO 2 ] 0 the peaks were observed experimentally at 1077 and 782 cm -1 , which match the harmonic coupled cluster frequencies of 1102 and 797 cm -1 ; note that the difference between experiment and theory (2-3%) is due to anharmonicity.We have computed the anharmonic CCSD(T) frequencies based on VPT2, which shows a very good match of the effect of anharmonicity on these frequencies (see Table 2).The only exception perhaps may be the 1015 cm -1 band of 3, which after anharmonic corrections is observed at 995 cm -1 with CCSD(T).
Few of the computed frequencies were observed experimentally, either due to too low intensity or not being able to assign these specifically to one of the manifolds of possible compounds that might be present in the experiments.For instance, in the reaction of AlF with O 2 the authors observed (at least) three species: 1, 3 and a third one involving four oxygens (assigned [12] as a di-superoxo species [FAl(O 2 • -) 2 ] 0 ).Based on isotope effects when treating all species with either 16 O 2 , 18 O 2 or a mixture of these, they were able to assign most of the observed peaks to one of these three species.a) shown in (parentheses) are the computed IR intensities (km•mol -1 ); b) [12]; c) CCSD(T)/aug-cc-pVQZ; d) ROKS-CCSD(T)/augcc-pVQZ on PBE orbitals; e) [13]; f) [14].
The geometric variables (bonds, angles) from the coupled cluster calculations are reported in Table 3.Based on the revised covalent radii by Alvarez and co-workers from 2008 [18], one would expect Al-F bonds of 1.78 Å, Al-O bonds of 1.87 Å, and Si-O bonds of 1.77 Å.Here, we observe Al-F, Al-O and Si-O bonds that are significantly shorter (0.15-0.20 Å).However, they are consistent with previous computational studies at lower levels of theory [13,19,20].Furthermore, the coupled cluster method with such a large basis set is often considered the 'gold standard' of computational chemistry.Indeed, the strong coherence between the computed and observed IR frequencies reinforce this notion.For this reason, the coupled cluster data provide a good reference set to be able to compare density functional approximations with.Density functional results.A total of 100 density functional approximations (DFAs) were tested against the coupled cluster reference data.Note that the same conditions were used as were employed for coupled cluster: using the aug-cc-pVQZ basis in the gas phase and without relativistic corrections.Hence, we looked at the bond lengths and angles (see Table 3), and frequencies.
There were four DFAs (PWPB95, mPWB1K, mPW1B95, PW6B95) that showed SCF problems during one or more of the optimizations and frequency calculations, which in all cases involved the 2 [FAlO 2 ] known numerical instability of the B95 correlation functional [21], which was solved by increasing the integration grid to (770,123) for the singlet systems.Further modifications were needed for the unrestricted 2 [FAlO 2 ] •+ system: changing the SCF procedure to EDIIS and loosening the convergence criteria for energy (10 -7 Hartree) and density (10 -5 atomic units) solved the issues there as well.
Most of the density functionals are able to predict well the angles (Table S3)  The picture changes again completely when looking at the vibrational frequencies (Table S4, Figure 2).Note that the DFT frequencies were computed at the CCSD(T) geometries; this allows to separate clearly the effect of accuracy for geometries (Tables S1-S3) from the accuracy for frequencies (Table S4).Here the hybrid functionals in all SciPost Chem.3, 001 (2024) the different guises are observed at the upper half of the list, with the O3LYP hybrid functional performing best (MAD 8.1 cm -1 ), followed closely by B97 (8.4 cm -1 ), and B97-1 (8.7 cm -1 ).The best composite is again r2SCAN-3c (9.4 cm -1 ).The best MGGAhybrid functional is now TPSSh (MAD 9.5 cm -1 ), TMHF is the best local hybrid again (MAD 10.1 cm -1 ), M06-L is the best metagga functional (MAD 10.2 cm -1 ), RC04 the best LDA functional (10.3 cm -1 ), HSE06 the best range-separated hybrid functional (MAD 11.1 cm -1 ), OPBE the best GGA functional (MAD 11.5 cm -1 ), wLH23tde the best rangeseparated local hybrid (MAD 13.1 cm -1 ), and PWPB95 the best double-hybrid functional (MAD 21.2 cm -1 ).Most of the list of 100 DFAs show MAD values for the vibrational frequencies that remain below 20 cm -1 , and 68 are within 10 cm -1 from the best performing one.
The average PROD value for all 100 functionals is rather high (0.32 Å 2 •°•cm -1 ), but this is distorted by the high value for HF-3c (2.477 Å 2 •°•cm -1 ).Instead, by taking the product of the averages for each of the four properties (Table 4), one obtains a PROD value of only 0.025 Å 2 •°•cm -1 , which is a more appropriate overall measure.Balance of exchange and correlation parts.For the dioxygen bonds, there is not a clear trend on the importance of either exchange or correlation.For instance, when comparing the deviations for O 2 distances by a variety of functionals based on the PBEc and Becke88x functionals, we observe a clear and large effect of the exchange part.Whereas PBE itself slightly underestimates the O 2 bond lengths (as estimated by inspecting the MD value in Table 5), the revPBE and RPBE functionals largely improve upon this.These latter two DFAs modify the exchange part of PBE.Likewise, OPBE, S12g and PBE1PBE, which modify the PBEx in different ways, show a systematic underestimation of the O 2 bonds by 0.030 Å (S12g) to 0.066 Å (PBE1PBE).This appears to indicate a major effect by the choice of exchange functional.However, by comparing the results of combining Becke88 exchange with different correlation functionals, similar effects are seen (ranging from -0.006 to +0.023 Å).Hence, both parts of the DFA formulation play an important role.

Conclusions
A new set of small dioxygen compounds is presented, for which high-level ab initio calculations were feasible at the CCSD(T)/aug-cc-pVQZ level.The resulting data could be used as reference in the development of new density functionals, or used to benchmark existing density functionals.The latter is done here, for a total of 100 density functionals coming from ten classes (LDA, GGA, MGGA, hybrid, MGGA hybrid, range-separated hybrid, local hybrid, range-separated local hybrid, double hybrid and composites).In general, the density functionals are remarkably accurate for bonds and angles, with average deviations of 0.027 Å for all bonds, 0.77° for angles; dioxygen bonds are more difficult to get right, with an average deviation from the coupled cluster of 0.058 Å, more than twice as large as the deviation for all bonds.The cause for this is difficult to pinpoint onto ingredients of DFAs or specific chemical bonding patterns.The frequencies from coupled cluster calculations are off by ca.10-20 cm -1 , where it should be noted that the well-known Hartree-Fock instability can significantly affect the frequency values for open-shell systems.
Basis sets not available in standard CFOUR (e.g.def2-TZVPD or aug-cc-pV(Q+d)Z) were obtained from the basis set exchange [30].The frozen core approach was used within the solving of the coupled cluster equations.The geometry optimization typically had SCF and coupled-cluster convergence criteria of 1.0•10 -10 , geometry converge of 1.0•10 -8 atomic units, leading to electronic energies accurate to at least 1.0•10 -8 Hartree.
The two-electron integrals were in all cases treated with the AOBASIS formalism to reduce disk space, in conjunction with the ECC program (with conventional UCCSD(T) for complex 2).
Additional CCSD(T) calculations using the same basis set and R(O)KS orbitals obtained by the PBE functional were performed with the MOLPRO program package (Version 2022.2) [31,32].All calculations employed the spin-unrestricted CCSD(T) routines (UCCSD(T)) to correctly handle the KS orbitals.
At the same level, also CASPT2 calculations were performed using a CASSCF wave function based on an active space containing the 2s and 2p orbitals of both O atoms.The initial guess for the active space was obtained using the atomic valance active space (AVAS) procedure [33] based on R(O)HF orbitals and both CASSCF and CASPT2 employed density fitting [34].

Figure 2 .
Figure 2. Mean absolute deviations (cm -1 ) of vibrational frequencies for the 100 density functionals compared to the CCSD(T)/aug-cc-pVQZ reference data.Shown from left to right are the classes of composites, LDA, GGA, MGGA, hybrid, MGGA-hybrid, rangeseparated hybrid, local hybrid, range-separated local hybrid and double hybrid.

Table 2 :
Effect of anharmonicity on computed vibrational data (cm -1 ) for dioxygen complexes a .

Table 3
a) not used for the deviations of density functional results (use of symmetry makes this value dependent on the other angle)

Table 5 .
Statistics a (Å) for deviations for O2 bonds by DFAs based on PBEc and Becke88 parts.