Determination of Formation Energies and Phase Diagrams of Transition Metal Oxides with DFT+U

Knowledge about the formation energies of compounds is essential to derive phase diagrams of multicomponent phases with respect to elemental reservoirs. The determination of formation energies using common (semi-)local exchange-correlation approximations of the density functional theory (DFT) exhibits well-known systematic errors if applied to oxide compounds containing transition metal elements. In this work, we generalize, reevaluate, and discuss a set of approaches proposed and widely applied in the literature to correct for errors arising from the over-binding of the O2 molecule and from correlation effects of electrons in localized transition-metal orbitals. The DFT+U method is exemplarily applied to iron oxide compounds, and a procedure is presented to obtain the U values, which lead to formation energies and electronic band gaps comparable to the experimental values. Using such corrected formation energies, we derive the phase diagrams for LaFeO3, Li5FeO4, and NaFeO2, which are promising materials for energy conversion and storage devices. A scheme is presented to transform the variables of the phase diagrams from the chemical potentials of elemental phases to those of precursor compounds of a solid-state reaction, which represents the experimental synthesis process more appropriately. The discussed workflow of the methods can directly be applied to other transition metal oxides.


Introduction
Due to their exceptional electronic, magnetic, or optical properties [1], transition metal oxide (TMO) compounds are key components of many modern technologies. Next to many other applications, TMOs are utilized as active anode and cathode materials in Li-and Na-ion batteries [2][3][4][5] and as electrodes in solid-oxide fuel and electrolyzer cells [6][7][8] (SOFC and SOEC), where they enable the catalytic reactions with oxygen. The decisive functional parameters of a TMO, such as the catalytic activity and the electronic conductivity in the case of a solid-oxide cell electrode, can be tuned by varying the stoichiometry, e.g., by external doping or the incorporation of intrinsic lattice defects [9][10][11]. A phase diagram with respect to elemental reservoirs determines the ranges of the experimental synthesis conditions within which the desired compound can be stabilized and within which its composition can be varied without unwanted competing phases being formed.
In order to derive such a phase diagram theoretically, formation energies need to be known of all compound phases, which consist of a subset of the involved elements. Highly accurate formation energies can, in principle, be calculated using the methods of density functional theory (DFT) [12]. However, the determination of formation energies of TMOs using common

•
The molecular state is the natural reference phase for oxygen. The well-known over-binding of the O 2 molecule in LDA and GGA [17] introduces an error in the formation energies of oxides. A correction scheme was proposed by Wang et al. [13], which builds on a comparison between the theoretical and experimental values of formation energies for a series of non-TM oxide compounds.
In the Materials Project (MP) database [18], for example, this scheme is applied and denoted by the term anion corrections.

•
The uncompensated electronic self-interaction imposed by approximate exchange-correlation functionals immanent in DFT methods, especially in the case of TMOs, where strongly correlated TM-d-electrons form the valence band, leads to incorrect total energies and underestimated band gaps. There are different approaches to cope with this inaccuracy, such as the use of hybrid functionals [19], a self-interaction correction (SIC) [20], or a Hubbard-U correction that acts on the d-electrons of the TM as an effective potential (DFT+U) [21].

