Hubbard-corrected oxide formation enthalpies without adjustable parameters

A density functional theory (DFT) approach to computing transition metal oxide heat of formation without adjustable parameters is presented. Different degrees of $d$-electron localization in oxides are treated within the DFT+$U$ approach with site-dependent, first-principles Hubbard $U$-parameters obtained from linear response theory, and delocalized states in the metallic phases are treated without Hubbard corrections. Comparison of relative stabilities of these differently treated phases is enabled by a local $d$-electron density matrix-dependent model, which was found by genetic programming against experimental reference formation enthalpies. This mathematically simple model does not explicitly depend on the Hubbard-corrected ionic species and is shown to reproduce the heats of formation of the Mott insulators Ca$_2$RuO$_4$ and Y$_2$Ru$_2$O$_7$ within $\sim$3% of experimental results, where the experimental training data did not contain Ru oxides. This newly developed method thus absolves from the need for element-specific corrections fitted to experiments in existing Hubbard-corrected approaches to the prediction of reaction energies of transition metal oxides and metals. The absence of fitting parameters opens up here the possibility to predict relative thermodynamic stabilities and reaction energies involving $d$-states of varying degree of localization at transition metal oxide interfaces and defects, where site-dependent $U$-parameters will be particularly important and devising a fitting scheme against experimental data with predictive power would be exceedingly difficult.


