Cerium Oxides without U: The Role of Many-Electron Correlation

Electron transfer with changing occupation in the 4f subshell poses a considerable challenge for quantitative predictions in quantum chemistry. Using the example of cerium oxide, we identify the main deficiencies of common parameter-dependent one-electron approaches, such as density functional theory (DFT) with a Hubbard correction, or hybrid functionals. As a response, we present the first benchmark of ab initio many-electron theory for electron transfer energies and lattice parameters under periodic boundary conditions. We show that the direct random phase approximation clearly outperforms all DFT variations. From this foundation, we, then, systematically improve even further. Periodic second-order Møller–Plesset perturbation theory meanwhile manages to recover standard hybrid functional values. Using these approaches to eliminate parameter bias allows for highly accurate benchmarks of strongly correlated materials, the reliable assessment of various density functionals, and functional fitting via machine-learning.

. The pseudopotentials employed are based on the frozen core approximation and the projector augmented wave method (PAW) [3], as specied in Table I. Those PAW potentials have been optimized for excited state properties (so-called GW POTCARs) including scalar relativistic eects. Furthermore, a norm-conserving (nc) PAW pseudopotential was used for Ce. Note that the pseudo-orbitals do not necessarily represent the charge of one electron, if non-norm-conserving PAW pseudopotentials are used, since the norm of the PAW pseudo-orbitals is not conserved in the core region. The calculation of expectation values thus requires so-called augmentation or compensation charges centered at the atoms. This correction is however missing for the unoccupied pseudo-orbitals, which leads to a poor description of unoccupied high-energy pseudoorbitals. In turn, errors are introduced in correlation energies which strongly depend on the high-energy orbitals. The importance of norm-conserving PAW pseudopotentials for correlation calculations was previously recognized for 3d metals in Ref. [27]. For s, p bonded materials, the eect is usually negligible. While the impact for lattice constants of ceria is also negligible, the reaction enthalpies change dramatically, see Tab   The k-mesh is specied by three numbers (k 1 × k 2 × k 3 ), corresponding to a uniform sampling of the Brillouin zone in each direction of the reciprocal lattice, including the center (Γ-point). Convergence of all calculated quantities with respect to the basis and k-mesh are discussed in sections V and VI.
Due to the presence of occupied f-orbitals, we include non-spherical one-center contributions to the gradient corrections (in GGA functionals like PBE or HSE03) inside the PAW spheres for all calculations. Additionally, the maximum angular momentum for augmenting the overlap densities in HF and correlation calculations is set to 6

II. GROUND STATES AND METASTABLE STATES OF SOLID Ce2O3
For open d-or f-shells materials, such as Ce 2 O 3 , self-consistent iterative algorithms (like conjugate gradient) are prone to getting stuck in local minima when applied to DFT+U , HF, or HSE functionals [9,22,32]. One physical reason is the near degeneracy of the multiplet structure of the f-orbitals. Another technical reason is the convexity of the Fock exchange functional and the Hubbard U term with respect to the partial orbital occupancies [32]. Note that the Fock exchange and the Hubbard U are responsible for a sizeable band gap, while LDA and PBE predict a metallic behavior with partial occupancies close to the Fermi level.
To address this issue we make use of the ramping method by Meredig et al. [32] in order to nd the correct semiconducting ground state at PBE + U and HSE03(α) level, based on the metallic ground state at the PBE level. In this method, the U or α parameter is gradually increased (starting from zero, i.e. PBE) while each self-consistent calculation starts with the orbitals and charge density from the previous calculation. In comparison to iterations starting from stochastic orbitals, we nd that the ramping method converges to the lowest ground states in a robust and reproducible way.

