General embedded cluster protocol for accurate modeling of oxygen vacancies in metal-oxides

Metal-oxides are a class of material with wide applications in fuel¹ and solar cells,² high-k dielectrics,³ and the catalysis industry.^4,5 As the most prevalent defect in metal-oxides, controlling the concentration of O vacancies (Ovs) in these systems underpins much of the major advances to their applications. The dominant quantity determining the Ov concentration is the Ov formation energy, E_Ov.

It is pivotal that E_Ov can be determined accurately. The Ov concentration can change by several orders of magnitude with small (⁠$∼0.1$ eV) changes in the E_Ov at a given temperature,⁶ resulting in drastic changes in thermodynamic, electronic, and optical properties of the metal-oxide. For example, such changes in the Ov concentration can change an insulating oxide into a photocatalyst⁷ or metal.⁸ E_Ov is also a measure of the reducibility of a metal-oxide system, with (thermal) reduction being a vital step in the thermochemical cycles for H₂O and CO₂ splitting.^9,10 Similarly, E_Ov has also been shown to correlate well with key catalytic properties, such as the adsorption and bond activation energies of various molecules^11–13 on metal-oxide surfaces.

The reliable experimental determination of E_Ov is very challenging as it sensitively depends on many factors¹⁴ (e.g., Ov concentration, presence of dopants, and crystallite size). As such, E_Ov has been mostly computed through electronic structure modeling, particularly with density functional theory (DFT). However, E_Ov is highly sensitive to the approximation of the exchange–correlation (XC) functional in DFT, oftentimes leading to large disagreements in their predictions. For example, in the (110) rutile TiO₂ surface,¹⁵ E_Ov can vary by more than 1.5 eV—around 50% of the absolute E_Ov—with similar discrepancies observed in MgO¹⁶ and other metal-oxide systems.^15,17 Additionally, DFT with the generalized gradient approximation (GGA) XC functional has been shown to be inadequate for modeling Ovs in transition metal-oxides, such as rutile TiO₂, predicting that the unpaired electrons produced during Ov formation are delocalized,¹⁸ while hybrid functionals that incorporate exact exchange predict localized electrons on adjacent Ti sites.¹⁹ There is a clear need for high-accuracy and well-converged reference values of the E_Ov for these systems.

The (electrostatic) embedded cluster approach²⁰ offers the potential to efficiently apply accurate methods to study the Ov. It limits explicit quantum-mechanical calculations to only a finite-sized quantum cluster, with the electrostatic interactions from the rest of the system approximated by point charges. Over the past few decades, there have been numerous applications of this approach to the Ov in metal-oxides.²¹ These studies have ranged from applying the basic electrostatic embedding that has been described^22–24 to more involved setups that utilize polarizable environments^25–28 (of varying degrees of sophistication) or couples the quantum cluster to the environment self-consistently via the “perturbed cluster” approach.²⁹ Regardless of the approach taken, the outstanding challenge is finding a quantum cluster that is sufficiently small, or even the smallest, which allows for inexpensive modeling at the reference level of theory while minimizing finite-size errors to converge results to the bulk (i.e., infinite size) limit. This process is difficult because the quantum cluster can take any chemical formula (i.e., size), and for each chemical formula, there can be many possible shapes, with widely differing convergence behaviors.^30,31

To circumvent the complexities with searching the entire size and shape space, converged clusters are normally selected from a set of clusters generated from chosen design rubrics.^24,32 The most common rubrics are to keep quantum clusters stoichiometric and spherical.^16,33 Identifying suitable quantum clusters that follow these design principles requires significant time investment and manual input, paired with lots of trial-and-error. As such, only a handful of clusters are normally created, making it difficult to affirm the quality of a cluster and study multiple crystal systems simultaneously. Additionally, recent work^32,34 has shown that these rubrics lead to poor convergence with cluster size.

There is growing interest in approaches that can lessen the manual labor involved. These approaches range from building clusters using layers (based on unit cell multipoles³⁵ or coordination spheres³⁶) to using building blocks of either the unit cell^31,34,37,38 or fully coordinated ions.^37,39 While the quantum cluster series generated from these approaches have been shown to converge properties of an ideal crystal [e.g., nuclear magnetic resonance (NMR) constants,³⁴ bandgaps,³⁷ and optical spectra^39,40], their extension to the study of surfaces or point defects such as the Ov is not clear and they still suffer from important deficiencies. For example, the positioning of the building blocks in the building block approach still requires manual definition. With the layered approach, the number of atoms within each layer increases significantly compared to the previous, providing little granularity in the sizes sampled. This leads to large converged quantum clusters that are not conducive for high accuracy calculations, which typically exhibit steep system size scaling.

To summarize, a systematic and general approach for designing a quantum cluster set, which provides good granularity in the sizes and shapes being sampled while ensuring rapid convergence, is currently lacking. In this work, we propose a quantum cluster design protocol (named SKZCAM after the authors’ initials), which achieves these qualities for computing accurate Ov formation energies in metal-oxides. The core enabling development is to put the control of the shape and size of a cluster into a robust and flexible framework using the radial distribution function (RDF) to divide metal cations into “shells.” The O anion positions for each cluster arise naturally from the metal cation positions based on the criteria that all dangling bonds are removed from the metal cations. The result is a process that requires no manual intervention to generate clusters of high granularity for virtually any metal-oxide crystal system and surface termination while converging rapidly with size.