Introduction
The Kohn-Sham formulation [1] of density functional theory (DFT) [2] is widely used to model the electronic structures of bulk solids, surfaces, and molecules [3][4][5]. With a reasonable trade-off between accuracy and computational complexity, the possibility to treat the electronic structures of solids as Bloch states of the Kohn-Sham single-particlelike Hamiltonian [6] renders DFT and extensions the de-facto standard for electronic structure simulation in computational materials science [7][8][9][10][11].
Generalized gradient approximations (GGA) and meta-GGAs to electronic exchange and correlation (XC) in Kohn-Sham DFT only depend locally on the charge density, its gradient and Laplacian, and the Kohn-Sham kinetic energy density and yield reasonable cohesive energies with typical errors in comparison to experiment of 0.2 eV per atom for those solids [12,13] where self-interaction errors due to these semilocal approximations are not too large. In materials with strong electron localization, these errors become large, and semi-local DFT thus typically yields significantly larger per atom errors for relative stabilities of e.g. transition metal (TM) oxide phases [14]. The failure of DFT also manifests itself in a qualitatively wrong description of the electronic structure of Mott insulators, where the semi-local, mean-field treatment of XC incorrectly yields a metallic ground state (see Fig. 1).
One approach to reducing self-interaction errors is to admix a fraction of Hartree-Fock exchange energies computed from the Kohn-Sham orbitals to the semi-local XC functional [15]. The computational cost for evaluating the non-local Fock-operator generally increases the cost for this hybrid DFT method by a factor of the system size in comparison to semi-local DFT. In contrast to unscreened Hartree-Fock, which unphysically predicts diverging Fermi velocity and vanishing density of states at the Fermi level for the homogeneous electron gas [16], screened hybrid DFT [17] can treat metallic phases. Magnetic metals are however described worse, e.g., with the widely used HSE06 screened hybrid functional [18] than with semi-local DFT, with overestimated magnetic moments of Fe, Co, and Ni [19,20], incorrectly predicting bcc and hcp ground state structures for Co and Ni [21], respectively, and describing surface reaction energetics worse in particular for magnetic TM surfaces [22]. HSE06 simulations were nevertheless shown to improve on transition metal oxide formation energies over GGA-DFT [23], showing the importance of the reduction of self-interaction errors due to localized electrons in the oxides.
Another, computationally more efficient approach to addressing the self-interaction errors is the GGA+U approach [24][25][26], which does not lead to an overall increase in computational complexity over semi-local DFT. In the GGA+U method, the Kohn-Sham Hamiltonian is supplemented with a Hubbard-like term, which depends on the local projections of the Kohn-Sham single-particle density matrix at ionic sites with strong electron localization. The Hubbard term depends thus explicitly on the Kohn-Sham occupation numbers and can e.g. open up a particle-hole excitation gap in the case of Mott insulators. The strength of the Hubbard term is given by the parameter U , Hubbard-corrected oxide formation enthalpies without adjustable parameters  . Already at half this parameter strength, Fe shows an incorrect suppression of majority spin density of states at the Fermi level E F and tendency towards spin saturation (the magnetic moment is increased by about 30% compared to the U = 0 case). At full U I strength, the magnetic moment increase is 40%. For V 2 O 3 , on the other hand, GGA incorrectly predicts a ferromagnetic metallic phase as the ground state (bandstructure plot folded into the smaller irreducible Brillouin zone of the antiferromagnetic phase), while GGA+U correctly predicts an antiferromagnetic Mott-insulating phase.
which can either be considered as an adjustable, to be fitted parameter or be computed from first principles [27][28][29]. In recent approaches, machine learning techniques were employed to find U -parameters to mimic hybrid DFT bandstructures and bandgaps with GGA+U [30] or to improve agreement with experimental lattice constants, magnetic moments, and bandgaps [31].
A drawback of the GGA+U method is that the ground state energies of Hamiltonians with different Hubbard parameters U at a given ionic site cannot be compared to each other. The compromise of choosing average U -parameters for each transition metal ionic species (irrespective of e.g. differences in oxidation state) leads to an improvement of oxidation energies of transition metal oxides over GGA-DFT [14,32,33]. For molecules, meaningful total energy differences in the case of firstprinciples, reactant-dependent U -values can be computed by integrating the forces acting on the ions along a reaction path [34]. For the prediction of oxidation energies of bulk TM oxides with first-principles U -parameters, the unit cell volume-dependence of the U -parameters can be computed, and integrals over the volume changes required to match the U -values of different oxides can then be used to estimate reaction energies with the oxides at their respective equilibrium volumes [35]. Oxidation energies were found to be further improved using this method rather than using average U -parameters, while oxide formation enthalpies could not be improved [35]. A general problem for either an average U -parameter approach or relying on volume change integrals is that first-principles methods to computing U yield nonzero values also for metallic phases. Delocalized metallic states are however better described without supplemental Hubbard terms (see Fig. 1 for an example of the unphysical effects due to application of a U -term to a metallic bandstructure). Improvement over TM oxide formation energy predictions from screened hybrid DFT was achieved by treating the oxide phases at the GGA+U level with average Uparameters for the TM species and the metallic phases within the GGA [36]. For each Hubbard-corrected TM species, an energy offset was fitted to experimental formation enthalpies, so that GGA+U total energies plus offset can be compared to the GGA total energies for the metallic phases. These fitted offsets allow the resulting corrected total energy differences between GGA+U and GGA calculations to be interpreted as reaction energies. In a similar approach, elemental reference energies were fitted directly in Ref. [37] against experimental formation enthalpies with TM oxides treated at the GGA+U level with average U -parameters, and formation energies predicted with this approach were used to fit a multivariate linear regression model of formation energies exclusively against atomic properties [38]. With the availability of sufficient experimental reference data, the approach of combining GGA and GGA+U total energies can be further refined by fitting offset energies not only for different TM species, but also for different oxidation states and different ligand coordination [39]. While experimental calorimetry reference data for thermodynamics of stoichiometric bulk TM oxides can be used here to parametrize these element-or even system-specific corrections, determining the required offset energies for Hubbard-corrected TM sites near e.g. point defects or interfaces will be difficult.
Therefore, a newly developed method free of adjustable parameters is presented here that combines the treatment of d-electrons of varying degrees of localization within the GGA+U approach with site-specific, first-principles U -parameters from linear response for TM oxides with GGA calculations optimal for metallic phases. To correct for the offset between GGA+U total energies with different U -parameters and GGA (i.e. U =0) total energies, genetic programming with experimental TM oxide formation energies as references is applied to find a mathematical model that does not explicitly depend on the Hubbard-corrected TM species. The resulting, relatively simple mathematical model does thus not have adjustable, system-specific parameters. Examples of the applicability of the model for TM oxides with Hubbard-corrected TM species not contained in the training data are shown.

Method
The method consists of computing DFT and DFT+U total energies for metal and oxide phases, respectively, as described below and using a genetic programming model for finding an energy offset for the Hubbard-corrected oxide total energies, such that physically meaningful energy differences between metals and oxides can be obtained (Fig. 2).  Figure 2. Schematic of the presented method to computing TM oxide heat of formation without adjustable parameters. TM oxide phase total energies are computed using GGA+U with sitedependent, first-principles U -parameters, and metallic phases and the O 2 gas phase reference are treated without Hubbard corrections. The O 2 total energy reference is corrected adding a constant offset of ∼0.84 eV as described in the main text. Based on energy offsets from a genetic programming model, the total energies from GGA+U with site-dependent U -parameters are adjusted, so that comparison to the GGA total energies of the metallic phases becomes possible and reaction energies involving both oxides and metals can be computed.