III. THE IMPACT OF THE ELECTRONIC STRUCTURE ON THE U DEPENDENCE
The Mott-Hubbard Hamiltonian is designed to capture localization in strongly-correlated materials stemming from on-site Coulomb interaction. The literature suggests however that the 4f band in Ce 2 O 3 is more complex in nature, presenting an unusually high degree of itinerancy. [12,16,17] Any such eects fall outside the on-site Hubbard correction and have to be captured by DFT. Alternatively, from the perspective of an (Anderson) impurity model itinerancy is best captured via hybridization with the surrounding Fermi liquid [13].
De Medici et al. has pointed out that both eects, i.e. localization and hybridization, are competitive in nature [8].
The authors demonstrated that at 0, K temperature the Mott transition is fully arrested due to the Kondo eect. In their Dynamic Mean-Field model any non-vanishing amount of hybridization suces. This impediment is tempered at higher temperatures, where given a suciently high U value, the Mott transition can re-establish itself.
We propose that an analogous scenario transpires with DFT+U, due to the vanishing electron temperature. In the case of pure cerium there has been a better appreciation for the dual nature of Ce(4f ) with an ongoing discussion into its impact on the γ-α phase transition [23]. The current consensus is that while the volume collapse may be explained by the Kondo eect at nite temperatures [18,25,42]

IV. THEORY BEHIND THE CORRELATION METHODS MP2, RPA, AND RPA+RSOX
We express the electron-electron interaction energy E ee in terms of mean-eld spin-orbitals, where the composite indices i, j, k, l, ... ∈ occ and a, b, c, d, ... ∈ unocc label the occupied and unoccupied mean-eld spin-orbitals, t ab ij are double excitation amplitudes and the electron-repulsion integrals read The well-known MP2 approximation to the correlation energy E c , Similar observations were made in Refs. [7,44]. corresponds to amplitudes dened by t ab ij = ab|ij /(ε i + ε j + ε a − ε b ), where Hartree-Fock spin-orbitals χ and energies ε are used. The RPA correlation energy is computed by means of the independent-electron response function P 0 (iν), with the Coulomb kernel r|v C |r = |r − r | −1 and As a basis Kohn-Sham spin-orbitals χ and energies ε from the HSE03 functional are used. In terms of the notation using amplitudes t ab ij , the RPA correlation energy corresponds to the neglect of the exchange-like term, i.e. t ab ij − t ab ji → t ab ij with amplitudes implicitly dened by where a sum over only k, l ∈ occ and c, d ∈ unocc is understood. This can be interpreted as a simplication of the amplitude equation for coupled-cluster singles and doubles, restricting to the direct particle-hole ring-like diagrams [38]. The solution of (6) can analytically be expressed as an innite sum over all ring-like Goldstone diagrams, The connection to Goldstone diagrams can be found in Figure 2 and Refs. [34,41]. The inclusion of −t ab ji is commonly denoted as RPA+SOSEX (RPA+second-order screened exchange), however, low-scaling algorithms are not yet available but possible [36]. Instead, we correct the missing exchange-like correlation in the RPA with a renormalized second-order approximation where f ij = i|f |j are matrix elements of the Fock matrix f in the Kohn-Sham spin-orbital basis. The denominator corresponds to a renormalization, where the ε's are the Kohn sham orbital energies. This renormalization is motivated by Epstein-Nesbet perturbation theory on top of the Kohn-Sham Hamiltonian [20]. It can be interpreted as an innite resummation of (here) second-order exchange-like diagrams, as depicted in Figure 4 in Ref. [21]. In fact, the bare SOX contribution (i.e. ∆ ab ji = 0) suers from the underestimated band gap of the Kohn-Sham basis. This would lead to a vast overcorrection of the exchange-like correlation, which is corrected by the renormalization ∆ ab ji in rSOX. The SOX correction to RPA accounts for the rst-order exchange-like contribution in the DFT basis. For the renormalized diagram, rSOX, we refer the reader to Figure 4 of Ref. [20].