This protocol is used to obtain small converged clusters for studying the E_Ov in the bulk and common surface planes of rock salt MgO and rutile TiO₂, two technologically relevant metal-oxide systems with properties that sensitively depend on the Ov. These clusters (involving fewer than 650 correlated electrons for all systems) enable the accurate local natural orbital (LNO) variant of the coupled cluster method with single, double and perturbative triple excitations [CCSD(T)]^41–45 to be applied to compute E_Ov. We use our computed values to benchmark common DFT XC functionals, ranging from GGAs up to double-hybrids (DHs), and find that all studied XC functionals underestimate E_Ov, with this error decreasing, on average, along the rungs of Jacob’s ladder.

We describe the embedded cluster and supercell approaches used to generate the necessary structures to evaluate E_Ov for bulk rock salt MgO and its (001) surface and bulk rutile TiO₂ and its (110) surface in this section. This quantity is calculated as

where E[P-MO] and E[D-MO] are the total energies of the pristine and defected metal-oxide systems, respectively. The defected system is created from the pristine system by removing an O atom without further geometrical optimization. Both systems are treated in the closed-shell singlet state. E[O] is the total energy of an O atom in the (unrestricted) triplet state and E_bind is the molecular binding energy of an O₂ molecule, which has been computed for various levels of theory in Sec. S1 of the supplementary material. The final two terms sum up to (half) the energy of an O₂ molecule, thus defining E_Ov under O-rich conditions.

The definition given in Eq. (1) computes the unrelaxed O vacancy formation energy E_Ov. While the use of unrelaxed Ov structures precludes a direct comparison to experiment when relaxation effects are significant, it provides an upper bound to the relaxed E_Ov and a valid reference when comparing DFT values to high accuracy methods at the same geometry. In any case, relaxation effects have been found to be negligible for rock salt MgO.^46,47 Additionally, while this effect is significant for rutile TiO₂,⁴⁸ it is hard to quantify accurately due to a strong dependence on the spin-state (see Sec. II C) and chosen DFT XC functional (see Sec. S2 of the supplementary material), making it most appropriate to evaluate the unrelaxed E_Ov for the purposes of this work.

The electrostatic embedded cluster approach, illustrated in Fig. 1(a), used in this study features a quantum(-mechanical) cluster centered around the O vacancy. To model the long-range electrostatic potential from the rest of the material, this cluster is surrounded by a sphere and hemisphere of point charges of radius 30 and 40 Å for the bulk and surface, respectively; formal point charges have been placed at the metal and O ion crystallographic positions. In the vicinity of the quantum cluster (<7 Å), the positive point charges are “capped” with the effective core potential (ECP) of the corresponding metal ion, taken from the Stuttgart/Cologne group,^49,50 to avoid electron leakage from the bonded O ions at the boundary of the quantum cluster. The chosen radii of the different regions (see Sec. S3 of the supplementary material) can converge E_Ov to 0.01 eV. The placement of the point charges and ECPs alongside the atoms of the quantum cluster was constructed using py-ChemShell.⁵¹ 

DFT calculations were performed in ORCA³⁵ version 5.0 and MRCC⁵² 2020, with the latter interfaced to libXC.⁵³ The def2-SVP, def2-TZVPP, and def2-QZVPP Weigend–Ahlrichs⁵⁴ basis sets were used throughout this paper, with the standard def2-JK^55,56 fitting basis set used for Coulomb and exchange integrals. Convergence tests indicate that the def2-SVP and def2-TZVPP basis sets exhibit maximum errors of 0.34 and 0.02 eV with respect to the def2-QZVPP basis set (see Sec. S4 of the supplementary material).

Localized orbital correlated wave-function theory (cWFT) calculations were performed with the LNO-CCSD(T) and local Møller–Plesset perturbation theory (LMP2) implementations of Nagy et al.^44,57 in MRCC using the “normal” LNO thresholds. The MP2 contribution to the B2PLYP⁵⁸ XC functional was also evaluated using the LMP2 implementation of MRCC. Complete basis set (CBS) extrapolation parameters for the def2-TZVPP and def2-QZVPP pair, CBS(TZVPP/QZVPP), taken from the work of Neese and Valeev,⁵⁹ were used for the Hartree–Fock (HF) and correlation energy components of the cWFT total energy. Oxygen basis functions were placed at the Ov site, and only valence electrons were correlated. Convergence tests indicate that these settings can give accuracy to within 0.1 eV (see Sec. S4 B of the supplementary material). Deficiencies due to frozen-core or basis set size were accounted for through a correction computed on small tractable clusters featuring a “reduced frozen-core”⁶⁰ (e.g., Ne for Ti and He for Mg) and basis set involving cc-pwCVnZ^61,62 and aug-cc-pVnZ⁶³ on the metal and O ions, respectively (see Sec. S5 of the supplementary material), which has been extrapolated⁵⁹ for n = T and Q. The def2-nZVPP-RI auxiliary basis sets^64,65 were used for the local cWFT calculations with the def2 basis sets, while the automatic auxiliary basis functions of Stoychev et al.^66,67 were generated for the correlation-consistent (cc) basis sets.