DFT calculations
GGA and GGA+U simulations are performed with the VASP code [40] using plane-wave basis sets with a cut-off energy of 600 eV and representing ionic cores by projectoraugmented wave (PAW) [41] frozen-core pseudopotentials [42]. Brillouin zones are sampled with Γ-point centered, regular grids with a spacing of at most 0.02Å −1 , and Gaussian smearing with a width of 0.05 eV is used to determine Kohn-Sham occupation numbers. The widely used PBE functional [43] is employed for the GGA to electronic exchange and correlation. Forces acting on ions and stress on the unit cells are simultaneously relaxed, with residual forces of less than 0.02 eV/Å per atom and residual pressures of less than 0.02 GPa, respectively. The total energy of the O 2 molecule is described poorly by GGA approaches, which is typical for multiply bonded molecules [27,44]. Empirical approaches fit the total energy of O 2 to experimental reaction energies, such as oxidation energies of solids [14] or the formation of H 2 O [45]. Here, instead of empirically fitting to the experimental formation energy of H 2 O, atomization energies as computed with the first-principles Weizmann-4 computational thermochemistry protocol (W4) [46] are referenced. Using the difference of the atomization energies of O 2 and H 2 O the total energy of O 2 used in this study is computed as (3) and such that Tr ρ σ = m=− ρ σ mm is the trace over the orbital degrees of freedom of the projected density matrices (2 +1 = 5 in the case of d-orbitals). σ denotes the electron spin degree of freedom (spins assumed to be collinear in this study), and N I = σ N σ I = σ Tr ρ σ I is the total projected occupation at site I.Ū I andJ I are the spherically averaged screened on-site Coulomb repulsion and exchange interactions at site I, respectively. Since Eq. 6 depends on the difference between these interactions at each site I, it is convenient to introduce an effective Hubbard parameter U I with Hubbard-corrected oxide formation enthalpies without adjustable parameters 7 The Kohn-Sham system is exposed to the following additional potential due to the Hubbard correction at site I: The Hubbard corrections (6) effectively constitute an energy penalty for non-integer eigenvalues of ρ I,σ . Considering Eq. 8 in the eigenbasis of ρ I,σ , the potential is positive for projected occupation eigenvalues less than 1/2, favoring a lowering of the eigenvalue towards zero. For eigenvalues larger than 1/2, the potential is negative thus favoring increasing the eigenvalue towards one. Accordingly, the Hubbard corrections with their explicit occupation number dependence can open a Mott-Hubbard gap in the Kohn-Sham spectrum of an otherwise metallic system in a single-particle bandstructure picture.

Hubbard-U from first principles
Rather than treating the effective Hubbard parameters U I as empirical fitting parameters, the linear response approach to determining U I by Cococcioni and de Gironcoli [28] was implemented here on top of the existing DFT+U code in VASP. DFT total energies spuriously show a non-linear dependence on fractional electron numbers (fractional particle numbers correspond to quantum mechanical ensemble averages, and the total energies should show linear behavior between pairs of adjacent integer particle numbers). The U I can be chosen to compensate spurious curvature with respect to changes in the Hubbard site occupations N I . Instead of constraining Hubbard site occupations to determine these curvatures (which is difficult with a plane-wave basis set), the linear response approach adds perturbations to the potentials (8):Ṽ I,σ mm = V I,σ mm + α I δ m,m . The functional (3) is Legendre transformed such that the energy is minimized with respect to the potential perturbations α I instead of the Hubbard-site occupation numbers, and the Hubbard parameters U I are then determined as [28]: χ IJ = ∂N I /∂α J is a density response function of the total occupations N I of each Hubbard site I to perturbations from the α J for any of the Hubbard sites J and is computed here with finite differences in α J , letting the electronic structure selfconsistently adjust to the perturbed Hubbard potentials. Not letting the Kohn-Sham potentials self-consistently adjust to the perturbed Hubbard potentials, but only computing the Kohn-Sham bands and thus the Kohn-Sham density matrix at the fixed perturbed potentials leads to the perturbed Hubbard site occupations N 0 I . The finite difference χ 0 IJ = ∆N 0 I /∆α J is used to compute the corresponding bare response function. While the updated Kohn-Sham potentials can self-consistently screen the effect of the perturbed Hubbard potentials, the bare response accounts for single-particle bandstructure-only rehybridization between the Hubbard sites and their surroundings. The resulting curvature in total energy due to these latter rehybridizations is not a measure of a self-interaction error to be corrected for with a Hubbard-term, and Eq. 9 thus removes this curvature from the computation of the U I .

