Michael John
Craig
*ab,
Felix
Kleuker
a,
Michal
Bajdich
*b and
Max
García-Melchor
*a
aSchool of Chemistry, CRANN and AMBER Research Centres, Trinity College Dublin, College Green, Dublin 2, Ireland. E-mail: craigmi@tcd.ie; garciamm@tcd.ie
bSUNCAT Center for Interface Science and Catalysis, SLAC National Accelerator Laboratory, Menlo Park, California 94025, USA. E-mail: bajdich@slac.stanford.edu
First published on 8th May 2023
Herein we report a method to extract formation energies from oxidation states, which we call FEFOS. This new scheme predicts the formation energies of binary oxides through analyzing unary oxide formation energies as a function of their oxidation states. Taking averages of fitted quadratic equations that represent how elements respond to oxidation and reduction, the weights of these averages are determined by constraining the compound to be neutral. The application of FEFOS results in mean absolute errors of ca. 0.10 eV per atom when tested against Materials Project data for oxides with general formulas A1−zBzO, A1−zBzO1.5, and A1−zBzO2 with specific coordinations. Our FEFOS method not only allows for the prediction of binary oxide formation energies with low variance and high interpretability, but also compares well with state-of-the-art deep learning methods without being biased by training data and the need for large resources to compute it. Finally, we discuss the potential applications of the FEFOS method in tackling the problem of inverse catalyst design.
The notion of an oxidation state (OS) dates back to 1835 when Wohler described elements which can form different oxides.7 It has since seen widespread use as a useful chemical concept to describe chemical bonding and elemental charge in compounds. How this quantity relates to other chemical notions such as covalency or ionicity remains controversial,8–10 even so, the notion of the formal oxidation state is a powerful and prominent concept in chemistry.
In this work, we develop a simple scheme to derive formation energies from oxidation states (FEFOS) to forecast standard formation energies of binary oxides from the Materials Project (MP) database,11 and compare our predictions to the Roost method. This approach takes a well-known phenomenon in solid oxide chemistry – namely, that cations in binary oxides balance their respective charges to their preferred oxidation states – and uses this knowledge to derive a model for the enthalpy of formation of binary oxides based on the enthalpy of formation of unary oxides. Notably, the FEFOS method only requires two parabolas per element and oxide stoichiometry, leading to an interpretable, extensible method with comparable errors to those shown by the state-of-the-art DL approach implemented in Roost. The physical implication of this is that a primary driving force towards the formation of binary oxides lies in the balancing of charges. Furthermore, the interpretability of FEFOS is demonstrated by describing two distinct cases of alloying SnO2 and NbO2 with earth-abundant elements and relating the performance of our model in terms of the preference of specific elements for a specific OS. Lastly, we present how the FEFOS model could be used to enable inverse design in catalysis, specifically for the oxygen evolution reaction (OER). We envision this new method to inform experiments on the propensity of metal oxides towards mixing, the stability of the resulting oxide(s), as well as the optimal concentration of the constituent elements, thus accelerating inverse catalyst design.
For each element, M, and each oxidation state of interest 2y, we construct two quadratic equations, fox,MOy(x) and fred,MOy(x), which take the general form:
fox/red,MOy(x) = a1,MOyx2 + a2,MOyx | (1) |
In Fig. 1a we fit the above quadratic equations over the standard formation energies for the oxidation of Mo, fox,MoO2(x) = 0.022x2 + 0.003x, and the reduction of Cr, fred,CrO2(x) = 0.319x2 − 0.638x, where we use CrO2 and MoO2 as our (0, 0) point and fit through the points which we define as the convex hull in the direction of oxidation or reduction. Given Cr2O3 is the chromium-based oxide with the lowest formation energy, this curve is a mathematical representation of the fact that Cr′s preferred OS is 3+ and is symmetric about this point so that the y-axis is zero (i.e. the same as CrO2) for OS 2+. On the other hand, the fox,MoO2(x) parabola is very shallow and concave, having been fitted through Mo8O23 and MoO3.
![]() | ||
Fig. 1 a) Parabolas for the reduction of Cr and the oxidation of Mo formed by fitting through ground state oxides found in MP, with fitted datapoints indicated by the stoichiometry of the material. The x-axis represents the absolute OS change with respect to the reference MO2 oxide (|x|). Recall that for an oxide AaOc, we have that ![]() |
The general idea is to combine the above quadratic equations to represent binary formation energies in the restricted composition space of Mo1−zCrzO2 with {z ∈ R| 0 ≤ z ≤ 1}. To maintain charge balance, we constrain the average OS of any material AaBbOc to be . We then enforce that the mixed compositions have the same fixed average OS as the initial pure oxides, e.g. in the case of MO2-type oxides this OS is 4+. If we now use o and r to represent the oxidation and reduction of a material A1−zBzOy, (where y = {1, 1.5, 2}) to keep the OS constant when mixing we need:
o(1 − z) − rz = 0, | (2) |
For instance, a Sb0.67Ni0.33O2 could exhibit Sb5+ and Ni2+ under this scheme, so that the overall OS remains 4+. In real systems, however, such ideal oxidation mixing eventually breaks down, and further compensation occurs via formation of defects, i.e. oxygen and metal vacancies.18 Yet, as we demonstrate in the following, our scheme is sufficient to forecast binary oxide formation energies at fixed stoichiometries.
To mathematically define a combination energy, ΔC, we use the quadratic equations fox,AOy(x), fred,AOy(x), fox,BOy(x), and fred,BOy(x) for a hypothetical material A1−zBzOy. These are used to test both combinations of oxidizing and reducing either of the constituent elements, and by minimizing these curves we arrive to:
ΔC = minx(fox,AOy(xz) + fred,BOy(x(1 − z)), |
fred,AOy(x(1 − z)) + fox,BOy(xz)) | (3) |
Then, to represent a mixing formation energy with respect to segregation against CrO2 and MoO2, we define ΔD(A1−zBzOy) with:
ΔD(A1−zBzOy) = 0 for z = {0, 1} | (4) |
ΔD(A1−zBzOy) = ΔC(A1−zBzOy) − (1 − z)ΔC(A1B0Oy) − zΔC(A0B1Oy) | (5) |
ΔEf(A1−zBzOy) ∼ ΔD(A1−zBzOy) + (1 − z)ΔEf(AOy) + zΔEf(BOy) | (6) |
We demonstrate the application of this method to predict ΔD(A1−zBzOy) for mixing Cr and Mo in the MO2 stoichiometry, as seen in the convex hull phase diagram constrained to the Mo1−zCrzO2 phase space in Fig. 1d. Importantly, we observe almost no deviation between our prediction and the observed Mo0.5Cr0.5O2 material in the MP database (ID mp-1213733). Therefore, our FEFOS method holds promise as a rapid means of garnering a forecast for the minimum energy over a span of concentrations for a given stoichiometry, as well as the OS of the respective elements. In the case outlined in Fig. 1, we can read off the predicted OS of Cr and Mo by noting the value of x which minimizes the curves in Fig. 1b and multiplying that value by the concentration of the other element. For example, we see that for z = 0.5, the minimum weighted average in eqn (3) is for xtot. = 1.86. Multiplying this change by z gives that the change in OS should be 0.93 relative to the MO2 unary oxide. Therefore Cr and Mo should have OS values of 3.07 and 4.93, respectively; such non-integral values for OS have been previously proposed to be observable.19
We now gather quadratic equations for elements across the periodic table and compare our approximated ΔEf to the values reported by MP for any binary oxide which is composed of elements within the considered set and with an M:
O ratio of 1
:
2, 1
:
1.5 and 1
:
1. The method used to construct these parabolas and their coefficients, as in eqn (1), for each considered element can be found in the ESI.†
Our FEFOS model is constructed with the notion of a constant average OS across the composition of a binary oxide, so we use a filtering condition which requires: i) A and B are within 2.5 Å to six oxygen atoms, and ii) oxygen atoms are taken to be O2− anions which are within 2.5 Å to 3, 4, or 6 metal atoms, if the binary oxide to be represented is of type A1−zBzO2, A1−zBzO1.5, or A1−zBzO, respectively. If each A and B is octahedral, then the total anionic charge is 12, the overall charge is then 4, 3 or 2 depending on whether there are 3, 4 or 6 metal atoms bound to the A and B atoms. This way we attempt to constrain the overall OS across the elements A and B, where both are elements that exhibit a unique OS for every site in the oxide. The upper bound of 2.5 Å is based on prior studies of metal–oxygen bond length distributions.20 We then gather every MP entry which fit these criteria and compare FEFOS forecasted ΔD(A1−zBzO(1,1.5,2)) values against the MP value.
A parity plot of our predictions for ΔEf(A1−zBzO1.5) against the observed MP is shown in Fig. 2a, where we can see a low error between forecasted and MP-calculated points despite the simplicity of the FEFOS model. MAE values of forecasted and MP-reported formation energies of 0.105 and 0.171 eV per atom are shown for the case where we do and do not filter by coordination, respectively. To highlight the need to filter by coordination, we include the unfiltered cases as pink circles, where we observe a higher error than if we restrict our comparisons to specific coordinations (blue circles). This is expected, since controlling for coordinations should reduce the possibility that each element has distinct OSs, which our model does not represent as it only considers that A and B have constant OSs. To confirm that the remaining structures still represent a diverse set of oxides, in Fig. S1† we have analyzed the distribution of space group numbers in the filtered and unfiltered groups, where we note that while the span of space groups is lower than the unfiltered group, it is not concentrated within any one subset of space groups.
Fig. 2b illustrates the result of including the forecasted ΔEf values against the MP entries where the filtering condition on oxide stoichiometry is applied for each OS. These data include the blue datapoints for A1−zBzO1.5, shown in Fig. 2a, as well as red and green circles for A1−zBzO2 and A1−zBzO oxides, respectively. Here we see comparable MAE values as seen for A1−zBzO1.5 oxides of 0.095 and 0.124 eV per atom for A1−zBzO2 and A1−zBzO oxides, respectively. This underscores the versatility of the FEFOS method described in this work as it can be applied to distinct classes of binary oxides with errors comparable to the best available DL methods. To demonstrate this, in the following we compare the performance of FEFOS to the aforementioned Roost method used to approximate formation energies, which reported a test set of 25663 entries from the Open Quantum Materials Database21 with a MAE of 0.03 eV per atom when trained with 230
959 datapoints.4
The developers of Roost thankfully accompanied their report with a spreadsheet of values including the composition, target variable and prediction for each of the 25663 test datapoints, and therefore we can condition the composition so that we only compare values for A1−zBzO(1,1.5,2) oxides. They also provided a learning curve, where they reported the error in their test set when the DL algorithm had been trained on progressively more datapoints. In Fig. 3, we compare both the test MAE across all data to ensure we reproduce the reported value, and the MAE for all oxides of the same type considered by our FEFOS model, showing these values as a function of the number of training datapoints. Notably, the Roost error is higher for oxides we have studied in this work, as may be expected given the higher degree of complexity within this class of materials. In addition, while the error with Roost becomes lower than with FEFOS when the DL algorithm is trained on more than 100
000 datapoints from the OQMD,21 our model has far lower complexity and more interpretability.
![]() | ||
Fig. 3 Learning curves of the Roost method4 for both the as-reported datapoints (red circles), which include the entire OQMD formation energy test set, and datapoints for OQMD entries which have the form A1−zBzO(1,1.5,2) (blue circles). The overall MAE from our FEFOS method from Fig. 2 is shown as a horizontal dashed line. |
![]() | ||
Fig. 4 Constrained convex hull diagrams showing the predicted ΔD from eqn (5) across the compositions spaces a) Sn1−zBzO2 and b) Nb1−zBzO2 after alloying them with earth-abundant elements, B. The datapoints denote MP-calculated values which the FEFOS model attempts to capture. Predictions are plotted across composition space, not only comparing them with MP data, as in Fig. 2. |
The elements which show the lowest ΔD(Sn1−zBzO2) at low concentrations (z < 0.1) are Nb and W, which have their lowest formation energies for OSs 5+ and 6+, respectively. We attribute this to the fact that all of the tested 3d elements except Ti are most stable at lower oxygen concentrations than is seen in rutile structures (M:
O = 1
:
2), so their fred,MO2(x) equation must out-compete fox,SnO2(x) = 0.77x. On the other hand, for Nb and W the relevant Sn equation is fred,SnO2(x) = 0.057x2 + 0.12x which exhibits a shallower slope. SnO2 is known to form oxygen vacancies more readily than tin vacancies,22 so this mathematical construction agrees with that observation, as well as a recent paper reporting the doping of SnO2 with Nb, W and Ta for application as a transparent conducting oxide.18 As an aside, we highlight that if ΔD(Sn1−zBzO2) < 0 it implies that the material is stable with respect to segregation to SnO2 and BO2, but not necessarily against segregation to other competing phases.
Unlike SnO2, Fig. 4b shows how doping NbO2 with first-row transition metals result in negative ΔD(Nb1−zBzO2) values due to their preference for reduction with respect to MO2 (excluding Ti). This is accompanied by the lowest formation energy for Nb being seen for Nb2O5, implying that Nb prefers to be oxidized from 4+ to 5+. It is also worth noting that, except for Mn, the positions of the minima of our curves appear to align well with the MP-calculated values (circles), which we could reasonably assume are close to the true minima. For example, the Cu, Co, and Ni points are clustered around the Nb (z = 0.33) region, attributable to these metals' affinity for OS 2+, while Fe and Cr have minima closer to Nb (z = 0.5) due to the ‘preference’ of these metals for OS 3+. In contrast, the second and third-row transition metals Mo and W show positive ΔD(Nb1−zBzO2) values since they are competing against Nb4+ reduction.
Interestingly, we find that for MnNb2O6 (MP ID: mp-22100), an older calculation not indexed by MP, is lower in energy. This appears to be due to a process by which MP does not index calculations with old input settings. In this case, the initialized magnetic moments for the Mn atoms were less than 1μB in the newer calculations compared to 5μB in the old ones, despite the expected high-spin state of Mn2+ ions. Notably, when the older, lower energy is used, the minimum aligns more closely with that forecasted by our FEFOS method, indicating that this tool is also useful in detecting outliers within materials databases, and demonstrating its predictive power. We also note that this calculated material matches an ICSD entry despite MP giving it an energy above hull of 0.25 eV per atom. In fact, if we take the older calculation as the ground truth, the formation energy is lowered by 0.29 eV per atom and it becomes below the hull.
In previous work, we showed that the OER descriptor for hypothetical molecular complexes is OS-dependent and that this quantity was not correlated between oxidation states.21,23 Each of these factors are important in thinking about mixing elements at a given M:
O ratio with OER catalysis in mind. As shown in Fig. 5a, an ideal scenario for OER would be one whereby A and B elements in AOx and BOx are non-ideal catalysts, but their reduced and oxidised counterparts, respectively, are. This is illustrated in Fig. 5a where A is reduced and B is oxidised by some unspecified degree, and in Fig. 5b we depict a scenario where this leads their OER activity descriptors to get closer to the ideal value. This picture could explain why doping Ir and Ru oxides with elements which have an affinity for OS 3+ or lower improves OER activity,22,24–30 since the OS of Ir and Ru atoms increases. Indeed, computational results have shown that IrO2 and RuO2 have lower OER activity descriptors than the ideal value,31,32 and that when these metals are in a higher OS it approaches the ideal value. Precisely this explanation from first principles calculation was used to justify the activity of low concentration Ir sites hosted in a WO3 oxide.22 Conversely, Sargent et al. have observed that the activity of NiFe and CoFe oxyhydroxides can be improved by doping them with high-valent elements Mo5+, Ta5+, Nb5+, Re6+ or W6+,13 and this enhancement in catalytic rate was justified on the basis that the Fe atoms were reduced relative to the undoped catalysts.
The spate of data in computational materials databases has led to the creation of machine learning-based techniques to forecast formation energies, which have seen some success. However, these methods suffer severe limitations since they require huge amounts of data to become useful and their predictions may be biased towards previously discovered classes of materials. Furthermore, understanding the outputs of machine learning models is challenging. In contrast, the FEFOS model only requires information about unary oxides, thereby allowing for the automatic screening of any potential combination of elements for a given oxygen content. Another advantage of FEFOS is that the source of stabilisation/destabilisation can be easily interpreted with respect to an element's propensity to be oxidised or reduced.
On a related note, the prediction of the stability of perovskites (with general formula ABX3) has relied on the use of the Goldschmidt tolerance factor,20t, defined as:
![]() | (7) |
![]() | (8) |
Our method is currently restricted to oxides which contain only octahedral sites. Future work will study whether oxides with mixed coordination numbers could be modelled with FEFOS using coordination-specific parabolic equations and consider the relative affinities for oxidation states in more than one type of coordination (e.g. octahedral and tetrahedral sites in a spinel oxide). The applicability of the FEFOS method to such systems relies on the notion that phase diagram surfaces can be drawn as quadratic equations, but the shape of these surfaces could conceivably be more appropriate as a piecewise linear function, or step functions in the case of p-block elements like Sn, which are only known to exist as Sn2+ or Sn4+. Finally, at low concentrations of dopant, the OS of the host material is only slightly shifted, meaning that the resolution of the phase diagram close to the minimum of the convex hull determines the accuracy of predictions.
The pymatgen structure data, coordination data and formation energy data are contained within the binary files unary_data.p, binary_data.p and binary_pairing_data.p within the data_gather folder. To fit the quadratic equations, we set the y-axis as the shift with respect to the reference oxide formation energy for the element unary oxide at each of the relevant oxygen stoichiometries. This reference formation energy was determined to be the oxide with the lowest energy for that element at that stoichiometry. If no such oxide existed in the database, we do not create quadratic equations for that element.
The reference formation energies – where they existed in MP – for each of the 68 elements can be found at: https://github.com/michaelcraiger/fefos/tree/main/supplementary_data in the file named reference_oxides.csv.
The fitted quadratic equations for fox,MO, fred,MO, fox,MO1.5, fred,MO1.5, fox,MO2, and fred,MO2 in eqn (1) for each of the considered elements can be found under the filename parabola_coeffs.csv, which is a spreadsheet containing sheets which each coefficient, as well as separated.csv files with the same data in MO_parabs.csv, M2O3_parabs.csv and MO2_parabs.csv.
All the resultant fitted quadratic equations are plotted along with the associated MP data points for every considered element in the ESI.†
The coefficients were then used to solve eqn (6) such that formation energy predictions for Fig. 2 could be made. The results of each prediction both with and without filtering oxide structure is found under the file fefos_predictions.csv.
To collect Roost test errors for Fig. 3, we collated and parsed the relevant information from ref. 4, with the .csv files contained in the roost_test_data folder of the repository.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cy00107e |
This journal is © The Royal Society of Chemistry 2023 |