Periodic supercell calculations with DFT were performed in the Vienna Ab initio Simulation Package (VASP).^68,69 These calculations serve to define the positions for constructing embedded cluster calculations and to provide reference E_Ov values. All systems used structures optimized at the R2SCAN⁷⁰ DFT level as these agree well with experimental lattice parameter values (see Sec. S6 of the supplementary material). The bulk rutile TiO₂ calculations employed a 192 atom $(22×22×4)$ supercell, and the bulk rock salt MgO calculations used a 64-atom (2 × 2 × 2) supercell, both with a (2 × 2 × 2) Γ-centered Monkhorst–Pack k-point sampling. The (001) MgO surface and (110) TiO₂ surface calculations both employed an asymmetric five-layer slab with the top two layers allowed to relax to form the pristine surface. A p(2 × 4) and (6 × 6) supercell was used for the TiO₂ and MgO surfaces, respectively, each computed with a (2 × 2 × 1) Monkhorst–Pack mesh and at least 12 Å of vacuum. For the TiO₂ surface, a correction to the p(2 × 6) supercell size, at the Perdew–Burke–Ernzerhof (PBE)⁷¹ level, was applied for the PBE0⁷² hybrid DFT calculation (see Sec. S7 of the supplementary material). An energy cutoff of 500 eV was used for all four systems, with small core projector augmented wave (PAW) potentials on the metal cations, leaving 12 and 10 valence electrons for Ti and Mg, respectively. The standard PAW potential was used on the O anion.

When constructing embedded cluster systems, the optimized bulk unit cells were repeated into supercells with dimensions larger than the embedded clusters to allow for them to be cleaved out. For the surface systems, the five-layer supercell slab was concatenated (along the surface normal direction) with more than 24 additional layers (taken from the bulk) before being repeated into supercells for embedded cluster calculations.

The unrelaxed Ov systems of rock salt MgO and rutile TiO₂, in both their bulk and common surface planes, all feature a well-controlled closed-shell singlet state.^26,46,47,73 While this spin state is appropriate for the MgO systems (since relaxation effects are negligible), there are significant discrepancies regarding the spin-state of the relaxed Ov structures of rutile TiO₂ within the literature, both from experimental and computational studies.

Experimental results from electron paramagnetic resonance (EPR) and infrared spectroscopy have detected that the excess electrons, believed to originate from Ovs, localize on Ti ions of the system, forming a deep bandgap state^74–76 and pointing toward an open-shell spin state.^77,78 However, these experiments also conflict with the high mobilities predicted by electrical measurements,⁷⁹ suggesting a shallow n-type donor,⁸⁰ which favors delocalized electrons.

Computational studies on the Ov of rutile TiO₂ have indicated various possible spin-states, ranging from the closed-shell singlet^73,81 to the open-shell triplet^26,82 and singlet^48,83,84 states. For the open-shell states, the formation of polarons⁸⁵ (i.e., localized electrons on Ti ions) brings additional complications. Depending on the localization site of the polarons and degree of polaronic distortion, there can be a wide range of formation energies.^19,84,86

A recent benchmark study using CCSD(T) by Chen et al.⁸¹ has found that the closed-shell singlet state is still the most stable state with the inclusion of relaxation effects. The validity of CCSD(T) for studying the Ov in TiO₂ was confirmed from preliminary full configuration interaction quantum Monte Carlo (FCIQMC) calculations, which reveal the single-reference character of the Ov state. However, it should be noted that Chen et al. did not directly model polarons due to the large cluster sizes required,³² which could potentially change the conclusions. We expect that the advances detailed in this work may provide the possibility to resolve these open questions/discrepancies for the TiO₂ system in future studies.

Designing a systematic and general set of quantum clusters for metal-oxides is a complex process; we have identified and summarized the key challenges involved in Fig. 1(b). Two highly important and interdependent factors are the size and shape of the quantum cluster. For metal-oxides, these two factors are predominantly controlled by the metal cations in the quantum cluster since these systems can be considered as formed from polyhedrons (e.g., TiO₆ octahedron for TiO₂) centered around the metal cations—highlighted in light blue for bulk rutile TiO₂ in Fig. 1(b). The choice of studied quantum cluster sizes is normally a manual process, often involving large arbitrary changes in sizes along the series. It is important that this process can be made systematic, while sampling clusters in small increments, so that small—even the smallest—converged clusters can be found for the property being studied, particularly when applying expensive cWFT methods. For a given size, the metal cations can take up virtually any spatial arrangement (i.e., shape) and E_Ov can vary significantly depending on the chosen shape, as shown in Sec. S8 A of the supplementary material. Under the electrostatic embedding scheme, the quantum cluster can be chosen to take up any charge. Thus, for a given metal cation configuration (i.e., size and shape), there can be virtually any number of O anions; if the ratio of O to Ti ions is more (less) than the stoichiometric ratio, then a negatively (positively) charged quantum cluster is formed. Here, the same challenges in deciding the number and spatial arrangement of O anions exist as the metal cations. If any of the factors discussed above are performed in an inconsistent manner, the convergence of the quantum cluster series becomes non-monotonic, as shown by the red line in Fig. 1(b). In this work, we propose the SKZCAM approach for constructing a set of quantum clusters, which shows fast and systematic convergence toward the bulk limit [blue line in the schematic graph of Fig. 1(b)]. We define a rigorous framework for deciding the metal cation configurations in the quantum cluster series, with the O anion configuration naturally arising based on a separate robust rubric.

