The role of decomposition reactions in assessing first-principles predictions of solid stability

The performance of density functional theory (DFT) approximations for predicting materials thermodynamics is typically assessed by comparing calculated and experimentally determined enthalpies of formation from elemental phases, {\Delta}Hf. However, a compound competes thermodynamically with both other compounds and their constituent elemental forms, and thus, the enthalpies of the decomposition reactions to these competing phases, {\Delta}Hd, determines thermodynamic stability. We evaluated the phase diagrams for 56,791 compounds to classify decomposition reactions into three types: 1. those that produce elemental phases, 2. those that produce compounds, and 3. those that produce both. This analysis shows that the decomposition into elemental forms is rarely the competing reaction that determines compound stability and that approximately two-thirds of decomposition reactions involve no elemental phases. Using experimentally reported formation enthalpies for 1,012 solid compounds, we assess the accuracy of the generalized gradient approximation (GGA) (PBE) and meta-GGA (SCAN) density functionals for predicting compound stability. For 646 decomposition reactions that are not trivially the formation reaction, PBE (MAD = 70 meV/atom) and SCAN (MAD = 59 meV/atom) perform similarly, and commonly employed correction schemes using fitted elemental reference energies make only a negligible improvement (~2 meV/atom). Furthermore, for 231 reactions involving only compounds (Type 2), the agreement between SCAN, PBE, and experiment is within ~35 meV/atom and is thus comparable to the magnitude of experimental uncertainty.


INTRODUCTION
The design and discovery of new materials are being rapidly accelerated by the growing availability of density functional theory (DFT) calculated property data in open materials databases, which allow users to systematically retrieve computed results for experimentally known and yet-to-be-realized solid compounds. [1][2][3][4][5][6] The primary properties of interest are the optimized structure and corresponding total energy, E, with, for example, ~50,000,000 compiled structures and energies available via the NOMAD repository. 7 Given E for a set of structures, one can routinely obtain the reaction energy, Erxn, to convert between structures. E for a compound is typically compared with E for its constituent elements to obtain the formation enthalpy, ΔHf, which provides the thermodynamic driving force at zero temperature and pressure for stability of a given structure with respect to its constituent elements: where E is the calculated total energy of the compound (Aα1Bα2…), αi the stoichiometric coefficient of element i in the compound, and Ei the total energy (chemical potential) of element i. ΔHf computed by

Equation 1
is typically compared to ΔHf obtained experimentally at 298 K with varying degrees of agreement depending on the density functional and compounds (chemistries) under investigation. 2,3,[8][9][10][11][12][13] However, ΔHf is rarely the useful quantity for evaluating the stability of a compound. The reaction energy for a given compound relative to all other compounds within the same composition space is a more relevant metric for accessing stability, where the reaction with the most positive Erxn is the decomposition reaction. 11,14,15 For example, for a given ternary compound, ABC, the relevant space of competing materials includes the elements (A, B, and C), all binary compounds in the A-B, A-C, and B-C spaces, and all ternary compounds in the A-B-C space. The stability of ABC is obtained by comparing the energy of ABC with that of the linear combination of competing compounds with the same average composition as ABC that minimizes the combined energy of the competing compounds, EA-B-C. The decomposition enthalpy, ΔHd, is then obtained by: ΔHd > 0 indicates an endothermic reaction for a given compound ABC forming from the space of competing compounds, A-B-C; the sign notation that ΔHd > 0 indicates instability is chosen to be commensurate with the commonly reported quantity, "energy above the hull", where ΔHd also provides the energy with respect to the convex hull but can be positive (for unstable compounds) or negative (for stable compounds). A ternary example was shown for simplicity, but the decomposition reaction and ΔHd can be obtained for any arbitrary compound comprised of N elements by solving the N-dimensional convex hull problem.
For the high-throughput screening of new materials for a target application, stability against all competing compounds is an essential requirement for determining the viability of a candidate material. 15 In this approach, compounds are typically retained for further evaluation (more rigorous calculations or experiments) if ΔHd < γ, where the threshold γ commonly ranges from ~20 to ~200 meV/atom depending on the priorities of the screening approach and the breadth of materials under evaluation. [16][17][18][19][20][21] The success of high-throughput screening approaches thus depends directly on the accuracy of ΔHd, which is typically obtained using DFT with routinely employed approximations to the exchange-correlation energy.
Nevertheless, despite the intimate link between stability predictions and ΔHd, new approaches (e.g., the development of improved density functionals and/or statistical correction schemes) are primarily benchmarked against experimentally obtained ΔHf. Here, we show that the decomposition reactions that are relevant to stability can be classified into three types, and that the ability of DFT-based approaches to predict ΔHd for each type relative to experiment is the appropriate determinant of the viability of that method for high-throughput predictions of compound stability.