V. BASIS SET CONVERGENCE AND THERMODYNAMIC LIMIT OF THE REACTION ENERGIES
Since we are using a periodic code, the reaction enthalpies primarily depend on two parameters, the basis set size (plane-wave cuto ENCUT) and the k-mesh to sample the Brillouin zone (BZ). Considering the high computational cost of correlation energy calculations, a deliberate and ecient convergence scheme has to be considered in order to reach the complete basis set and thermodynamic limits.
All presented enthalpies correspond to zero-temperature and neglect the pV term. Hence, the reaction enthalpies ∆H are simply given by Helmholtz free energy dierences, e.g. for the reduction 2CeO 2 Here F denotes the Helmholtz free energy at T = 0, while (s) and (g) indicate the bulk and gas phase, respectively. Since the pV term is neglected, (g) corresponds to an isolated molecule. In Sec. VIII we present a scheme to calculate ground state energies of isolated molecules using the periodic vasp code.
Since both, periodic (solids) and aperiodic (isolated molecules) systems are involved, we apply the following scheme in order to nd converged reactions energies.
1. Converge only the periodic part (F Ce2O3(s) −2F CeO2(s) ) with respect to the number of k-points using a reasonable basis set size (ENCUT = 631 eV). Results can be found in Tab. IV, V, and VI for the HSE03, EXX+RPAc, and HF+MP2c levels of theory, respectively.
2. Converge the reaction energy ∆H (i.e. periodic and molecular part) with respect to the basis set size employing the previously found k-mesh for the periodic part. For the molecular part, use a relatively small cell size (about L = 6.5 or 5.0 Å for mean eld or correlation energies, respectively) and a 2 × 2 × 2 k-mesh (see Sec. VIII for justication). Results can be found in Tab. X, IX, and XI for the HSE03, EXX/RPAc, and HF/MP2c level of theory, respectively.
3. Extrapolate the cell size of the molecule to the case of an isolated molecule (as described in Sec. VIII) using the previously found basis set size and correct the molecular part of the energy from the previous step. Results Tables X, IX, and XI, denoted as ∆H(L → ∞).

are incorporated in the
The correlation methods employ an additional auxiliary plane-wave cuto (ENCUTGW in vasp) for the calculation of two-electron repulsion integrals. It is rst restricted to ENCUTGW = 2/3 ENCUT. In a follow-up step, both RPAc and MP2c extrapolate the internal auxiliary basis set, ENCUTGW → ∞ [24,27,37]. This way the impact of ENCUT on the correlation energy is mainly left to unconverged mean-eld orbitals. Note that the internal cuto extrapolation is not used in the rst step. The absolute correlation energies of step 1 and 2 must thus not be compared.
The referred tables allow for a robust estimation of the numerical error due to the nite k-meshes and basis set sizes. The error of the box size extrapolation for obtaining isolated molecules is estimated via the quality of the t, see for instance Figure 4 and 7.