The RDF of the number of metal cations as a function of distance from the Ov is used to generate the metal cation configurations in the quantum cluster series. Using bulk rutile TiO₂ as an example in Fig. 2(a), we see that the metal Ti cations arrange as shells (with the first three given distinct colors) of symmetry-related equidistant cations around the Ov, denoted by the gray sphere within a black circle. In an RDF plot, given at the bottom panel of Fig. 2(b), these Ti cation shells will appear as distinct peaks. Starting from the first metal cation shell/peak found in the RDF plot, metal cation configurations of systematically increasing size can be created by adding subsequent metal cation shells, with the first six metal cation configurations for this series visualized in Fig. 2(b). The RDF shells/peaks are completely controlled by the crystal structure, point defect site, and surface termination of the system, requiring no manual input. Furthermore, it provides good granularity in the size (i.e., number) of metal cations sampled in the quantum cluster series. At the same time, a variety of shapes are studied, ranging from spherical to cuboidal, all while ensuring symmetry about the point defect is maintained.

The metal cations serve as the base for the subsequent O anion configuration in each cluster. It is most physical to consider only O sites that ensure at least one bond to a nearby metal cation, shown as translucent red spheres for the quantum clusters in Fig. 2(b). In principle, any number of these O anion sites can be used; we take the unambiguous approach of placing O anions in all shown positions, ensuring no dangling bonds on the metal cations. This is the most chemically intuitive approach because it is equivalent to fully coordinating all of the metal cations for bulk systems. For surface systems, some metal cations at the surface are not fully coordinated due to the nature of the surface termination. As a result, the O anion configuration is completely determined by the metal cation configuration, which is, in turn, controlled by the RDF of metal cations around the point defect. We note that because the ratio of O anions to metal cations exceeds the stoichiometric ratio with this choice, the resulting quantum clusters are all negatively charged.

Beyond being fully systematic and general while providing good granularity, this SKZCAM approach also converges rapidly with cluster size. As shown by the blue markers in Fig. 2(c), E_Ov becomes converged to within 0.05 eV of the bulk limit (illustrated by the gray error bars) for the 22 Ti ion cluster (consisting of the first ten RDF shells/peaks), with all subsequent clusters staying converged. As a comparison, we have also constructed stoichiometric neutral clusters for a subset of the Ti cation configurations, as shown by the red markers. These clusters are the most common type within the literature, and we see a poor non-monotonic convergence toward the bulk limit, requiring clusters larger than the studied range of sizes. Notwithstanding the slow convergence, there is also significant ambiguity in constructing stoichiometric clusters for a given metal cation configuration since there are more possible O anion sites than allowed by stoichiometry. We give E_Ov for six possible O anion configurations in the three Ti ion quantum cluster (visualized in Sec. S8 B) in Fig. 2(c)—there is a wide range in E_Ov of 1.0 eV. For larger (stoichiometric) cluster sizes, there will be more possible O anion configurations; we have calculated E_Ov for only one such O anion configuration at larger cluster sizes based on including O anions closest to the Ov.

While the converged 22 Ti ion cluster found is already quite small, it may be beneficial to search for smaller converged clusters, particularly for performing expensive cWFT calculations. The SKZCAM approach provides a robust framework for defining the shape of a quantum cluster based on metal cation shells, which serve as the building blocks of the cluster. The 22 Ti ion cluster consists of the first ten metal cation RDF shells/peak, with the removal of the furthest (tenth) shell leading to an unconverged 20 Ti ion cluster. By systematically removing closer RDF shells from the 22 Ti ion cluster (see Sec. S9 of the supplementary material), we identify a smaller 18 Ti ion cluster formed from removing the seventh metal cation shell [see the green marker in Fig. 2(c)]. Removal of any subsequent shells from this 18 Ti ion cluster leads to large changes in E_Ov.

The SKZCAM approach outlined in Sec. III A is applied to create systematic sets of quantum clusters to study E_Ov at the DFT level and beyond for bulk rock salt MgO and its (001) surface, as well as bulk rutile TiO₂ and its (110) surface in Fig. 3.