Relevant reactions for determining the stability of compounds
The decomposition reactions that determine ΔHd fall into one of three types: Type 1a given compound is the only known compound in that composition space, the decomposition products are the elements, and thus ΔHd = ΔHf (Fig. 1, left); Type 2a given compound is bracketed (on the phase diagram) by compounds and the decomposition products are exclusively these compounds (Fig. 1, center); and Type 3a given compound is not the only known compound in the composition space, is not bracketed by compounds and the decomposition products are a combination of compounds and elements ( Fig. 1, right).
For a given compound, one of these three types of decomposition reactions will be the relevant reaction for evaluating that material's stability. Notably, these decomposition reactions apply to both compounds that are stable (vertices on the convex hull, ΔHd ≤ 0, Fig. 1, top) and unstable (above the convex hull, ΔHd > 0,

Fig. 1, bottom).
As it pertains to thermodynamic control of synthesis, Type 2 reactions are insensitive to adjustments in elemental chemical potentials that are sometimes modulated by sputtering, partial pressure adjustments, or plasma cracking. Any changes to the elemental energies will affect the decomposition products and the compound of interest proportionally, and therefore, while ΔHf for all compounds will change, ΔHd will be fixed. This is in contrast to Type 1 reactions which become more favorable with increases in the chemical potential of either element. The thermodynamics of Type 3 reactions can be modulated by these synthesis approaches if the elemental form of the species whose chemical potential is being adjusted participates in the decomposition reaction, i.e. the compound must be the nearest (within the convex hull construction) stable, or lowest energy metastable, compound to the element whose chemical potential is being adjusted. 22  Fig. 2, right). As the number of unique elements in the compound, N, increases it becomes increasingly probable that other compounds will be present on the phase diagram and the decomposition will therefore be dictated by these compounds.

Functional performance on formation enthalpy predictions
The decomposition reactions determining compound stability that are Type 1 are the least prevalent among Materials Project compounds (~3%) suggesting that ΔHd rarely equals ΔHf, especially for N > 2 (< 1% of these compounds). Despite this, the primary approach currently used to benchmark first-principles thermodynamics methods is to compare experimental and computed ΔHf. We compared experimentally obtained ΔHf from FactSage 23 to computed ΔHf using the generalized gradient approximation (GGA) density functional as formulated by Perdew, Burke, and Ernzerhof (PBE) 24 and using the strongly constrained and appropriately normed (SCAN) 25 meta-GGA density functionals for 1,012 compounds spanning 62 elements (see Fig. S1 for the prevalence of each element in the evaluated compounds).
Importantly, this reduced space of compounds with experimental thermodynamic data decompose into the full range of Type 1 (37%), 2 (22%), and 3 (41%) reactions. However, we first only analyzed ΔHf for all compounds to establish a baseline for subsequent comparison to ΔHd. On this set of 1,012 compounds, the mean absolute difference (MAD) between experimentally determined ΔHf (at 298 K) 23 and calculated ΔHf, nominally at 0 K and without zero-point energy (ZPE), was found to be 196 meV/atom for PBE and 88 meV/atom for SCAN (Fig. 3a). In addition to a reduction in the magnitude of residuals by ~55%, the distribution of residuals is nearly centered about 0 for SCAN in contrast to PBE which consistently understabilizes compounds relative to their constituent elements (particularly diatomic gases), leading to predictions of ΔHf that are too positive by ~200 meV/atom. Unlike PBE, SCAN has been shown to perform well for a range of diversely bonded systems 25 27 and suggests that the disagreement between experiment and theory should not be lower than ~30 meV/atom on average because this is the magnitude of uncertainty in the experimental determination of ΔHf.
A potential source of disagreement between experimentally obtained and DFT-calculated ΔHf is the incongruence in temperature, where experimental measurements are performed at 298 K and DFT calculations of ΔHf are computed at 0 K, typically neglecting the effects of heat capacity from 0 K to 298 K as well as ZPE. These contributions are typically assumed to be small based on the results obtained for a limited set of compounds. 32 This assumption is robustly confirmed here for 647 structures where the vibrational and heat capacity effects on ΔHf are shown to be ~7 meV/atom on average at 298 K (Fig. S4).
Notably, at higher temperatures, the effects of entropy are significant and should be considered for accurate stability predictions at elevated temperature. 33