•
In case that the DFT+U method is used to obtain a corrected total energy of a TMO, the formation energy contains an error if the total energies of the elemental reference phases were calculated by LDA or GGA, as it is generally done for the elemental phases and for compounds not containing the TM elements. The error can be systematically corrected by applying a method worked out by Jain et al. [14], where the total energies from the DFT+U calculations are shifted by a constant amount per TM atom. The approach is used in the MP database [18], denoted by the term advanced corrections, to obtain formation energies of TM-containing compounds.
This paper aims at providing a complete picture of the derivation of formation energies and phase diagrams of TMOs based on DFT+U and the proposed correction schemes by evaluating, discussing, and generalizing the approaches proposed in the literature. After a description of the employed methods in Section 2, we first revisit the correction of the oxygen reference energy (Section 3.1). The approach described in the literature [13] is extended by taking into account a larger set of non-TM oxides from the alkali, alkaline earth, boron, and carbon groups of the periodic table of elements in their ground state structures. The uncertainty of the resulting correction value is quantified and taken into account in the following calculations. Explicitly for Fe-containing oxides, we review in Section 3.2 the commonly used scheme [14] to correct for the error arising when comparing total energies from the DFT+U and DFT calculations. The underlying assumption, that the deviation between calculated and experimentally determined formation energies linearly approaches zero with the decreasing TM content, is checked by applying a general linear function and including ternary compounds with lower TM contents in the procedure. We performed the analysis for different U values acting on the 3d orbitals at Fe sites, which directly led us to a way of finding an optimal U value by minimizing the total mean squared deviation between the experimental and calculated formation energies for the considered compounds (Section 3.3). The widely used U Fe value of approximately 4 eV [13,[22][23][24][25][26] is confirmed, and it provides reasonable band gap predictions for the considered Fe oxides (Section 3.4).
Using the corrected formation energies, we determined the phase diagrams for LaFeO 3 , NaFeO 2 , and Li 5 FeO 4 , which are technologically promising as materials for electrodes in SOFC/SOEC devices [10,11], as positive electrodes in Na-ion batteries [27,28], and as multi-redox active cathodes in Li-ion batteries [29], respectively. In the cited literature, these compounds are described to be synthesized by solid-state reaction routes. This is a well-established technique to process ternary and higher-component oxide compounds, where oxidic precursor phases are mechanically mixed and exposed to an oxygen atmosphere of variable temperature and pressure [29]. In order to reflect this process in phase diagrams, we present these in Section 3.5 not in the conventional way with respect to elemental reservoir energies of the nongaseous components but with respect to the chemical potentials of oxidic solid precursor compounds and molecular oxygen gas. The appendix provides a table listing experimental reference values and calculated formation energies (Appendix A), a heuristic explanation of the dependence of oxidation energies on U (Appendix B), a formal comparison of the method for determining an optimal U described in Section 3.3 and the approach using oxidation energies described by Wang et al. [13] (Appendix C), and a detailed description of the transformation of variables in the phase diagrams (Appendix D).

DFT(+U) Calculations of Total Energies
Total energies of elemental and compound phases were calculated using the method of density functional theory (DFT). The plane-wave-based DFT code VASP [30] was used for this purpose, with strict convergence criteria that ensure accurate results. A cutoff energy of 600 eV was set for the plane-wave basis functions describing the valence electrons. The partial occupancies were set according to the linear tetrahedron method with Blöchl corrections [31], except for the O 2 molecule, for which Gaussian smearing was applied with a width of 0.05 eV. The interaction with the core electrons was modeled by projector-augmented-waves (PAW) pseudopotentials [32]. The generalized gradient approximation (GGA) of Perdew et al. [33] was chosen for the exchange-correlation (xc) functional. Electronic self-consistency loops were stopped when the energy difference between two steps was less than 10 −5 eV, and the structures were relaxed until the minimal force component acting on an atom was below 10 −4 eV/Å. Brillouin zone integrations were performed on grids with about 40 3 /V k-points, with V denoting the initial supercell volume in Å 3 . Cell volumes were optimized by total energy minimization, as implemented in VASP. All calculations were performed for cells with periodic boundary conditions. For the calculations of compounds containing Fe, the Hubbard-U correction of Dudarev et al. [34] was applied to the Fe 3d orbitals in order to account for the artificial self-interaction and, concomitantly, the too-weak localization of the strongly correlated 3d electrons. Spin-polarization was taken into account, and for the following phases, the initial magnetic moments were set up according to the known ground-state spin configurations: Fe (ferromagnetic (FM)), FeO (antiferromagnetic (AFM)), Fe 2 O 3 (AFM), Fe 3 O 4 (ferrimagnetic), and LaFeO 3 (AFM). All other considered iron compounds were set up in FM spin configurations.

Correction of the O 2 Over-Binding
Equation (1) defines the formation energy of one formula unit of a binary oxide compound A l O n with respect to the elemental phases of a nonoxygen component A and oxygen: E total denotes the total energy per formula unit of the compound. The total energies per atom of the elements in their ground states are expressed by the chemical potentials µ (0) . Equation (1) can straightforwardly be adapted to systems containing multiple nonoxygen components, e.g., A l A l . The natural reference state of oxygen is the O 2 molecule, with the energy If the energy of an isolated O atom is set to zero, µ (0) (O 2 ) is the binding energy of the molecule. Its absolute value is known to be overestimated in DFT using LDA or GGA, corresponding to a too-strong binding [17]. Wang

Energy Correction Method for the Combined DFT and DFT+U Results
In addition to the O 2 over-binding error, a systematic inaccuracy is introduced if one calculates the formation energies of oxide compounds containing a transition metal element M, e.g., binary systems M m O n or general compounds A l M m O n , by comparing the total energies of the oxide compounds calculated by DFT+U to the total energies of the elemental phases, which are generally determined by DFT. The application of U on the 3d orbitals of the transition-metal atoms M in the compound shifts the reference energy. Hence, the results of calculations involving M obtained by the two different methods cannot directly be compared. As described by Jain et al. [14], the value of this error can again be systematically estimated by comparing calculated to experimental formation energy values for a set of compounds {A l M m O n }. The concept behind this approach is that the formation energy differences per atom, between the combined DFT and DFT+U results and the experimental values are expected to depend linearly on the fraction x M := m/N fu of the transition metal M in the formula unit, which contains a total number of N fu = l + m + n atoms. For the case of iron, i.e., M = Fe, the binary oxide compounds FeO, Fe 2 O 3 , and Fe 3 O 4 were taken into account in the original work [14], and their energies were calculated using a given, independently determined U Fe value of 4 eV [13]. The data points ∆e were fitted by a line through the origin, ∆e fit,0 (x Fe ) = m 0 · x Fe , justified by the argument that, without any transition metal atoms in the compound, and considering an O 2 over-binding correction as described in Section 2.2, there should be zero deviation between the calculated and experimental results. This opens up the possibility to correct the DFT/DFT+U combination error for any iron-containing oxide compound by using the so-determined slope m 0 of the best-fit line.
In this work (Section 3.2), we compare this method to a more general approach, where the best-fit line is not forced through the origin, i.e., ∆e fit (x Fe ) = m·x Fe + c. It is shown that the line approximately approaches the origin if both binary iron oxides and representative ternary compounds are taken into account. In addition, we tested this method for five different values of U Fe in order to discuss its generality and applicability. A U value is derived, for which the corrected values are, on average, as close as possible to the experimental results (Section 3.3). The oxygen correction derived in Section 3.1 was applied, and the effect of its uncertainty on the DFT/DFT+U correction is discussed.

Derivation of Phase Diagrams
The synthesis process of a material can be imagined as an exchange of elemental components and precursor compounds between reservoirs and the forming phase. Assuming thermodynamic equilibrium, the total energy of a solid phase-say , A l M m O n -can be expressed as the stoichiometric sum of the reservoir energies per atom-namely, their chemical potentials µ. Using Equation (1), this leads, in the case of elemental reservoirs, to a formation energy: Here, ∆µ i := ∆µ i − µ (0) i , i.e., the chemical potentials are referenced to the energies of the elemental ground state phases. This formulation allows to relate the theoretical definition to experimental synthesis conditions, which can be described as "rich" in a component X, if ∆µ X is close to zero, and as "poor" in X, if ∆µ X has a large negative value. In the X-rich case (∆µ X = 0), the elemental phase X forms in its ground state with energy µ (0) X , imposing the constraints ∆µ X < 0 on the chemical potentials of all involved elements for the formation of a single compound phase.
Analogously, to prevent the formation of competing, unwanted phases from a subset of the provided components (e.g., A l O n ), the stoichiometric sum of the corresponding chemical potentials must be lower than the formation energy of the wanted phase (e.g., l ∆µ A + n ∆µ O < E form (A l O n )). The phase diagram of A l M m O n is then defined by the space of chemical potentials fulfilling Equation (3) together with every possible constraint of such a form. For a ternary compound, it can be visualized by a 2-dimensional graph with two of the chemical potentials as axes and the third one sketched by contour lines or color coding. In such a plot, lines represent the constraints confining the stability region.
The chemical potential of a gaseous phase, such as O 2 , can be formally expressed as a function of the temperature and pressure of the gas [35]. If the synthesis takes place under a gas atmosphere, this links the chemical potential directly to accessible experimental processing parameters. However, it is difficult to give a quantitative interpretation of the chemical potentials of substances not provided as a gas during the synthesis beyond statements of conditions being "rich" or "poor" in the respective component.
In order to reflect the solid-state reaction process in the phase diagram, where often precursor compounds are mechanically mixed and heat-treated rather than nongaseous elemental phases, the axes need to be transformed to the chemical potentials of such compound phases and ∆µ O . The formalism is explicitly demonstrated in the appendix (Appendix D) for the example of the synthesis of LaFeO 3 from La 2 O 3 and Fe 2 O 3 . It can be adopted to describe various processing routes for arbitrary compounds.

O 2 Correction
Using the computational DFT settings described in Section 2.1, the total energies were calculated for binary oxide compounds composed of elements from the main groups I to IV of the periodic table: A 2 O with A = Li-Cs, AO with A = Be-Ba, A 2 O 3 with A = Al-Tl, and AO 2 with A = Si-Pb. For those compounds, we chose the most stable structures according to experimental observations at ambient conditions (see Table A1 in Appendix A). We calculated the formation energies according to Equation (1) from the total energies of the corresponding elemental ground-state phases (the chemical potentials µ (0) ). In Figure 1, the results are plotted against experimental values of the standard enthalpies of the formation taken from References [36][37][38][39]. With the exception of Tl 2 O 3 and PbO 2 , the points can be represented by a line, y = x + b, which confirms the generality of the approach of Wang et al. [13], where six arbitrarily chosen compounds were taken into consideration for this analysis.
Excluding Tl 2 O 3 and PbO 2 , a fit of the data points leads to b = 0.64 eV, which can be regarded as an average energy error per O atom in the simulations. The derived b is close to the value given by Wang et al. of 0.68 eV/O [13] and to 0.7 eV/O, which is used in the (MP) database [18]. The correction can now be applied by adding b to the chemical potential of oxygen: Performing the correction in this way reflects that the deviation of the formation energies from the experimental values has its origin in an inaccuracy in µ (0) (O) rather than in E total (A l O n ). However, the same results are obtained for subtracting the product nb from E total and using the unchanged µ (0) (O), as it is done in the MP approach (denoted there as "MP anion correction"). From our results, clear trends of the data points belonging to different groups of the periodic table cannot be identified, except for a clustering of group II compounds at the lower formation energy values (blue squares) and a slight but systematic offset of the group IV compounds from the best-fit line (yellow upward triangles).
Such an offset was reported by Wang et al. for SiO 2 [13] as well and explained by the highly covalent character of the Si-O bond. As in the case of the O 2 molecule, the energy of such a bond is, to some extent, affected by the over-binding tendency of oxygen, leading to partial error cancellation and a formation energy of the compound closer to the experimental value. As apparent from our data, this argument can be applied for GeO 2 and SnO 2 as well. Mixed bonding configurations in the oxides of the heavy elements Pb and Tl could be the reason for an even stronger error cancellation in the case of the formation energies of PbO 2 and Tl 2 O 3 . As an alternative or complementary approach to explain the deviation of these points from the best-fit line, we calculated the energies of the compounds PbO 2 and Tl 2 O 3 and reference phases Pb and Tl by taking into account spin-orbit coupling. This led to energy shifts of the order of 0.1 eV/O in the direction of the best-fit line.
Materials 2020, 13, x FOR PEER REVIEW 6 of 21 uncertainties in the experimental data may contribute to the uncertainty of the correction value as well. Note that the measurements are carried out at ambient conditions, i.e., at room temperature and atmospheric pressure. On the other hand, the calculated formation energies are derived at zero temperature and pressure. An energy contribution ∆ (O ; , ) can, in principle, be added to the ground state chemical potential of the oxygen in Equation (1). It can be formally derived from statistical mechanics considering the entropy effects from the different degrees of freedom of the O2 molecule [35]. However, such a correction must not be additionally included in the method described above, since it is already compensated for by the shift , which encompasses both the O2 overbinding error inherent in DFT, as well as energy contributions due to differing external conditions. Analogously, this argumentation is valid for energy shifts ∆ ( ; , ) of the nonoxygen components and ∆ (A O ; , ) of the compound of interest. In order to quantify separately the effects from O2 over-binding and different external conditions, all of these contributions need to be explicitly calculated, e.g., by considering volume changes with the temperature and applying an equation of state or by performing thermodynamic integrations of heat capacities [15]. An alternative way to deal with the effect of O2 over-binding in DFT was recently reported by Gautam and Carter [40]. They applied the strongly constrained and appropriately normed (SCAN) approximation for the exchange-correlation functional and obtained formation energies of main group oxide compounds, which agreed quite well with the experimental values without the need to perform any post-processing corrections. With our study, though, we aimed to evaluate and discuss common, widely applied DFT methods, which are implemented in many of the available codes and which, for example, were used to generate large sets of data as available in the MP database [18].

Energy Correction for the Combined DFT and DFT+U Results
Figure 2a depicts the differences ∆ defined in Equation (2) between the formation energies per atom derived by using the results from the DFT and DFT+U calculations and the corresponding The mean deviation of all considered data points from the best-fit line is ±0.09 eV/O, and the cited values used by other research groups are within this range as well. This range of b values has to be taken into account in further calculations using the O 2 correction value. Next to numerical fluctuations, which we reduced as much as possible by applying strict convergence criteria, uncertainties in the experimental data may contribute to the uncertainty of the correction value as well.
Note that the measurements are carried out at ambient conditions, i.e., at room temperature and atmospheric pressure. On the other hand, the calculated formation energies are derived at zero temperature and pressure. An energy contribution ∆µ(O 2 ; T, p) can, in principle, be added to the ground state chemical potential of the oxygen in Equation (1). It can be formally derived from statistical mechanics considering the entropy effects from the different degrees of freedom of the O 2 molecule [35]. However, such a correction must not be additionally included in the method described above, since it is already compensated for by the shift b, which encompasses both the O 2 over-binding error inherent in DFT, as well as energy contributions due to differing external conditions. Analogously, this argumentation is valid for energy shifts ∆µ(A; T, p) of the nonoxygen components and ∆E total (A l O n ; T, p) of the compound of interest. In order to quantify separately the effects from O 2 over-binding and different external conditions, all of these contributions need to be explicitly calculated, e.g., by considering volume changes with the temperature and applying an equation of state or by performing thermodynamic integrations of heat capacities [15].
An alternative way to deal with the effect of O 2 over-binding in DFT was recently reported by Gautam and Carter [40]. They applied the strongly constrained and appropriately normed (SCAN) approximation for the exchange-correlation functional and obtained formation energies of main group oxide compounds, which agreed quite well with the experimental values without the need to perform any post-processing corrections. With our study, though, we aimed to evaluate and discuss common, widely applied DFT methods, which are implemented in many of the available codes and which, for example, were used to generate large sets of data as available in the MP database [18].

Energy Correction for the Combined DFT and DFT+U Results
Figure 2a depicts the differences ∆e defined in Equation (2) between the formation energies per atom derived by using the results from the DFT and DFT+U calculations and the corresponding experimental standard enthalpies of the formation of the binary compounds FeO, Fe 2 O 3 , and Fe 3 O 4 and the ternary compounds Li 5 FeO 4 and NaFeO 2 . The phases were set up in their experimentally observed ground-state structures (see Table A1 in Appendix A).
Materials 2020, 13, x FOR PEER REVIEW 7 of 21 experimental standard enthalpies of the formation of the binary compounds FeO, Fe2O3, and Fe3O4 and the ternary compounds Li5FeO4 and NaFeO2. The phases were set up in their experimentally observed ground-state structures (see Table A1 in Appendix A). The ∆ values are plotted against , the fraction of Fe atoms in the compound. The ternary compounds were included to add more data points for the lower values. They were chosen based on the availability of the experimental data and such that they do not contain any further transition metal elements. The data points ∆ were fitted by a straight line for each of the five considered U values for Fe between 3.1 eV and 7.5 eV. A linear fit is well-justified for most of the data points, except for those belonging to Fe3O4 and FeO obtained with the higher U values. The lines ∆ ( ) have slopes and -intercepts , which themselves exhibit linear trends as a function of U: Figure 2b. The given ranges of the -intercepts of ( ) and ( ) The ∆e values are plotted against x Fe , the fraction of Fe atoms in the compound. The ternary compounds were included to add more data points for the lower x Fe values. They were chosen based on the availability of the experimental data and such that they do not contain any further transition metal elements. The data points ∆e were fitted by a straight line The lines ∆e fit U (x Fe ) have slopes m U and y-intercepts c U , which themselves exhibit linear trends as a function of U: m U = m(U) ≈ 0.24·U + (1.02 ∓ 0.04) eV and c U = c(U) ≈ 0.02·U − (0.13 ± 0.04) eV. m(U) is displayed in Figure 2b. The given ranges of the y-intercepts of m(U) and c(U) stem from the uncertainty limits ±∆b in the O 2 over-binding correction described in the previous section. By changing b within these limits, the energy differences ∆e shown in Figure 2a are consistently shifted to higher or lower values between the limits, which, for the sake of clarity, are exemplarily depicted by error bars only for the data points corresponding to U = 3.1 eV. The presented numbers determining ∆e fit U (x Fe ) are, of course, sensitive to the incorporation of additional ternary Fe-containing compounds in the analysis. However, this effect is of minor influence, as long as the data points of additional compounds do not deviate more strongly from the lines than those already considered.
∆e fit U (x Fe ) shifts the calculated formation energy per atom of an arbitrary oxide compound with a fraction of x Fe Fe atoms closer to its experimental value. Following the definition of the formation energy (Equation (1)), this correction can be formally applied by changing the total energies calculated with DFT+U at a given U, according to: As a proof of principle, correcting the energies E total U (LaFeO 3 ) in this manner leads to formation energies very close to the experimental values. This can be seen in Figure 2a, where the open symbols corresponding to LaFeO 3 , which were not included in the fitting processes, closely match the values of ∆e fit U (x Fe = 0.2) for all considered U values.
In order to compare this generalized approach for deriving the correction energy arising in the combined DFT and DFT+U calculations to the method described by Jain et al. [14], we repeated the whole analysis using best-fit lines through the origin, marked by the dashed lines in Figure 2a. For the higher U values, the lines hardly deviate from the lines ∆e fit U (x Fe ). The slopes, which are the only correction parameters in this approach, are derived as m 0 (U) = 0.29·U + (0.70 ∓ 0.13) eV (depicted in the Figure 2b). Their range of uncertainty that again originates in the uncertainty limits ±∆b of the O 2 over-binding correction, is enlarged, as compared to the general method.

Determination of U
The data points ∆e and, therefore, the deviations δ := ∆e − ∆e fit , follow linear trends with U for each of the five iron compounds included in the fitting procedure. Using the corresponding fit parameters, the average of the squared deviations, known as the mean squared error MSE, can be obtained as a function of U: The square root of this expression generally serves as a measure for the quality of the fits and, specifically, here for the quality of the correction scheme. It is displayed in the Figure 2c for both types of fitting discussed in Section 3.2. A difference is visible only for U values below 4.5 eV. Minimization leads to U Fe = (3.95 ± 0.17) eV and U Fe = (4.20 ± 0.37) eV for the general fits and the fits through the origin, respectively. The given ranges depicted by bars originate again from the uncertainty limits of the oxygen correction value. Around these minima, the application of either of the two correction schemes reproduces the experimental formation energies on average with an accuracy of about 0.03 eV per atom. The general fit leads to slightly more accurate results, but this difference can be regarded as insignificant, considering the uncertainties in the approaches.
A U value for the 3d-electrons of Fe of 4.0 ± 0.1 eV was also derived by Wang et al. [13] following a different but related approach. The authors calculated oxidation reaction energies E (r) between the binary iron oxides FeO, Fe 2 O 3 , and Fe 3 O 4 . No elemental Fe phase is involved in these reactions. Therefore, E (r) could be determined by balancing only the uncorrected total energies obtained by DFT+U and the previously corrected chemical potential of oxygen. The change of E (r) with U, which can be understood heuristically (see Appendix B), led to an apparently optimal U value where the results for E (r) best matched the experimental reaction energies. E (r) depends linearly on the formation energies of the involved compounds. Since these values also enter the method described in this section, the approach of Wang et al. may be interpreted as being formally equivalent to the method described above if only the binary oxides were included in the fitting. This is, however, not the case, as is shown in detail in Appendix C.
A U Fe value of approximately 4 eV is widely reported in the literature to lead to band gaps and reaction energies of iron oxides in good agreement with experimentally derived values [13,[22][23][24][25][26]. It has to be noted that, in all of the cited studies, as in this work, the DFT code VASP was used together with the standard set of PAW pseudopotentials provided in the VASP package. In contrast, Xu et al. [41] reported optimal values of U Fe for the iron oxides, which were derived with the DFT code Quantum ESPRESSO (QE) [42] by the linear response method [43]. With U Fe between 3.47 eV and 4.10 eV, the values they obtained using pseudopotentials from the original QE library were in acceptable agreement with the VASP values. However, considerably different U Fe values between 5.21 eV and 6.07 eV resulted from the application of QE with a different set of pseudopotentials [41], namely those from the GBRV high-throughput library [44]. This shows that one has to be careful in proposing and adopting universal element dependent U values, since they can strongly depend on the choice of the pseudopotential and, additionally, may differ between DFT codes due to different implementations of the Hubbard-U correction.

Band Gaps of the Iron Compounds
The usual approximations of DFT, LDA, or GGA insufficiently describe the correlated electrons in 3d-orbitals of transition metals, which, together with the oxygen 2p-orbitals, form the valence band edge in transition metal oxide compounds. This leads to incorrect, generally underestimated, band gaps, which can be opened up by applying U to localize the 3d-orbitals, making the electronic structure less "metal-like". Accordingly, an optimal U value can be found by tuning U until the experimental band gap value is reproduced in the band-structure calculations [22]. To compare this approach to the one described in Section 3.3, we calculated the band gaps of the considered compounds for the different U values of iron. As shown in Figure 3, the band gaps exhibit the expected opening with increasing U Fe , following, in most cases, almost linear trends. For Fe 3 O 4 , a distinct bend of the line around U ≈ 5.3 eV is apparent. This is a significant feature that is visible, albeit much less pronounced, for the other compounds as well. Very strict convergence criteria were applied to ensure that this effect is not an artifact of numerical fluctuations but, rather, a characteristic feature of the DFT+U approach. On a small scale, such a bending was also observed in the total energies and the quantities deduced from them. This can be seen, for example, by closely examining the deviations of the data points from the lines in Figure 2b. However, this is considered as insignificant in the approximate analysis conducted there.
The following experimental band gap energy data were taken from the literature: for  Figure 3. Fe 3 O 4 is known to be electronically conductive in the cubic state, so there should be no band gap, which is reproduced by the DFT+U calculations for U below about 3 eV. A small band gap of 0.07 eV is reported only for the monoclinically distorted structure forming below the Verwey transition temperature [50]. The complex electronic structure of Fe 3 O 4 of mixed octahedrally coordinated Fe 3+ and tetrahedrally coordinated Fe 2+ ions cannot be sufficiently described by the same, single U value that is applied for the octahedrally coordinated Fe ions in the same trivalent charge state in the other compounds. For Li 5 FeO 4 and NaFeO 2 , no experimental band gap values are reported, to the best of our knowledge. Hoang et al. [51] derived a band gap of 4.4 eV for Li 5 FeO 4 using hybrid-functional DFT calculations with nonoptimized mixing or screening parameters. This value is considerably larger than the DFT+U results obtained here. For NaFeO 2 , only one DFT+U calculation was found in the literature [52] with a band gap value of 1.5 eV for U Fe = 4 eV, in agreement with our results.

Phase Diagrams of LaFeO3, Li5FeO4, and NaFeO2
In order to derive the phase diagrams for the ternary oxide compounds LaFeO3, Li5FeO4, and NaFeO2, we calculated ground-state formation energies for a set of binary oxides and ternary iron oxides containing the elements La, Li, or Na, respectively. A U value for Fe of 4 eV was applied, and the corrections determined in Sections 3.2 and 3.3 were adopted. The choice of compounds was based on phases listed in the MP database [18], with energies close to the convex hull. The way the phase diagrams are presented in the following reflects the solid-state processing routes that were reported for the considered materials.

LaFeO3
LaFeO3 was experimentally synthesized by a mechanical mixing treatment of La2O3 and Fe2O3 powders and calcination under an oxygen atmosphere [53,54], according to: The phase diagram with variables ∆ (La O ) (∆ (Fe O )) and ∆ (O) is shown in Figure 4, with an alternative y-axis assigned to the pressure of the oxygen gas at a temperature of 1400 K. In addition to the binary iron oxides and La2O3 (space group 3 ), the phase La3FeO6 ( 2 ) was taken into account. However, it turned out that the latter is not relevant for the phase diagram of LaFeO3, since the corresponding phase separation line lies outside of the displayed region.
Since the stability region is horizontally confined by the precursor compounds of the reaction,

Phase Diagrams of LaFeO 3 , Li 5 FeO 4 , and NaFeO 2
In order to derive the phase diagrams for the ternary oxide compounds LaFeO 3 , Li 5 FeO 4 , and NaFeO 2 , we calculated ground-state formation energies for a set of binary oxides and ternary iron oxides containing the elements La, Li, or Na, respectively. A U value for Fe of 4 eV was applied, and the corrections determined in Sections 3.2 and 3.3 were adopted. The choice of compounds was based on phases listed in the MP database [18], with energies close to the convex hull. The way the phase diagrams are presented in the following reflects the solid-state processing routes that were reported for the considered materials.

LaFeO 3
LaFeO 3 was experimentally synthesized by a mechanical mixing treatment of La 2 O 3 and Fe 2 O 3 powders and calcination under an oxygen atmosphere [53,54], according to: The phase diagram with variables ∆µ(La 2 O 3 ) (∆µ(Fe 2 O 3 )) and ∆µ(O) is shown in Figure 4, with an alternative y-axis assigned to the pressure of the oxygen gas at a temperature of 1400 K. In addition to the binary iron oxides and La 2 O 3 (space group Ia3), the phase La 3 FeO 6 (Cmc2 1 ) was taken into account. However, it turned out that the latter is not relevant for the phase diagram of LaFeO 3 , since the corresponding phase separation line lies outside of the displayed region.
involving these elements. Heifets et al. [55] calculated and discussed the phase diagrams of LaFeO3 with respect to the elemental chemical potentials Δ and Δ . The phase diagram that is derived there using the experimental input data shows chemical potentials Δ between −3.2 and −2.8 eV along the Δ = 0 line in agreement with our results in Figure 4. This is expected, since the calculated formation energies we used for deriving the phase diagram are close to the experimental values because of the applied correction procedure.
Although a representation of the phase diagram in the form of Figure 4 is closer to the reality of the actual synthesis process than a representation with respect to elemental reservoir energies [56][57][58][59], it remains a difficulty of adjusting and quantifying or, in general, interpreting Δ and Δ in an experimental setup. In analogy to a gas phase, "rich" conditions correspond to a high chemical reactivity, which can, for example, be realized by a high density of reactive surfaces, i.e., small and densely packed grains.

Li5FeO4
The experimental synthesis of Li5FeO4 was achieved by tempering a powder mixture of Li2O and elemental Fe at 1000 °C and a low pressure of about 1 mPa [60], following the reaction: The corresponding phase diagram is shown in Figure 5, with variables ∆ (Li O) and (O ) for = 1300 K. Since Fe is involved in the reaction, its chemical potential is of relevance as well. This energy of the second reaction partner cannot be represented by an alternative horizontal axis, as in the case of LaFeO3 (see Figure 4), due to the existence of an additional species (Li) in the reaction (see detailed explanation in Appendix D). Instead, ∆μ(Fe) was made visible via color coding in the stability region. Next to Li2O (space group 3 ) and the binary iron oxides, the phases taken into account were Li2O2 ( 6 ), Li2FeO2 ( ), Li2FeO3 ( 2), and LiFeO2 ( 3 ). Based on the scale of the horizontal axis, Li5FeO4 can be expected to form only under very rich Li2O conditions. The experimentally applied pressure of 10 atm corresponds to rather rich Fe Since the stability region is horizontally confined by the precursor compounds of the reaction, the formation of an intermediate phase during the synthesis of LaFeO 3 is not expected. Molecular oxygen and metallic Fe can form at very high and at very low gas pressures, respectively. Lines corresponding to metallic La, FeO, or Fe 3 O 4 do not cross the stability region, indicating that no thermodynamic condition can be realized under which LaFeO 3 would be in equilibrium with either of these phases. The functional dependencies between ∆µ La , ∆µ Fe , ∆µ O , ∆µ La 2 O 3 , and ∆µ Fe 2 O 3 (see Appendix D) can be applied to determine the chemical potentials of the metallic elements at each point in the phase diagram, which are, e.g., needed for calculating the point defect formation energies involving these elements. Heifets et al. [55] calculated and discussed the phase diagrams of LaFeO 3 with respect to the elemental chemical potentials ∆µ La and ∆µ Fe . The phase diagram that is derived there using the experimental input data shows chemical potentials ∆µ O between −3.2 and −2.8 eV along the ∆µ Fe = 0 line in agreement with our results in Figure 4. This is expected, since the calculated formation energies we used for deriving the phase diagram are close to the experimental values because of the applied correction procedure.
Although a representation of the phase diagram in the form of Figure 4 is closer to the reality of the actual synthesis process than a representation with respect to elemental reservoir energies [56][57][58][59], it remains a difficulty of adjusting and quantifying or, in general, interpreting ∆µ La 2 O 3 and ∆µ Fe 2 O 3 in an experimental setup. In analogy to a gas phase, "rich" conditions correspond to a high chemical reactivity, which can, for example, be realized by a high density of reactive surfaces, i.e., small and densely packed grains.

Li 5 FeO 4
The experimental synthesis of Li 5 FeO 4 was achieved by tempering a powder mixture of Li 2 O and elemental Fe at 1000 • C and a low pressure of about 1 mPa [60], following the reaction: The corresponding phase diagram is shown in Figure 5, with variables ∆µ(Li 2 O) and p(O 2 ) for T = 1300 K. Since Fe is involved in the reaction, its chemical potential is of relevance as well. This energy of the second reaction partner cannot be represented by an alternative horizontal axis, as in the case of LaFeO 3 (see Figure 4), due to the existence of an additional species (Li) in the reaction (see detailed explanation in Appendix D). Instead, ∆µ(Fe) was made visible via color coding in the stability region. Next to Li 2 O (space group Fm3m) and the binary iron oxides, the phases taken into account were Li 2 O 2 (P6 3 mmc), Li 2 FeO 2 (Immm), Li 2 FeO 3 (C2), and LiFeO 2 (Fd3m).
Materials 2020, 13, x FOR PEER REVIEW 12 of 21 conditions as well. During the vaporization of Li5FeO4, the precipitation of solid LiFeO2 was observed [60], accompanied by gaseous Li and O2. This is in agreement with the phase diagram, where LiFeO2 terminates the narrow stability region on one side. The Li gas cannot solidify, since the corresponding line (∆ (Li) = 0) lies far outside the Li5FeO4 region at unreasonably low oxygen gas pressures below 10 atm. In contrast to the situation of LaFeO3, the phase diagram indicates no equilibrium between Li5FeO4 and a binary iron oxide. For example, the line ∆μ(Fe O ) = 0, which is parallel to ∆μ(LiFeO ) = 0, is located at ∆μ(Li O) = −0.18 eV outside of the displayed region. This impedes a direct solid-state route, such as Li2O + Fe2O3  2 Li5FeO4, which, at rich Fe2O3 conditions, would lead to the formation of LiFeO2 instead.

NaFeO2
NaFeO2 was experimentally synthesized by milling the precursor compounds Na2O2 (sodium peroxide) and Fe3O4 and exposing the mixture to a temperature of about 900 K [27,28]. Due to stoichiometric constraints, oxygen has to be released during this reaction: 3 Na2O2 + 2 Fe3O4 ⟶ 6 NaFeO2 + O2.
The phase diagram is presented with variables ∆ (Na O ), ∆ (O), and (O ) for (O ) = 900 K in Figure 6, with the chemical potential of Fe3O4 being depicted by a color code. Next to Na2O2 (space group 6 2 ) and the binary iron oxides, the phases taken into account were Na2O ( 3 ), NaO2 ( ), Na2FeO3 ( 1 ), Na3FeO3 ( 2 / ), Na4FeO4 ( 1 ), and NaFeO2 ( 3 ). Further phases consisting of Na, Fe, and O turned out to be irrelevant for the stability discussion of NaFeO2.
Under Fe3O4-rich conditions, NaFeO2 only forms at low pressures and Na2O2-poor conditions (lower left part of the stability region). If the atmospheric pressure is too high, the oxidation of Fe3O4 to Fe2O3 can occur, and if it is too low, a reduction to metallic Fe is expected. Under Na2O2-rich

NaFeO 2
NaFeO 2 was experimentally synthesized by milling the precursor compounds Na 2 O 2 (sodium peroxide) and Fe 3 O 4 and exposing the mixture to a temperature of about 900 K [27,28]. Due to stoichiometric constraints, oxygen has to be released during this reaction: The phase diagram is presented with variables ∆µ(Na 2 O 2 ), ∆µ(O), and p(O 2 ) for T(O 2 ) = 900 K in Figure 6, with the chemical potential of Fe 3 O 4 being depicted by a color code. Next to Na 2 O 2 (space group P62m) and the binary iron oxides, the phases taken into account were Na 2 O (Fm3m), NaO 2 (Pnnm), Na 2 FeO 3 (P1), Na 3 FeO 3 (P2 1 /c), Na 4 FeO 4 (P1), and NaFeO 2 (R3m). Further phases consisting of Na, Fe, and O turned out to be irrelevant for the stability discussion of NaFeO 2 .
Materials 2020, 13, x FOR PEER REVIEW 13 of 21 diagram, two alternatively plausible processing routes can be proposed: Na3FeO3 + Fe2O3  3 NaFeO2 and Na2O + Fe2O3  2 NaFeO2. The latter does not include the complication of preliminarily having to synthesize Na3FeO3, but it can only be realized under not-too-rich Na2O conditions. A reaction of metallic Fe with sodium superoxide (NaO2) appears possible as well, but due to the high reactivity of NaO2 with water and a costly production process, this route is presumably inexpedient from a practical point of view.

Summary and Conclusions
In the first part of this work, we compiled, reconsidered, and reevaluated a set of correction methods previously described in the literature and applied in a frequently accessed materials database (Materials Project, MP), which deals with inaccuracies of the usual LDA or GGA calculations of DFT in deriving accurate formation energies of transition metal oxide compounds. Common to these methods is the incorporation of physically justified and generally applicable parameters for systematic error corrections, which can be tuned to minimize an average deviation of computed and experimentally determined energies for a variety of compounds. We come to the following conclusions: • The method to correct the formation energy error due to the over-binding of the O2 molecule remains valid if a larger set of non-transition metal oxide compounds from the groups I to IV of the periodic table is considered instead of the previously chosen smaller subset. The magnitude Under Fe 3 O 4 -rich conditions, NaFeO 2 only forms at low pressures and Na 2 O 2 -poor conditions (lower left part of the stability region). If the atmospheric pressure is too high, the oxidation of Fe 3 O 4 to Fe 2 O 3 can occur, and if it is too low, a reduction to metallic Fe is expected. Under Na 2 O 2 -rich conditions (upper right part of the stability region), NaFeO 2 only forms at high oxygen pressures and Fe 3 O 4 -poor conditions, and the phases NaO 2 and Na 4 FeO 4 are likely to emerge. Based on the phase diagram, two alternatively plausible processing routes can be proposed: Na 3 FeO 3 + Fe 2 O 3 → 3 NaFeO 2 and Na 2 O + Fe 2 O 3 → 2 NaFeO 2 . The latter does not include the complication of preliminarily having to synthesize Na 3 FeO 3 , but it can only be realized under not-too-rich Na 2 O conditions. A reaction of metallic Fe with sodium superoxide (NaO 2 ) appears possible as well, but due to the high reactivity of NaO 2 with water and a costly production process, this route is presumably inexpedient from a practical point of view.

Summary and Conclusions
In the first part of this work, we compiled, reconsidered, and reevaluated a set of correction methods previously described in the literature and applied in a frequently accessed materials database (Materials Project, MP), which deals with inaccuracies of the usual LDA or GGA calculations of DFT in deriving accurate formation energies of transition metal oxide compounds. Common to these methods is the incorporation of physically justified and generally applicable parameters for systematic error corrections, which can be tuned to minimize an average deviation of computed and experimentally determined energies for a variety of compounds. We come to the following conclusions: • The method to correct the formation energy error due to the over-binding of the O 2 molecule remains valid if a larger set of non-transition metal oxide compounds from the groups I to IV of the periodic table is considered instead of the previously chosen smaller subset. The magnitude of the energy correction we derived (0.64 eV) agrees within 0.1 eV with the reported values, which reflects the uncertainty range we determined for the approach. Experimental enthalpies of the formations for the considered compounds can be reasonably reproduced within 0% to 5%, except for Tl 2 O 3 and PbO 2 , for which the values deviate more. A possible explanation is given.

•
For the binary iron oxide compounds, we confirmed that it is well-justified to correct the error arising in combining the results of DFT and DFT+U calculations by adding a value proportional to the Fe content to the formation energy. While it was originally derived considering only the binary oxides for a fixed U value for Fe, we strengthened and generalized the scheme by (a) taking into account ternary compounds, (b) not a priori constraining the correction value to be zero for hypothetical compounds without Fe, and (c) considering different values of U.

•
Our U-dependent correction value offers a new possibility to determine an optimal U value for Fe for which the experimental formation energies are reproduced best. With this approach, we confirmed the frequently used value of U ≈ 4 eV, which additionally turned out to reproduce experimental band gaps of the considered compounds within 0.3 eV.
In the second part of this work, by applying the above-mentioned correction schemes, we calculated the formation energies for a set of binary oxide and ternary iron oxide compounds containing La, Li, or Na in order to derive the phase diagrams for LaFeO 3 , Li 5 FeO 4 , and NaFeO 2 with respect to the reservoir energies (chemical potentials) of elements and/or precursor compounds. Our representation of the phase diagrams corresponds closely to actual solid-state synthesis routes reported in the experimental literature. This allows for motivating the experimental adjustments of pressures, for predicting the phase transformations upon changing the synthesis conditions, and for explaining the phases found after vaporization, which are demonstrated for the described compounds. In addition, alternative synthesis routes can be proposed and compared.

Appendix A : List of Calculated and Experimental Formation Energies
The corresponding reaction energy per formula unit can be expressed as: The application of U changes the total energies. We define and analogously δ Fe 2+ (for FeO) as the differences of the total energies per Fe atom in the compounds for a fixed difference between a higher (U > ) and a lower U value (U < ). By applying higher U values, the correlated electrons, which become too delocalized by LDA or GGA, are forced to a higher localization, which increases the energy of the system due to Coulomb repulsion. Therefore, δ Fe 3+ > 0 and δ Fe 2+ > 0. With higher oxidation states, the number of correlated electrons is increased, so U has a stronger influence on the change in the energy: δ Fe 3+ > δ Fe 2+ . A comparison of the reaction energies for different U values then leads to: This implies, that oxidation reaction energies increase with the increasing U, which is confirmed by our results for iron oxides and by the results of Wang et al. [13] for vanadium, chromium, manganese, iron, cobalt, nickel, and copper oxides.

Appendix C : Comparison of Methods for the Determination of U
In this section, two conceptually different approaches are explicitly described to derive an optimal U value by comparing the experimental and calculated formation and oxidation energies of iron oxides. The first method is equivalent to the procedure worked out in Section 3.2, except that now only the energies of the binary iron oxide phases are taken into account. The function to be minimized with respect to U is: with ∆e as defined in Equation (2). For the fitting, the linear functions ∆e fit,0 through the origin are chosen, since the difference between ∆e fit,0 and the general linear fit ∆e fit is insignificant in the U region where f is minimal (see Figure 2a). The sum runs over Fe 2 O 3 , Fe 3 O 4 , and FeO. With x Fe ≡ x, the fit lines can be expressed as: By expressing the energy differences as linear functions of U: ∆e i (U) = a i U + b i , only the parameters x i , a i , and b i enter the function f (U). Minimization leads to: with α := i x i a i /x 2 i and β := i x i b i /x 2 i . In the approach by Wang et al. [13], the energies of the oxidation reactions (r1) 2 if i = 1, 2, 3 stands for Fe 2 O 3 , Fe 3 O 4 , and FeO, respectively. Since this relation holds for the calculated, as well as for the experimental, values, it can be formulated for their respective differences ∆e (r1,r2) and ∆e 1,2,3 . Note that changing the energies e 1,2,3 by any correction proportional to the number of Fe atoms involved cancels out and has no effect on e (r1,r2) , so it is irrelevant whether or not the correction of the combined DFT and DFT+U data described in Section 3.2 is considered here. Therefore, we can apply the linear dependencies ∆e i (U) = a i U + b i with the same set of parameters as above in the first approach, which leads to Equation (A7).
Materials 2020, 13, x FOR PEER REVIEW  17 of 21 and ∆ , , . Note that changing the energies , , by any correction proportional to the number of Fe atoms involved cancels out and has no effect on ( , ) , so it is irrelevant whether or not the correction of the combined DFT and DFT+U data described in Section 3.2 is considered here. Therefore, we can apply the linear dependencies ∆ ( ) = + with the same set of parameters as above in the first approach, which leads to Equation (A6). The function to be minimized to find the closest match between the calculated and experimental oxidation energies is: ( ) = ∆ ( ) ( ) + ∆ ( ) ( ) . (A9) With the definitions ⃗ ≔ ( , , ) and ⃗ ≔ ( , , ), the U value where is minimal can be found as: with the matrix as defined in Equation (A7). From the data points ∆ ( ) shown in Figure 2a, one obtains ⃗ = (0.132, 0.153, 0.105) eV/U and ⃗ = (0.198, 0.209, 0.523) eV, leading to = 4.00 eV and = 4.05 eV for the first (Equation (A6)) and second approach (Equation (A9)), respectively. The functions ( ) and ( ) are shown in Figure A1, indicating that the two approaches are indeed different from each other, even though they use the same set of calculated and experimental input data and lead to very close minima.

Appendix D: Phase Diagram with Respect to the Precursor Compounds
This section explicitly describes for the example of LaFeO3 (LFO), how a phase diagram can be derived with respect to reservoir energies of binary compounds and molecular oxygen instead of the usual representation considering only the elemental phases. The latter corresponds to a hypothetical solid-state reaction La + Fe + (3/2)O2  LaFeO3. In a more realistic synthesis scenario, stable oxide compounds for LFO, e.g., La2O3 and Fe2O3, are mixed by mechanical and/or chemical treatment and exposed to an oxygen atmosphere [53,54]. In an analogy to Equations (1)  The function to be minimized to find the closest match between the calculated and experimental oxidation energies is: (A9) With the definitions → a :=(a 1 , a 2 , a 3 ) and → b :=(b 1 , b 2 , b 3 ), the U value where g is minimal can be found as: with the matrix T as defined in Equation (A8). From the data points ∆e i (U) shown in Figure 2a Figure A1, indicating that the two approaches are indeed different from each other, even though they use the same set of calculated and experimental input data and lead to very close minima.
according to the reaction 4 Li 2 O + Fe −→ Li 5 FeO 4 + 3 Li . The procedure then follows the steps described above. However, in contrast to the situation for LaFeO 3 and Equation (A11), Equation (A16) contains two free variables. Therefore, choosing ∆µ(Li 2 O) as the horizontal (x)-axis of the phase diagram does not allow for an alternative horizontal axis x = ∆µ(Fe) or x = ∆µ(Li). Instead, these quantities can be represented by height functions in the phase diagram, and, e.g., illustrated by color coding, as shown in Figures 5 and 6.