For the MgO systems and TiO₂ bulk, there is systematic convergence toward the bulk limit, approximated from a supercell calculation. We consider convergence reached at the point when E_Ov starts to plateau, with the cluster located at this point being the smallest converged cluster (SCC). We find an SCC with 38 Mg ions and 17 Mg ions for the MgO bulk and surface systems respectively. Beyond these SCC sizes, E_Ov of larger clusters is less than 0.02 eV from the bulk limit, as illustrated by the gray error bars. Although the six Mg ion quantum cluster for bulk MgO is within 0.02 eV of the bulk limit, this agreement is fortuitous at the DFT level because at higher levels of theory (LMP2), it differs by > 0.3 eV from larger converged cluster sizes (see Sec. S11 A of the supplementary material). For TiO₂ bulk, the SCC with 18 Ti ions shows an error of 0.05 eV with respect to the bulk limit, with this error decreasing slowly with cluster size subsequently. This slow convergence could arise because the electrostatic embedding setup used here does not explicitly include (long-range) polarization effects; more sophisticated setups^87–89 including polarizable force fields beyond the quantum cluster could potentially improve this convergence.

The series of quantum clusters for the MgO systems were generated using the approach described in Sec. III of adding RDF shells of increasing distance from the Ov. We do not expect for there to be smaller converged clusters (in terms of the number of Mg cations) beyond the sampled series for these systems due to their high degree of ionic bonding and cubic symmetry. For rutile TiO₂, where there is a degree of directional covalent bonding⁹⁰ and anisotropy, there may be even smaller clusters. As discussed in Sec. III A, the SKZCAM approach provides the flexibility to find these smaller converged clusters, as has been done to find the 18 Ti ion SCC in TiO₂ bulk.

The (110) rutile TiO₂ surface requires a separate discussion due to its complex electronic structure. As seen in Sec. S12 of the supplementary material, the “noisy” nature of its convergence arises from well-behaved odd–even oscillations in E_Ov when Ti ions are added along specific crystallographic directions of the surface. Such odd–even size oscillations commonly appear in rutile TiO₂ surface^91–93 calculations. We find that the amplitude of these oscillations is correlated with the degree of self-interaction error (SIE) in the electronic structure method, being weaker in PBE0 compared to PBE, and completely absent in methods (e.g., HF and LMP2) that do not suffer from self-interaction error. The 21 Ti ion cluster was selected as the SCC on the basis of good agreement with the bulk limit at the PBE and PBE0 levels (<0.03 eV for the latter) and having reached the convergence plateau for HF and MP2 (see Sec. S11 B of the supplementary material). The gray 0.05 eV error bar in Fig. 3(d) indicates the average error of the clusters larger than (and including) the chosen SCC.

We find evidence that the SCC determined at an appropriate DFT level also leads to converged clusters—not necessarily the smallest possible—at other levels of theory, including cWFT methods. As seen in Fig. 4 for bulk rutile TiO₂, the SCC (predicted from PBE calculations) shows small changes that are less than 0.03 eV (indicated by the gray error bars) in E_Ov compared to larger quantum clusters at all studied levels of theory, from HF to PBE0 and LMP2. We observe the same behavior in the MgO systems and TiO₂ surface (see Sec. S11 of the supplementary material).

While there is good agreement between all levels of theory in the quantum clusters larger than (and including) the SCC from the DFT calculation, their behavior can vary significantly for smaller unconverged clusters. To bypass the steep scaling of cWFT methods, it is common to employ the $ΔLLHL$ approach within the literature^16,94,95 to produce reference quantities. Here, the difference between a high-level (HL) and a low-level (LL) theory (e.g., DFT) for a small tractable quantum cluster is added as a correction to the same low level theory computed at the bulk and basis set limits. The implicit underlying assumption is that the difference between the HL and LL methods stays the same regardless of cluster size. The differing size convergence of the various levels of theory suggests that such an assumption could exhibit uncontrolled errors if performed for arbitrarily small crystals.

Recent advances in localized orbital variants of CCSD(T) [e.g., LNO-CCSD(T),⁴⁵ DLPNO-CCSD(T),⁹⁶ and PNO-LCCSD(T)⁹⁷] have extended the remit of coupled cluster methods. To put these advances into perspective, for the def2-QZVPP basis set on a node equipped with 72 CPU cores, a PBE and PBE0 single-point calculation on the 17 Mg ion SCC for the MgO surface took 1.5 and 2.5 h, respectively, while this time increased to only 7 h for LNO-CCSD(T). Such a calculation would be far outside the reach of canonical CCSD(T), with the largest studied being a cluster with 6 Mg ions and 9 O anions.¹⁶ These developments, alongside the identification of relatively small and converged clusters, enable us to compute E_Ov at the local CCSD(T) level with large basis sets. The O vacancy formation energy, E_Ov, computed with LNO-CCSD(T) in the O-rich limit for the four systems is 7.68 ± 0.15 eV (MgO bulk), 7.18 ± 0.15 eV [MgO (001) surface], 6.39 ± 0.15 eV (TiO₂ bulk), and 5.55 ± 0.15 eV [TiO₂ (110) surface]. The decrease in E_Ov moving from bulk to surface can be expected due to the lowered coordination around the O anion sites on the surface.