Optimizing elemental reference energies
Various approaches have been developed to improve the PBE prediction of ΔHf by systematically adjusting the elemental energies, Ei, of some or all elemental phases. 2,3,[8][9][10] In the fitted elemental reference energy scheme, the difference between experimentally measured and calculated ΔHf is minimized by optimally adjusting Ei by a correction term, δμi: To quantify the magnitude of errors that can be resolved by adjustments to the elemental reference energies, we applied Equation 3 to ΔHf computed with PBE and SCAN ( Fig. 3b) with all elements considered in this optimization (these approaches are denoted in this work as PBE+ and SCAN+, respectively). Fitting reference energies for PBE approximately halves the difference between experiment and calculation and centers the residuals (MAD = 100 meV/atom). Because the difference between experiment and SCAN is less systematic, fitting reference energies improves SCAN errors substantially less than it improves PBE, and only reduces the MAD by ~20% (MAD = 68 meV/atom).  showing that SCAN significantly improves the prediction of ΔHf over PBE. MAD is the mean absolute difference; RMSD is the root-mean-square difference; R 2 is the correlation coefficient; N is the number of compounds shown; μ is the mean difference; σ is the standard deviation. A normal distribution constructed from μ and σ is shown as a solid curve. b) For the same compounds, a comparison of PBE and SCAN with experiment using fitted elemental reference energies for the calculation of ΔHf (PBE+ above; SCAN+ below) showed that for Type 1 reactions fitted elemental reference energies significantly improve the prediction of ΔHf, especially predictions by PBE. These results are provided in Table S1 (for elemental energies) and Table S2 (for compound data). The chemical dependence of these results is shown in Fig. S2a and the distribution of ΔHf,exp is provided in Fig. S5a.

Decomposition reaction analysis
While the improved construction of the SCAN meta-GGA density functional and the use of fitted reference energies ameliorates errors associated with the insufficient description of the elements and thus improves the prediction of ΔHf considerably relative to PBE, the effects these approaches have on the prediction of thermodynamic stability -i.e., ΔHdhave not yet been quantified. We used ΔHf obtained from experiment, PBE, and SCAN for the 1,012 compounds analyzed in Fig. 3 to perform the Ndimensional convex hull analysis to determine the decomposition reaction and quantify ΔHd. For 646 compounds that decompose by Type 2 or 3 reactions, the MAD between experimentally measured and DFT-computed ΔHd is substantially lower than for ΔHf -~60% lower for PBE and ~30% lower for SCAN ( Fig. 4). Notably, the decomposition reaction that results from using experiment, PBE, or SCAN is identical in terms of the competing compounds and their amounts for 89% of the 1,012 compounds evaluated.
For 231 Type 2 decomposition reactions where compounds compete only with compounds and fitted reference energies thus have no influence on ΔHd, SCAN and PBE are found to perform comparably with MADs of ~35 meV/atom compared with experiment. This agreement between theory and experiment using either functional approaches the "chemical accuracy" of experimental measurements (~1 kcal/mol = 22 meV/atom for binary compounds) and is similar to the difference in ΔHf between two experimental sources evaluated in this work (30 meV/atom). A previous study of the formation energies of 135 ternary metal oxides from their constituent binary oxides found that PBE with a Hubbard U correction specifically fit for transition metal oxides achieved a MAD of 24 meV/atom with experiment. 11 The formation of compounds with greater than two elements (ternaries, quaternaries, etc.) from their corresponding binaries is sometimes used as an approximation for ΔHd. 34 computing the decomposition reaction energy using total energies or formation enthalpies is equivalenttherefore the results with (Fig. 4b) and without (Fig. 4a) [25][26][27]36 The discrepancies between the structures and polymorph energy orderings predicted by PBE and SCAN with experiment may also contribute to the reported differences between the approaches.  Fig. S2b and the distribution of ΔHd,exp for Type 2 reactions is provided in Fig. S5b. b) For the same compounds, a comparison of PBE and SCAN with experiment using fitted elemental reference energies for the calculation of ΔHd (PBE+ above; SCAN+ below) showing identical results as (a) due to a cancellation of elemental energies for these Type 2 decomposition reactions. c) A comparison of experimentally measured and DFT-calculated ΔHd (PBE above; SCAN below) for 415 compounds that undergo Type 3 decomposition showing similar performance between PBE and SCAN in predicting ΔHd. d) For the same compounds, a comparison of PBE and SCAN with experiment using fitted elemental reference energies for the calculation of ΔHd (PBE+ above; SCAN+ below) showing that including fitted elemental reference energies does not significantly improve the prediction of ΔHd for Type 3 decomposition reactions. Annotations are as described in the Fig. 3 caption. The chemical dependence of the MAD between theory and experiment for Type 3 reactions is shown in Fig. S2c and the distribution of ΔHd,exp for Type 3 reactions is provided in Fig. S5c.