Ab initio U -parameters in practice
In the present implementation, the unperturbed electronic structure is computed first, then the perturbation due to α J is added, and the Kohn-Sham bands are updated yielding the bare response χ 0 IJ for all sites I, and finally this same simulation is continued allowing the Kohn-Sham potentials to adjust yielding the response χ IJ for all sites I. The Spglib code [48] is used to identify symmetrically equivalent Hubbard sites and the corresponding symmetry operations, so that simulations only need to be performed for perturbations at nonequivalent sites and the full χ and χ 0 response matrices can be obtained using symmetry operations.
To decouple the effect of the perturbation from its periodic images, the TM oxide unit cells are repeated here into supercells separating the periodic site images by at least ∼8Å. Coccocioni and de Gironcoli [28] find that extending the Hubbard sites by a fictitious background site for charge neutrality of the perturbations improves convergence with respect to supercell size. The response function matrices χ and χ 0 are thus extended by one row and one column, and the additional matrix elements are determined by the charge neutrality requirement of the sums over each row and each column of the matrices yielding zero. With this construction, the extended response matrices now have a zero singular value and cannot be inverted. This zero singular value is due to the fact that adding a constant offset to all Hubbard site potentials and the background does not lead to any total energy curvature. Here, the Moore-Penrose pseudo inverse is used instead of the matrix inverses in Eq. 9, such that this zero singular value leads to a zero singular value for the pseudo inverse, correctly not affecting the computed U I .
After the U I have been computed in supercells, GGA+U optimizations of the primitive cells with respect to atomic forces and unit cell stress are performed. Based on the optimized structures new supercells are created to re-compute the U I and then to re-optimize the primitive cells again. Performing this relaxation and U I determination cycle twice led to mean average TM oxide lattice constant changes < 10 −2Å , and thus no further structural-dependence self-consistency cycles of the U I have been performed.
Generally, the U I are computed here from GGA+U calculations with all Hubbard parameters set to zero, yielding the spurious energy curvatures to correct for. Some of the considered TM oxides are incorrectly predicted as metallic at the GGA level. For these oxides, GGA+U simulations are performed at a range of Hubbard parameters U input large enough to open a bandgap, and, following the idea of Ref. [49], the computed U I values are linearly extrapolated to U input = 0 (see Table S1 for a list of all computed U I values, indicating which TM oxide cases are based on this extrapolation technique).
The values of the Hubbard-parameters U I determined via Eq. 9 have a dependence on the choice of the projectors of the Kohn-Sham density matrix on the Hubbard sites. The VASP code uses the projections of the Kohn-Sham pseudo wavefunctions onto the PAW projector functions, which are non-zero only inside the augmentation spheres around the ions [50,51]. The DFT+U implementation in the Quantum Espresso planewave pseudopotential DFT code, e.g., offers further choices for the definition of the projection onto the Hubbard sites, such as the d-orbitals of free atoms with or without orthogonalization with respect to the atomic orbital overlap between Hubbard sites [52]. Using the density functional perturbation implementation [53] in the Quantum Espresso suite to computing χ and χ 0 with ultrasoft pseudopotentials [54] from the GBRV dataset [55], the computed U for bulk Fe is about 7.0 eV and 5.5 eV with and without projector orthogonalization, respectively. Similarly, U -parameters obtained from constraining local occupations were found to depend strongly on TM ion muffin-tin sphere radii [56]. These dependencies show the importance of determining the parameters U I with the same projectors as are being used in the DFT+U simulation.
Similarly large U -parameters for bulk Fe based on constrained DFT are reported in Ref. [57], and also the VASP PAW implementation (without inter-site projector overlap due to non-overlapping PAW augmentation spheres) finds a large value of about 7.2 eV. Even only applying half this Hubbard-parameter strength in a GGA+U simulation of bulk Fe leads to unphysical suppression of majority spin density of states around the Fermi level (see Fig. 1) and the magnetization is spuriously driven towards saturation. As the corresponding d-band shifts would deteriorate the description of, e.g., transition metal surface chemistry [22], this example highlights that GGA+U simulations of metallic systems with first-principles U -parameters should be avoided because of potentially unphysical predictions.
A promising direction for (at least partial) compensation of unphysical effects on metallic electronic structures in DFT+U is the DFT+U +V approach [58], where also inter-site Hubbard corrections, determined from the off-diagonal elements in (9), are applied. The DFT+U +V method is implemented in Quantum Espresso, and for the above orthogonalized projector case with an on-site Hubbard parameter of 7.0 eV, the V -parameter between nearest neighboring sites is only ∼0.3 eV, barely affecting the electronic structure qualitatively and both GGA+U and GGA+U +V calculations yielding a spuriously saturated, integer magnetization of ∼3 Bohr magnetons per Fe atom. Generally, the DFT+U +V method can be extended to such off-diagonal Hubbard terms between d-and p-projections, also on the same site [58], which could improve suppression of unphysical effects of Hubbard-U on metallic systems, but here only onsite Hubbard-corrections based on d-projections are considered and implemented. All results presented in the following treat metallic phases without Hubbard corrections.