Error bars have been added to the LNO-CCSD(T) values reported above. These have been conservatively estimated at 0.15 eV. This estimate comprises of errors arising due to the following: (i) basis set size and frozen core treatment, (ii) local approximation thresholds, and (iii) quantum cluster finite size errors. Tests on smaller clusters (see Sec. S4 B of the supplementary material) show that the (i) CBS(TZVPP/QZVPP) basis and frozen core treatment (Ar on Ti and Ne on Mg) chosen give good agreement (⁠$∼0.1$ eV) with respect to a larger basis set and small frozen core (Ne on Ti and He on Mg). Based on the deviations observed in small clusters, our best estimates of E_Ov were corrected for the bias due to the basis set and frozen core treatment (see Sec. S5 of the supplementary material). The (ii) LNO threshold settings chosen have been validated against canonical CCSD(T) to give small errors of <0.04 eV (see Sec. S13 of the supplementary material). Finite size errors are expected to be less than 0.05 eV—the typical error found for DFT with respect to the bulk limit (in Fig. 3) and from LMP2 calculations against larger clusters (see Sec. S11 of the supplementary material). Given the fast basis set convergence and all-electron nature of the embedded cluster DFT calculations, finite size errors (0.05 eV) are expected to be the dominant source of error.

On top of enabling the accurate LNO-CCSD(T) method to be applied, the embedded cluster approach virtually allows all DFT XC functionals to be applied at low cost, including double-hybrids—normally computationally inaccessible for solid-state periodic DFT codes. We use our obtained LNO-CCSD(T) reference values to evaluate the performance for a range of DFT XC functionals, with their deviation from LNO-CCSD(T) E_Ov values plotted in Fig. 5. These errors, alongside the computed E_Ov, for all of the studied methods are summarized in Table I. For all four systems, the studied XC functionals underestimate E_Ov with respect to LNO-CCSD(T). This error decreases, on average, when using XC functionals on higher rungs, with meta-GGAs, hybrids, and double-hybrids (DH) showing consistent improvement over the GGAs. The observed trends and variations of the XC functionals are quite consistent between the bulk and surface of the same system, but differ from one material to the next.

Out of the studied DFT XC functionals, the double-hybrid B2PLYP functional shows the best agreement to the LNO-CCSD(T) reference (mean absolute error, MAE, of 0.14 eV), with most points lying within the combined DFT and LNO-CCSD(T) error bars in Fig. 5. The ωB97X⁹⁸ hybrid and SCAN⁹⁹ meta-GGA XC functionals give the next best performance, both with an MAE of 0.23 eV. In the hybrid functionals, B3LYP¹⁰⁰ also shows good (MAE of 0.32 eV) and consistent performance across the four systems. On the other hand, PBE0 and HSE06¹⁰¹—two common hybrid functionals for metal-oxide systems—give variable performance, having large errors (of around 0.7 eV) for the MgO systems, which lowers (to around 0.2 eV) in the TiO₂ systems. The GGAs all severely underestimate E_Ov, with MAEs all exceeding 0.6 eV. The XC functionals in Fig. 5 are arranged based on the decreasing MAE within their respective Jacob’s ladder rungs. XC functionals that incorporate the Becke-88 (B) exchange and Lee–Yang–Parr (LYP) correlation functionals, such as BLYP and B3LYP, are one of the top performers in their respective rungs, with B2PLYP being the best overall studied XC functional, as discussed previously.

Lower-level cWFT methods, such as LMP2 and LNO-CCSD, are automatically generated in any LNO-CCSD(T) calculation, and they are also compared in Fig. 5. Both LMP2 and LNO-CCSD show excellent agreement to LNO-CCSD(T), with MAEs of 0.11 and 0.19 eV, respectively, and only the B2PLYP DFT XC functional has errors of a similar (small) size. As the errors of these three methods are all close to or within the error bars of the LNO-CCSD(T) values, it is not possible to ascertain which method performs better. We do find, however, that the good agreement of LMP2 appears to arise from fortuitous error cancellations, as elaborated in Sec. IV A.

The two key developments of this work are as follows: (i) a protocol for obtaining small converged clusters for performing reference calculations of the Ov at high levels of theory and (ii) the assessment of the accuracy of XC functionals for studying E_Ov. It is important that both these developments are properly contextualized, either from physical theory or comparison to the literature. For the latter development, we will try to rationalize the observed trends in performances of the XC functionals in Sec. IV A. For the former, we will compare our LNO-CCSD(T) reference values to those in the literature in Sec. IV B.

Seeking to understand the relative performance of DFT XC functionals in complex systems, such as metal-oxides, is not straightforward. In particular, E_Ov is a quantity that depends on many factors. Nonetheless, we believe that our results can reveal some useful insights on this property.