VIII. TREATMENT OF ISOLATED MOLECULES IN A PERIODIC CODE
For the reactions, we aim to nd the ground state energy (equivalent to the Helmholtz free energy at T = 0) of the isolated molecules H 2 O, H 2 , CO 2 , CO, and O 2 . In a periodic code, an isolated molecule can be simulated either by using a suciently large unit cell or by explicitly targeting the limiting case of an innitely dilute gas by means of predictive cell size extrapolations. In the rst case, periodicity is troublesome and the distance between periodic images has to be suciently large. In the second approach, periodicity is incorporated by interpreting the simulation as a gas of molecules (note, that gas is understood only as a technical term for periodically repeated molecules). Here, we chose the second approach for the following reason: the correlation energy varies inversely with the HOMO-LUMO gap of the mean-eld solution, while the HOMO-LUMO gap itself converges only slowly with the cell size. Hence, undesirable large unit cells have to be employed in order to achieve a converged correlation energy of an isolated molecule. Contrariwise, a predictive cell size extrapolation allows to restrict to smaller cells by employing tting functions, which cover the physical long-range interactions between molecules and the eect of the HOMO-LUMO gap.
In order to simulate a gas of molecules in a periodic code, the following technical aspects have to be considered. The number of interacting molecules is determined by the number of molecules in the Born-von Karman (BvK) cell. With computational eciency in mind, we increase this number by a uniform n k × n k × n k sampling of the Brillouine zone, being equivalent to a supercell containg N k = n 3 k molecules. The volume of the BvK cell is then given by Ω = N k Ω 0 , while Ω 0 is the volume of the unit cell containting only one molecule. The use of n k = 1 would prohibit any physical meaningful interaction between molecules, hence n k ≥ 2 is required. This can be seen by the construction of the electron-electron Coulomb interaction in reciprocal space: v(G) = 4π/G 2 with G = g + k where g is a reciprocal lattice vector and k is a sampling vector of the Brillouin zone (BZ). In real space the Coulomb kernel corresponds to the discrete Fourier series of v, v(r − r ) = 1 N k BZ k g v(g + k) e i(g+k)(r−r ) .
Equation (9) is equivalent to 1/|r − r | only in the thermodynamic limit, i.e. N k → ∞ and 1 The deviation is particularly substantial in the case of a Γ-only sampling of the BZ, i.e. N k = n k = 1 and thus The Coulomb potential in real space is then periodic along each lattice vector R and v(r −r ) takes only intramolecular interaction into account, while intermolecular interactions (with distances of the order of multiples of R) are mapped back to intramolecular interactions. In an intermediate regime, where the unit cell does not substantially exceed the spatial extend of the electron density of the molecule, the periodicity of v(r − r ) leads to a spurious mixing of intramolecular and intermolecular interaction. In this case, the BvK boundary conditions are simply unphysical.
In passing we note that this argument is not valid for extended molecules, where signicant intermolecular interactions can be of much smaller distance compared to the linear dimension of the BvK cell. A Γ-only approach can then still be sucient though. Since we are considering only small molecules, our approach of a predictive extrapolation of the intermolecular interactions requires a denser k-point sampling. This is particularly important to capture the dispersion interactions in correlation methods and will be discussed in the following subsections.
A further challenge arises from the g+k = 0 contribution in Eq. (9), which corresponds to the long-range interaction in real space, i.e. |r − r | → ∞. Although integrable in the thermodynamic limit, any nite k-mesh requires a special context-dependent treatment to avoid the numerical divergence. We discuss the solutions to this problem in the relevant subsections below.

HF, HSE, EXX
While the convergence of periodic DFT calculations of a molecule with respect to the linear unit cell dimension L = Ω 1/3 0 is rather cheap, the introduction of the Fock operator (as in HF or hybrid functional calculations) poses a particular challenge. This is due to the fact, that the Fock exchange operator suers from unphysical electrostatic interactions between BvK cells, when assuming periodic boundary conditions with a nite k-mesh (i.e. a nite volume of the BvK cell) [14]. Since the leading contributions of those interactions are of monopole-monopole character, the Fock exchange energy is divergent for any nite k-mesh. Formally, in a periodic code, the divergence of the exchange energy becomes apparent by the g + k = 0 contribution in Eq. (9). Note that the dierence between the classical electrostatic interaction (Hartree term) and the exchange interaction is that the former is largely compensated for by the positively charged nuclei, thereby preventing divergence.
Correcting the divergence from the exchange energy is thus mandatory. Among other known approaches [40], we consider two strategies to address this issue. (i) The rst is the probe-charge method (default in vasp) [33], which corrects for the spurious monopole-monopole interactions between BvK cells. This leaves us with non-diverging monopole-dipole interactions that are sill unphysical however. (ii) The second one is the Alavi-Spencer scheme [39], which applies a simple radial cuto to the Coulomb interaction in the Fock operator (HFRCUT=-1 in vasp). It essentially crops the interactions between BvK cells.
For range-separated functionals like the HSE03 hybrid functional, the situation is less demanding. By denition the Fock operator is restricted here to short-ranged interactions. These are enforced by multiplying the Coulomb interaction with the complementary error function, erfc(0.3 · |r − r |), which suppresses the interactions beyond 6.67 Å. The long-ranged exchange interactions are left to the quicker decaying PBE exchange functional. The exact exchange energy (EXX, i.e. the Hartree-Fock functional evaluated using HSE03 orbitals) is again exposed to the same divergence problem described above however and requires correction.
Technically, the physically meaningful range of the Coulomb interaction (10) is limited by the number of k-points in vasp, if the unit cell contains only one molecule. To accurately capture and extrapolate the electrostatic long-ranged electrostatic interactions, a k-point sampling beyond the Γ-only approach is necessary. This constitutes the rst step in our practical scheme for nding the limit of an innitely dilute gas (isolated molecule) at the HF, HSE03 and EXX level. (i) We increase the number of k-points until the ground state energy is converged to less than 1 meV, using a xed but reasonable large L (usually about 6 Å). This step ensures that the ground state energies are free from unphysical exchange interactions between BvK cells and that enough molecular images are incorporated (measured by the number of k-points) to accurately capture the intermolecular interaction energy per molecule. (ii) We extrapolate the ground state energy using E 0 + αL −3 + βL −6 , which accounts for intermolecular electrostatic monopole-dipole (L −3 ) and dipole-dipole interactions (L −6 ). We sample the data points over the range of L = 6 to 10 Å. The k-mesh has been determined in the previous step and remains xed throughout.
In these systems we use symmetry breaking orthorhombic unit cells, however, the parameter L = Ω 1/3 0 is still adopted as the linear cell dimension. Comparing approaches, the probe-charge method on the suers from the slowly decaying unphysical monopole-dipole interactions between BvK cells in the exchange energy. The Alavi-Spencer scheme on the other hand converges exponentially with respect to the number of k-points [39]. So even a 2 × 2 × 2 mesh is sucient for all considered molecules if L ≥ 6 Å, as illustrated in the left plot of Figure 3. We therefore regard the Alavi-Spencer scheme as superior for our setup. Indeed, the intermolecular HSE03 energies all nicely converge for a 2 × 2 × 2 k-mesh, as exemplied in the right plot in Figure 3. Similarly to HF, the EXX energy converges faster with the Alavi-Spencer scheme. Using this setting, the HF and HSE03 energy can then be extrapolated with an accuracy of 1 meV to the case of an isolated molecule (innitely dilute gas), see left panel of Figure 4.