TM oxide and metallic phases
To find a model enabling the energetic comparison of products and reactants with first-principles U -parameters rather than average U -values, experimental formation enthalpies of transition metal oxides requiring Hubbard corrections are used as training data with corrections as described below. The formation enthalpies of 66 such TM oxides containing the TM ions V, Cr, Mn, Fe, Co, Ni, and Mo are taken from Ref. [59]. A Debye model fitting technique introduced and benchmarked in Ref. [60] is used to determine T = 0 K enthalpies from the experimental values at 298 K. Neglecting thermal expansion of the solid phases, the heat capacity at constant pressure and the entropy are fitted numerically to match the T = 298 K entropies and heat capacities provided in Ref. [59] (c p (298 K) and S(298 K) for ZnFe 2 O 4 taken from Ref. [61] due to unphysical c p (298 K) in [59]), both for the TM oxides and the metallic reference phases, considering the Debye temperature T D and the pre-factor y as fitting parameters.
Here, the fitted T D is furthermore employed to provide an experimental estimate of the phonon zero-point energy E ZP per formula unit using that the Debye model density of states is quadratic in the phonon frequency ω below the Debye cut-off frequency ω D = k B T D / : N A is the number of atoms per formula unit, the reduced Planck constant, and k B Boltzmann's constant. The respective enthalpy changes 298 K 0 c p (T )dT and the relatively small E ZP are subtracted from the formation enthalpies for each of the solid product and reactant phases (for the metallic reactant phases considered here, the enthalpy of formation at 298 K is zero by definition).

O 2 reference
From the zero enthalpy reference of O 2 at standard conditions the integral 298 K 90 K c p (T )dT is subtracted (heat capacity data valid in this temperature range taken from Ref. [62]). For the remaining integral from 0 K to 90 K to be subtracted, a constant c p = 5R/2 (R: molar gas constant) is assumed. Furthermore, the zero-point energy of O 2 as computed with the W4 quantum chemistry protocol is subtracted [46].
Assuming some degree of cancellation of pressure times volume (pV ) terms between solid reactants and products [60], pV terms are neglected for the solid phases. The formation energy training targets are thus given as the experimental enthalpy differences between the TM oxide and metallic phases and O 2 with heat capacity integrals and zeropoint energies subtracted for all phases, as described above.

Ru oxide benchmarks
As test cases outside the training and validation data, two correlated Ru oxides with available experimental thermodynamic data are considered: Ca 2 RuO 4 and Y 2 Ru 2 O 7 . The experimental enthalpy of Ca 2 RuO 4 at 1140 K relative to 2 CaO and RuO 2 is taken from Ref. [63]. The formation enthalpies of CaO and RuO 2 and their entropies and heat capacity data for Debye model fits and c p (T ) integration are taken from Ref. [59], so that O 2 and solid Ca and Ru can be used as reference phases (with zero-point energies from Debye model fits to data in [59] subtracted from the Ca and Ru references and the same adjustments for O 2 as described above). The Debye temperature of Ca 2 RuO 4 is determined from the curvature in low-temperature c p (T )/T measurements in Ref. [64] to be T D ≈418 K. The entropy of Ca 2 RuO 4 is computed as S(298 K) = The experimental enthalpy of formation of Y 2 Ru 2 O 7 at 298 K is taken from Ref. [66]. Debye model fits for the corrections of Y and Ru are based on data in [59], and also O 2 is corrected analogously to the other oxide heats of formation. A low-temperature heat capacity fit yielding a Debye temperature of 449 K for Y 2 Ru 2 O 7 is taken from Ref. [67]. Low and moderate-temperature heat capacity data from Refs. [68] and [66], respectively, is used in the integral (11) to determine the entropy at 298 K of ∼216 JK −1 mol −1 and fit the parameter y of the Debye model for T = 0 extrapolation and subtraction of E ZP . The resulting benchmark Ru oxide formation energies and the 3d-TM and Mo oxide formation energies used for the model search are listed in Table S2 in the SI.