Systematic studies on a series of metal-oxides have shown that there is a correlation between the bandgap and E_Ov in these systems,^48,102 as is confirmed by the larger E_Ov in MgO (with a PBE bandgap of 4.83 eV in the bulk⁴⁷) over TiO₂ (1.88 eV PBE bandgap¹⁰³) in our results. Physically, this correlation arises because the removal of an O atom leaves two electrons (originally occupying the O 2p band), which must redistribute by occupying an empty band from the conduction band.¹⁰ This redistribution energy, and, in turn, E_Ov, would then be correlated with the size of the bandgap of the crystal system.

Within the same crystal system, we find that the predicted DFT bandgaps provide a rational basis for understanding the underestimation of E_Ov in most DFT XC functionals and the relative trends between the Jacob’s ladder rungs. More specifically, the errors in E_Ov with respect to reference methods should be related to the deviations of the XC functionals from the experimental or reference method bandgap. Due to the presence of self-interaction error¹⁰⁴ (SIE), semilocal functionals (meta-GGAs and below) will underestimate the bandgap, hence underestimating E_Ov, with meta-GGAs normally giving improved bandgaps over GGAs.¹⁰⁵ Hybrid functionals correct for some of the SIEs via the incorporation of exact exchange, with the possibility that bandgaps are overestimated.^99,106 For both the rock salt MgO and rutile TiO₂ systems, most of the studied hybrid functionals underestimate the bandgap,^37,107 hence underestimating E_Ov.

The bandgap can only serve as a general guide for the relative performance of DFT XC functionals, and there are functionals that do not follow this trend, suggesting that there are other important factors that influence E_Ov. The key example is the PBE0 functional in rutile TiO₂, which underestimates E_Ov despite predicting a bandgap¹⁰³ of 4.05 eV that overestimates experimental electronic bandgaps from photoemission experiments, typically in the range of 3.3–4.0 eV.^108–111 Additionally, despite predicting bandgaps that are 0.6 and 1.0 eV higher than HSE06 for bulk MgO^112,113 and rutile TiO₂,^73,103 respectively, these two functionals have very similar performances (both with MAE close to 0.5 eV) across the range of systems. Despite predicting worse bandgaps,¹¹⁴ the improved overall performance of the SCAN functional (with MAE of 0.23 eV) over many hybrid functionals is another example.

In GGAs, on top of the underestimated bandgap due to the SIE, a major contribution to its errors also arises from a poor description of the binding energy of the O₂ molecule—they tend to predict strong overbinding (e.g., 0.50 eV/atom for PBE). It is common within the literature to correct for this error by replacing the GGA O₂ binding energy in Eq. (1) with the experimental binding energy of 5.22 eV.¹¹⁵ With this change, the underestimation is decreased significantly, with MAE decreases of the range of 0.3–0.5 eV for the GGAs (see Sec. S14 of the supplementary material). XC functionals on higher rungs are less impacted, with improvements in the MAE of less than 0.2 eV for the meta-GGAs and hybrids. We also find that the better performance of LMP2 over (the more accurate) LNO-CCSD method is fortuitous, arising from its overbinding of 0.22 eV/atom (in Sec. S1 of the supplementary material), and it exhibits a larger MAE of 0.34 eV relative to the 0.11 eV of LNO-CCSD when using the experimental binding energy.

As the prototypical system for studying the Ov in metal-oxides, bulk MgO has been subject to several studies involving high level reference methods,^{16,29,47,116,117} and out of the four studied systems, it is the only system, to our knowledge, where E_Ov has been experimentally determined.¹¹⁸ Hence, bulk MgO makes for a good system to assess the accuracy of our obtained reference LNO-CCSD(T) values.

Experimental determination of E_Ov can be challenging due to the many experimental factors that can influence it. For MgO bulk, the single available value of 9.29 eV was obtained from additive coloring experiments by Kappers et al.^16,118 These experiments involve heating MgO crystals in Mg vapor under high temperatures and pressures. According to Smakula’s formula,¹¹⁹ the measured optical absorption spectra at these temperatures allow for the determination of E_Ov from the relative Mg vapor and Ov concentrations. As seen in Table II, this value differs by >1.5 eV from our LNO-CCSD(T) values and other high-level methods. This discrepancy is likely attributed to the uncertainties arising in the experiment. For example, the oscillator strength obtained by Kappers et al. differs significantly (>70%) compared to a previous experiment.¹²⁰ Richter et al.¹⁶ also attributed this discrepancy to thermal equilibrium not being reached in the experiment. Additionally, some of the assumptions in Smakula’s formula can be questionable for ionic solids, such as MgO.¹¹⁹ Beyond experimental uncertainties, some of this discrepancy can also arise due to the neglect of temperature effects in the static LNO-CCSD(T) calculations.

Given the above considerations, it is more appropriate to compare high-level references available for the E_Ov in MgO bulk. Ertekin and Grossman⁴⁷ evaluated E_Ov with DMC for MgO bulk, finding the DMC value to be 0.5 eV higher than PBE. These E_Ov values were computed under Mg-rich conditions, where E_Ov is defined as