DISCUSSION
For 1,012 compounds, we show that fitting elemental reference energies for both GGA (PBE) and meta-GGA (SCAN) density functionals improves computed formation enthalpies, ΔHf (Fig. 3). However, to accurately predict the stability of materials, it is essential to accurately compute the decomposition enthalpy, ΔHd, which dictates stability with respect to all compounds and elements in a given chemical space. ΔHd is computed by determining the stoichiometric decomposition reaction with the most positive reaction energy.
ΔHf is only relevant for the stability of compounds that undergo Type 1 decompositions, where the compound only competes with elemental phases and consequently, ΔHd = ΔHf. (Fig. 1). Furthermore, Type 1 decompositions occur for only 17% of binaries and almost never (< 1%) for non-binaries, as shown for the ~60,000 N-component compounds evaluated (Fig. 2). For this reason, ΔHf and the agreement between experiment and theory for ΔHf are rarely relevant to the stability of materials. However, for other applications such as the calculation of defect formation energies, ΔHf is the relevant materials property and the adjustment of calculated chemical potentials using the fitted elemental reference energy scheme may still have significant utility, especially when using PBE. The accuracy of ΔHf is also critical when only select compounds in a given chemical space are not well-described by a given functionale.g., when calculating the stability of peroxides with PBE and the correction developed by Wang et. al, 9 where O2 2− groups are overstabilized. 11,37 If a given error is not systematic for all compounds in a given chemical space, errors in ΔHf may propagate to the errors in ΔHd.
The stabilities of compounds that undergo Type 2 decompositions (63% of compounds tabulated in Materials Project) can be determined without any consideration of elemental energies. For these compounds, PBE and SCAN perform similarly and approach the resolution of experimental approaches to determining ΔHf (~30 meV/atom) (Fig. 4a, Fig. S3). Importantly, the performance metrics we provide are evaluated over a wide range of compounds and chemistries. For chemical spaces that are known to be problematic for a given approach (e.g., 3d transition metals for PBE), the error can significantly exceed the average difference reported here. 27,30 While the majority of compounds in the Materials Project compete with Type 2 decomposition reactions, this is not generally known when first evaluating a compound and so high-throughput screening approaches that typically survey a wide range of compounds will likely include analysis of Type 1 and Type 3 decomposition reactions that do require the calculation of elemental energies. Type 1 decompositions, which occur for binary compounds in sparsely explored chemical spaces, will be highly sensitive to the functional and elemental energies and SCAN improves significantly upon PBE for these compounds. Notably, fitting elemental reference energies for PBE still results in larger errors than SCAN and fitting reference energies for SCAN leads to only modest additional improvements. For Type 3 decompositions, which are ~10 more prevalent than Type 1 reactions in Materials Project, SCAN improves upon PBE by ~20% and the use of fitted elemental reference energies has almost no effect (~2 meV/atom on average) on either approach (Figs. 4c-d). Interestingly, considering the ~60,000 compounds in Materials Project (Fig. 2, left), a roughly equal fraction of Type 2 compounds are stable (48%) and unstable, yet only 37% of Type 3 compounds are stable. However, Type 3 compounds are more amenable to non-equilibrium synthesis approaches that allow for increased chemical potentials of the elements and thus potential access to metastable compounds. 22 In summary, we've shown that the decomposition reactions that dictate the stability of solid compounds

METHODS
Experimental values for ΔHf were obtained from the FactSage database 23 for 1,012 compounds as reported at 298 K and 1 atm. For each compound, the NREL Materials Database (NRELMatDB) 3 was queried for structures matching the composition within 50 meV/atom of the ground-state structure as reported in the database. If a given compound had no calculated structures tabulated in NRELMatDB, the procedure was repeated with the Materials Project database 1 . Structures containing potentially magnetic elements were sampled in non-magnetic, two ferromagnetic (high-and low-spin), and up to 16 antiferromagnetic configurations (depending on cell configuration) where the ground-state magnetic configuration was retained for each structure. Sampling was performed using the approach described by NRELMatDB. This process was also repeated for all 62 elements represented in the dataset with the exceptions of H2, N2, O2, F2, and Cl2 which were calculated as diatomic molecules in a 151515 Å box.
After magnetic sampling, 2,238 unique structures were found for the 1,012 compounds and 62 elements.
All structures were optimized with PBE and SCAN using the Vienna Ab Initio Simulation Package (VASP) 38,39 using the projector augmented wave (PAW) method 40,41 , a plane wave energy cutoff of 520 eV, and a Γ-centered Monkhorst-Pack k-point grid with 20|bi| discretizations along each reciprocal lattice vector, bi. The energy cutoff, k-point density, and related convergence settings were sufficient to achieve total energy convergence of < 5 meV/atom for all calculations. Pseudopotentials used for each element are provided in Table S1. For the calculation of phonons to compute thermal effects, the finite displacement method with 222 supercells as implemented in PHONOPY 42 was used with SCAN and an increased plane wave cutoff of 600 eV and further tightened convergence criteria for total energy convergence of < 1 meV/atom. These results are compiled in Table S3.

Data availability
All necessary data for reproducing this work is contained within, including the computed formation enthalpies, decomposition enthalpies, and chemical potentials that are provided in the Supplementary Information. Additional data are available from the corresponding author upon request.