Model search
The ansatz for oxide formation energies in the model search used here is based on an energy offset E off (U I , ρ I ) for every correlated site I: where the sum over ν accounts for the N metals metallic references calculated at the GGA level of theory and the factor β ν for the stoichiometric occurrence of these metals in the oxide. γ accounts for the stoichiometry of O in the oxide, and E(O 2 ) is the corrected value from Eq. 2. The offset at site I is chosen to be of the form f GP is a function of the sum of the eigenvalues of the projected density matrix for site I and of the sum of the squares of these eigenvalues, ensuring rotational invariance of the model. Sums of higher eigenvalue powers as additional arguments for f GP did not lead to improved models, and models depending on higher eigenvalue powers hence are not presented in the following.
As f GP has no explicit site-dependence, site indices I are dropped in the following from the density matrices and Hubbard parameters for simplicity, noting that evaluating a model always corresponds to a sum over all correlated sites I in the oxide. Starting with the initial guess corresponding to E off (U, ρ) equaling the density matrix-dependent part E U (U, ρ) (Eq. 6) of the DFT+U functional (3), genetic programming [69,70] is used to evolve generations of programs {f z GP } with the gplearn code [71]. As a measure of fitness for ranking the models f z GP , the Akaike information criterion (AIC) [72] is implemented: where L is the maximum of the log-likelihood function for the considered model and Φ is the number of free parameters. Under the assumption of n normally and independently distributed observations a i , the maximum log-likelihood of a corresponding Gaussian regression model can be written as [73]: b i are the model predictions. Choosing a different energy unit than eV for normalization in Eq. 18 merely leads to a constant offset of the AIC, not affecting the relative ranking of the models. With decreasingσ 2 , i.e., with increasing goodness of fit, the term −2L normal in the AIC (16) decreases. The term 2Φ, on the other hand, is an overfit penalty for the number of parameters of the model. The selection strategy is thus to prefer models with small AIC. The observations are experimental oxide formation energies and the computed values from Eq. 13 the model predictions, both normalized for each oxide as described in the following. The more atoms there are per formula unit, the larger typically the error in the predicted formation energy of the compound. Cohesive energies of solids that are well described within the GGA, e.g., have typical errors of 0.1 eV-0.2 eV per atom [12,13]. The considered oxides have a relatively wide range of ratios κ of number of atoms N A per formula unit to Hubbard-corrected sites N H per formula unit ranging from 2 to 8, and oxides with larger κ typically will have a larger GGA contribution to the error in the predicted formation energy per formula unit. If a i and b i were normalized by dividing by N H or by number of formula units, residual errors for oxides with larger ratio κ or larger formula units would be weighted more strongly in Eq. 18, respectively. The Table 1. Complexity weights w j for variables, numerical coefficients, and mathematical operations in the programs. Divisions have higher weight than additions and multiplications, to favor mathematical simplicity of the programs. The variables are var 1 = σ Trρ σ and var 2 = σ Tr(ρ σ ) 2 , both with a weight of one. No weight is assigned for taking the square root of these variables. Other powers of the variables can only be expressed as products (or divisions) with corresponding penalty weights for multiplications and variable occurrences.
variable integer non-integer coeff. add, subtract multiply divide f z GP could potentially be trained to compensate for relatively large GGA contributions to the errors, which is undesired. If a i and b i instead included a factor 1/N A , the E off (U, ρ) would have a small relative weight for oxides with large ratio κ. As a compromise, the experimental and computed oxide formation energies per formula unit are divided by √ N A , ensuring that all E off (U, ρ) have appreciable, balanced relative weight in the residuals a i −b i . With these normalized residuals a i −b i , the correspondingly normalized mean-absolute error (MAE) is defined: which will be presented in addition to the model AIC in Sec. 3. The complexity term 2Φ in Eq. 16 is generalized to not only account for and penalize the number of numerical parameters, but also mathematical operations: w j is a weight for a mathematical operation, variable, or numerical coefficient (see Tab. 1) and ϕ j the corresponding occurrence in the considered model (or program) f z GP . Using the resulting AIC as a fitness function evaluated on a random selection for each of the programs of half of the 66 TM oxide reference data, k-tournament selection [74] is employed as implemented in the gplearn code [71]. In this procedure, k programs out of the current generation (here k=25 out of a total population size of 100 per generation) are randomly selected, from which the program with lowest AIC is selected for crossover and mutation operations (operations named in analogy to biological genetic evolution). With 70% probability, the winner from a second k-tournament is used for crossover of terms between the two programs, or with 25% probability, random mutations are performed as described in Ref. [71], and the resulting new program becomes part of the next generation. With 5% probability, the tournament winner is promoted to the next generation unchanged. The initial generation consists of 30 copies of program (15) and 70 randomly generated f z GP , and the programs are evolved for 10 5 generations.
Beginning from the last generation, 10 5 further generations are evolved, where now all programs in each generation are optimized numerically (using the Nelder-Mead simplex algorithm [75,76] as implemented in scipy [77]) and simplified symbolically (using the sympy code [78]) in addition to the genetic operations. Fractions in programs are normalized to equivalent fractions so that the coefficient of the lowest power in the variables in the denominator is one. To constrain the search space for the final 10 5 generations, a numerically large penalty is added to the AIC for models that do not vanish for integer density matrix eigenvalues, i.e., where E U vanishes, the DFT and DFT+U energy functionals coincide, and thus no correction term for total energy comparison is required. A modified version of the gplearn code with unique weights for mathematical operations, coefficients, and variables and integrated symbolic and numerical optimization can be found on-line [79].