with E[MgO] and E[Mg] being the total energy of bulk MgO and Mg per formula unit, respectively. Compared to the O-rich limit in Eq. (1), this definition will not suffer from the poor O₂ binding description in E_bind. If the O₂ overbinding (0.50 eV/atom) is corrected in the PBE E_Ov, our LNO-CCSD(T) value also becomes 0.52 eV higher than PBE, agreeing with DMC values. $ΔPBECCSD(T)$ embedded cluster calculations by Richter et al.¹⁶ have obtained an E_Ov value of 6.85 eV, which is 0.83 eV smaller than our LNO-CCSD(T) estimate. In Sec. S15 of the supplementary material, we show that 0.2 eV of this difference can be attributed to differing lattice parameters and inclusion of structural relaxation around the Ov by performing $ΔPBELNO-CCSD(T)$ calculations on the same quantum clusters as Richter et al. Remaining differences could arise from the use of differing embedding environments, which our work (Fig. 4 and Sec. S11 A of the supplementary material), alongside others,⁸¹ have shown cWFT methods to be highly sensitive to, requiring careful convergence. The reference E_Ov value from the $ΔPBECCSD(T)EOv$ study was also found to be smaller than the PBE (with overbinding corrected) value by 0.09 eV, which goes against trends observed from high-level calculations on other metal-oxide systems^121,122 and the DMC study of Ertekin et al. Explicit MP2 calculations on a converged cluster²⁹ give a E_Ov of 7.13 eV, closer to our LNO-CCSD(T) value. The difference with respect to our LNO-CCSD(T) or LMP2 values arises due to the use of a small (double-zeta quality) basis set; larger basis sets or CBS extrapolations can increase E_Ov significantly by >0.5 eV with respect to a double-zeta basis set for MgO bulk (see Sec. S4 of the supplementary material).

To our knowledge, the only available high level reference calculation in the TiO₂ system comes from a DFT+GW study,¹²³ with the large number of electrons in Ti making application of methods such as canonical CCSD(T) highly expensive.⁸¹ Additionally, DMC has not been applied to study the Ov in TiO₂ potentially because the inclusion of a 3d transition metal brings additional complications with the starting trial wave-function¹²⁴ and the need to validate its pseudopotentials.^125,126

In this work, we discuss a systematic and general approach (named SKZCAM after the authors’ initials) for designing small converged quantum clusters for studying the O vacancy with the electrostatic embedded cluster method. When combined with localized orbital correlated wave-function theory methods, such as LNO-CCSD(T), this approach allows for an accurate determination of the Ov formation energy (E_Ov) at a computationally tractable cost. We applied this approach to compute E_Ov values for the bulk and common surface planes of rock salt MgO and rutile TiO₂ systems. A comparison of these reference values to common DFT XC functionals shows that the studied XC functionals underestimate E_Ov for all studied systems. We observe general improvements in the XC functional errors as Jacob’s ladder is ascended, which we find can be correlated with the improvements in the predicted bandgaps of the XC functionals. Of the XC functionals studied, the double-hybrid B2PLYP functional gives the best performance, with a mean absolute error within the error bars of the reference calculation. Other BLYP-based functionals, such as BLYP and B3LYP, are also found to perform well within their respective Jacob’s ladder rungs, alongside the meta-GGA SCAN and the hybrid ωB97X functionals.

Although this work focuses on the Ov formation energy, the simple and intuitive nature of the outlined protocol should allow its application to other chemical problems, including molecular adsorption, spectroscopic quantities, or complex (ternary) metal-oxide systems. Furthermore, the approach can be automated, making it amenable for integration into the existing high-throughput calculation frameworks. Such high-throughput frameworks¹²⁷ for point defects or molecular adsorbates are highly desirable. For example, it can be used to screen for target applications in catalysis^11,102 or to produce large reference databases to validate current XC functionals.¹²⁸ Beyond the applications, we have defined a rigorous and well-founded framework for controlling the shape of the embedded quantum clusters. This framework can be valuable in not only defining the quantum cluster in electrostatic embedded cluster approaches but also in other embedding approaches,^20,129–132 where multiple clusters, corresponding to different levels of theory, may have to be defined simultaneously.

See the supplementary material for a detailed compilation of the obtained results and further data and analysis to support the points made throughout the text.

The authors acknowledge resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (Capital Grant No. EP/T022159/1) and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk); computational resources granted by the UCL Myriad and Kathleen High Performance Computing Facility (Myriad@UCL and Kathleen@UCL) and associated support service; computational support from the UK Materials and Molecular Modeling Hub, which is partially funded by EPSRC (Grant Nos. EP/P020194/1 and EP/T022213/1); and computational support from the UK national high performance computing service, ARCHER 2. The access for both the UK Materials and Molecular Modeling Hub and ARCHER 2 was obtained via the UKCP consortium and funded by EPSRC Grant Ref No. EP/P022561/1.

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article and its supplementary material. The input and output files associated with this study and all analysis can be found on GitHub at benshi97/Data_Embedded_Cluster_Protocol. The scripts used to generate the quantum clusters with the SKZCAM approach are also provided in a separate GitHub repository at benshi97/SKZCAM-Protocol and Binder. The data analysis and scripts can be viewed interactively online through associated Jupyter notebooks (via Binder), with links provided in the corresponding GitHub repositories.