RPA and MP2
A more subtle task is to perform a cell size extrapolation of the correlation energy with predictive power to the limit of an isolated molecule. In the following, we derive a tting function to extrapolate the correlation energy with respect to the linear unit cell dimension L, using the MP2 method as an example. The MP2 correlation energy of distant molecules with non-overlapping electron densities can be written as a sum over intramolecular and intermolecular pairs: The Coulomb integrals 2v ij ab v ab ij in the inter pairs account for the London dispersion interaction, L −6 . Meanwhile, the L dependency of the intra pairs is mainly driven by the change of the energetic distance between the occupied and unoccupied energies in the denominator, (ε i + ε j − ε a − ε b ) −1 , as indicated by the HOMO-LUMO gap in Figure  5. We thus assume a behavior following (1 + γL −3 ) −1 = 1 − γL −3 + O(L −6 ) for the intra pairs, since the gap closes largely with L −3 . Thus, for large L the slow convergence of the intramolecular correlation dominates over the more rapid convergence of the intermolecular correlation. A predictive extrapolation of the MP2 correlation energy should therefore cover both the closing HOMO-LUMO gap (L −3 ) and the dispersion interaction (L −6 ): The L −3 term here is not to be confused with the monopole-dipole interaction in Section VIII 1. Actually, the densities ϕ * i (r)ϕ a (r) are monopole free in perturbation theory, since the occupied and unoccupied manifolds are mutually orthogonal. In the particular case of MP2, the exchange part −v ij ab v ab ji is of very local nature and therefore contributes only to the intra pairs. Its opposite sign slightly mitigates the impact of the closing HOMO-LUMO gap, whereas RPA correlation energies lack exchange terms.
As discussed in the previous section, the range of the Coulomb interaction is technically limited by the number of k-points, if the unit cell contains only one molecule. Using a Γ-only sampling of the BZ, the correlation energy is thus expected to show only the L −3 behavior of the closing HOMO-LUMO gap. A 2 × 2 × 2 sampling of the BZ however can at least capture the dispersion interaction between nearest neighbors. Although the missing interaction of more distant neighbours decays only with N −1 k , the extrapolation (13) is remarkably stable with respect to the k-meshes employed, see left plot in Figure 6. A denser k-mesh captures more intermolecular dispersion per molecule, thus increasing the ratio |β/α| in Eq. (13). This leads to a further ranged dominance of the L −6 term. A similar eect can be observed when a nite size correction (FS) for the g + k = 0 contribution to the correlation energy is employed in Eq. (9). We apply FS using kp-perturbation theory [11], the output of which vasp stores into the WAVEDER le. Apparently, the missing long-ranged intermolecular dispersion is overcorrected again, leading to a further ranged dominance of the L −6 term. The failure of the Γ-only approach is evident from Figure 5. The HOMO-LUMO gap converges there under a constant slope, while ner k-meshes reach the nal slope earlier. Thus the extrapolation (13) predicts a too small energy gap in the denominator for Γ-only, which in turn yields too strong correlation.
We conclude that using a 2 × 2 × 2 mesh and Eq. (13) suces to extrapolate the MP2 correlation energy to the limit of an innitely dilute gas for all molecules under consideration. The correction of the g + k = 0 contribution (WAVEDER) does not lead to signicantly more accurate extrapolations. Since our MP2 implementation that automatically extrapolates the basis set does not support the g + k = 0 correction, we neglected this correction in favour of the basis set extrapolation. Finally, the cell size extrapolated results are supported by the MP2 calculations using a truncated Coulomb interaction ( Figure 6). In passing we note that truncating the Coulomb kernel in correlation methods is meaningful only if the HOMO-LUMO gap of the underlying mean eld solution is fully converged, hence very large L have to be considered.
Likewise, the RPA correlation energies are extrapolated using Eq. (13), since the arguments apply to all perturbative correlation energy methods. Figure 7    2 CeO 2 (s) + CO(g) −−→ Ce 2 O 3 (s) + CO 2 (g) and derive the experimental reaction enthalpies, ∆ r H 0 K , from available experimental atomization enthalpies, ∆ at H 0 K , While zero-temperature atomization energies are available [1] for the gases (the second value is the ZPE correction), we derive the atomization enthalpies for the bulk cerium oxides using Thus the atomization enthalpies of both oxides are given by which nally allows us to calculate the experimental zero-temperature reaction enthalpies, corrected by ZPE: We note that the main uncertainty originates from inconsistent reports of the formation enthalpy of solid Ce 2 O 3 .
Our selected enthalpy corresponds to the careful selection by Konings et al. in Ref. [28]. Taking all their presented formation enthalpies of Ce 2 O 3 into account (ranging from −18.9 to −18.5 eV), results in the following range of zerotemperature reaction enthalpies: ∆ r3 H 0 K = 0.71.1 eV .  probe-charge@1 × 1 × 1 probe-charge@2 × 2 × 2 Alavi-Spencer@4 × 4 × 4 Alavi-Spencer@2 × 2 × 2 Figure 4: Convergence of the Hartree-Fock energy in the limit of an innitely diluted gas (L → ∞). To aid the comparison, the energy in the limit is shifted to 0 for all molecules. Left: validity of the t/extrapolation using the Alavi-Spencer scheme at 2 × 2 × 2 k-points. The dashed lines are tting functions of the form E ∞ + αL −3 + βL −6 . Right: comparison on a broader scale of (i) dierent k-meshes, and (ii) correction schemes to the spurious exchange interaction, for the case of H 2 O. The simplest setting (probe-charge@1 × 1 × 1) exhibits an excessively attractive behavior. It apparently converges to the same limit as the reference setting (Alavi-Spencer@4 × 4 × 4), though much slower, so that a high accuracy cannot be expected from an extrapolation. We employed a plane-wave cuto of ENCUT = 800 eV.   Figure 7: Validity of the t/extrapolation of the correlation energy to an innitely diluted gas. To aid the comparison, the energy in the limit is shifted to 0 for all molecules. Left: MP2 extrapolation using 2 × 2 × 2 k-points and ENCUT = 600. Right: RPA extrapolation using 2 × 2 × 2 k-points and ENCUT = 800. The dashed lines are tting functions according Eq. (13).