Results and discussion
From the constrained genetic programming search for models vanishing for integer density matrix eigenvalues, the programs with lowest AIC at a given complexity from the 10 5 generations are shown in Fig. 3. All models found satisfying this constraint contain the factor σ Trρ σ − Tr(ρ σ ) 2 , which vanishes for integer eigenvalues of ρ σ .
The simplest program is to subtract E U from the GGA+U total energy, yielding a relatively high ∆Ẽ MAE of about 0.5 eV. Multiplying this model with a pre-factor increases the complexity Φ to 6. The resulting model 1.17 · ( σ Trρ σ − Tr(ρ σ ) 2 ) lowers ∆Ẽ MAE by only 6 meV, while the AIC increases by almost 6, thus identifying this model as an overfit. Significant improvements of the AIC and ∆Ẽ MAE are achieved by programs further multiplying E off by σ Tr(ρ σ ) 2 or σ Tr(ρ σ ) 2 / σ Trρ σ , respectively. The program with lowest AIC and hence optimal goodness of fit vs. complexity trade-off yields the correction term where the site index I is not suppressed to emphasize that this model with site-specific U I Hubbard parameter strength is to be subtracted from GGA+U total energies for every Hubbard-corrected site I in the system. The next more complex model in Fig. 3 has a non-integer and numerically optimized coefficient in the denominator (leading to a complexity increase of 2), but the optimized coefficient is almost an integer, and the model is approximately equal to (21), and the AIC increases. The most complex model in Fig. 3 further has slightly different numerical coefficients for σ Trρ σ and ( σ Trρ σ ) 2 in the denominator, which increases the AIC, identifying a tendency towards overfitting. Eq. 21 is thus the model identified as optimal for computing offsets to GGA+U total energies for comparison at different U -values, including (U = 0) GGA total energies. E opt off only depends on the density matrices through the difference ∆ = σ Trρ σ − Tr(ρ σ ) 2 , and this dependence is depicted in Fig. 4a. The model can be considered a rescaling of E U . For differences 0 < ∆ < 1.36, E opt off is increased over E U , while for larger differences, it is decreased relative to E U .
In Fig. 4b, the performance of the model is compared to GGA+U calculations with average, fitted U -values. The average U -values are taken from the Materials Project [9] and are listed in Tab. S3. Constant energy offsets for each TM species are fitted to the binary 3d and Mo oxide heats of formation, following the procedure in Ref. [36]. With the seven different TM ions considered and with the empirical, fitted U -values and energy offsets, there are thus effectively 14 fitting parameters. E opt off with the first-principles, site-dependent U -values and thus significantly fewer degrees of freedom (through its functional form and numerical coefficient) leads to similar residual errors in comparison to the experimental benchmarks, showing generally good agreement (experimental and computed TM oxide heats of formation listed in Tab. S2). With the similar quality of prediction for these bulk oxide formation benchmark cases, the method presented here holds the promise of enabling predictions of TM oxide formation energies and other reactions energies between transition metals and their oxides where experimental benchmark data for fitting is lacking.
The heats of formation of the ruthenates Ca 2 RuO 4 and Y 2 Ru 2 O 7 serve here as well-defined benchmark cases; other Ru oxides, such as RuO 2 , CaRuO 3 , SrRuO 3 , and Sr 2 RuO 4 are not considered here, as these oxides are metallic [80] (Sr 2 RuO 4 is metallic at temperatures below ∼130 K and superconducting below ∼1 K [81]) and should thus not be treated at the GGA+U level of theory. As there are only two Hubbard-corrected Ru oxides considered here, no energy offset nor average U -value is fitted for Ru, but only   the first-principles U -value approach based on E opt off is applied to compute the Ru oxide heats of formation. Calculation of the experimental benchmark values is explained in Sec. 2.3.3.
Ca 2 RuO 4 is an insulator at ambient conditions [82]. With powder diffractionbased structural refinements of the P bca-phase at 295 K from Ref. [83] as a starting guess, the structure of Ca 2 RuO 4 is optimized iteratively at the GGA+U level with self-consistent U -parameters as described in Sec. 2.2.3. Antiferromagnetic order with a propagation vector of (1, 0, 0) (corresponding to the A-centered mode [83]) is considered (the energy of the B-centered mode with propagation vector (0, 1, 0) is found to be within ∼10 −4 eV/atom and is thus energetically almost indistinguishable from the A-centered phase within the accuracy of the simulations, allowing the lattice to relax for both the A-and B-centered modes). A Hubbard-parameter strength U Ru I ≈4.5 eV is calculated with Eq. 9. With the corresponding GGA+U total energy of Ca 2 RuO 4 and GGA total energies for the reactants and applying corrections (2) and (21), a heat of formation of −2.29 eV/atom is predicted, which is in good agreement with the experimental result of −2.35 eV/atom. Y 2 Ru 2 O 7 has pyrochlore structure (spacegroup F d3m [84]) and is an insulator with  Fig. 5 were compared energetically at the GGA+U level. The low-energy antiferromagnetic phase used in the following is depicted in Fig. 5 (dashed lines indicate the primitive magnetic cell containing four formula units). After iterative structural optimization and computation of the first-principles Hubbard-parameter strength, a parameter U Ru I ≈4.3 eV is found. The computed heat of formation is −2.36 eV/atom, comparing well to the experimental estimate of −2.44 eV/atom. This 0.1 eV/atom error is within the range of errors found in Ref. [14] for the heat of formation of oxides not containing transition metal ions (and thus not requiring Hubbard corrections) with a fitted O 2 gas reference, and within the range of GGA reaction energy errors involving metallic TM phases in general [12,13,22]. While oxidation state-specific TM references fitted to experiments can reduce per atom oxide reaction energy errors significantly [39], having established here an approach to correcting mixed GGA and GGA+U total energy differences with sitespecific, first-principles Hubbard parameters, which leads to errors typical for GGA-only estimates for reaction energies without localization errors, will enable oxide reaction energy predictions in the absence of experimental reference data.

Conclusion
The presented approach enables the prediction of transition metal oxide formation energies and other reaction energies involving oxide and metallic phases with Hubbard corrections computed from first principles.
Comparison of total energies from Hamiltonians with different Hubbard parameter strengths including U = 0 is facilitated with a site-projected density matrix-dependent model, that was found through genetic programming against experimental transition metal oxide formation energies. This mathematically simple model does not involve fitting parameters to be adjusted to transition metal ion types or oxidation states. Unlike previous Hubbard-corrected approaches to TM oxide heat of formation prediction requiring experimental benchmark data to fit energetic corrections and average Hubbard strengths, the method presented here can thus also be used to compute transition metal oxide reaction energies when such experimental fitting targets are lacking, which will be of particular importance for oxide interfaces and defect chemistry.  Table S2: Transition metal oxide formation energies per formula unit computed using the presented genetic programming (GP) model with site-dependent U -parameters, using fitted metallic references (metref) and site-independent U -parameters, and as measured experimentally (exp). Ca 2 RuO 4 experimental formation enthalpy from [3]; low-temperature and moderate temperature thermodynamic data from [4] and [5], respectively, used for extrapolation to 0 K. Y 2 Ru 2 O 7 formation enthalpy from [6] and low-temperature and moderate temperature thermodynamic data from [6][7][8]. Thermodynamic data for O 2 from [9]. Quantum chemistry corrections and zero-point energy of O 2 from [10]. Thermodynamic corrections for ZnFe 2 O 4 from [11]. All other experimental thermodynamic data from [12]. Zero-point energies E ZP subtracted from experimental energies. E ZP for solids from Debye model fits. See main article for details.

Acknowledgment
3. Average, site-independent U-parameters from the Materials Project  Table S3: Average fitted, site-independent (and thus also oxidation state-independent) U -parameters used in the Materials Project [13].