# Interpretable machine learning to understand the performance of semi local
density functionals for materials thermochemistry

Santosh Adhikari  Department of Chemistry and Biochemistry, University of
South Carolina, Columbia, SC-29208, USA  Christopher J. Bartel  Department of
Chemical Engineering and Materials Science, University of Minnesota,
Minneapolis, MN 55455, USA  Christopher Sutton<EMAIL_ADDRESS>Department of Chemistry and Biochemistry, University of South Carolina,
Columbia, SC-29208, USA

###### Abstract

This study investigates the use of machine learning (ML) to correct the
enthalpy of formation ($\Delta$Hf) from two separate DFT functionals, PBE and
SCAN, to the experimental $\Delta$Hf across 1011 solid-state compounds. The ML
model uses a set of 25 properties that characterize the electronic structure
as calculated using PBE and SCAN. The ML model significantly decreases the
error in PBE-calculated $\Delta$Hf values from an mean absolute error (MAE) of
195 meV/atom to an MAE = 80 meV/atom when compared to the experiment. However,
a similar reduction in the MAE was not observed for SCAN. Rather, the errors
from the ML model (MAE = 76 meV/atom) and SCAN (MAE = 85 meV/atom) were
observed to be comparable. To explain the substantial decrease in the error of
PBE $\Delta$Hf values and less so for SCAN, we employed partial dependence
plots (PDPs) of an interpretable model, specifically generalized additive
models (GAMs). The PDP+GAM approach allowed for an examination of the impact
of all 25 features on the errors associated with the PBE and SCAN $\Delta$Hf
values. For PBE, the PDP+GAM analysis shows compounds with a high ionicity
(I), i.e., $I$>0.22, have errors in $\Delta$Hf that are twice as large as
compounds having I < 0.22 (246 meV/atom compared to 113 meV/atom). Conversely,
no analogous trend is observed for SCAN-calculated $\Delta$Hfs, which explains
why the ML model for PBE can more easily correct the systematic error in
calculated $\Delta$Hfs for PBE but not for SCAN. Subgroup discovery (SGD) is
used to better understand the relationship between the electronic structure
features and the errors in the PBE calculated $\Delta$Hf. Out of these 25
features, SGD identifies the most reliable region or lowest error subgroup
(108 meV/atom) to be comprised of 368 compounds (out of 1011 total) containing
low charge transfer between atoms based on the selector I<=0.52,
$\zeta_{D}$<=0.43, and $P_{D}$<=0.02. Interestingly, although the literature
suggests PBE is reliable for intermetallics but less so for oxides and
halides, our analysis reveals intermetallics pose a challenge for PBE only
when the charge transfer is significant (I>0.22). Meanwhile, oxides and
halides may be described accurately by PBE for systems in which charge
transfer is relatively low (I < 0.22).

††preprint: APS/123-QED

##  I Introduction

Enthalpy of formation ($\Delta$Hf) and enthalpy of decomposition ($\Delta$Hd)
are the two critical factors for understanding the stability of materials [1].
Because $\Delta$Hd is calculated from the convex hull construction using
$\Delta$Hf for the set of relevant compositions, the development of effective
computational approaches for materials design primarily relies on the accurate
calculation of $\Delta$Hf. Density functional theory (DFT) approximations
based on generalized gradient approximation (GGA) are widely applied to
calculate $\Delta$Hf for solid-state materials. Although the accuracy of DFT-
calculated $\Delta$Hf can depend on the choice of functional, PBE [2] is the
most common GGA that has been used to several hundred of thousand of compounds
in open computational databases such as the Materials Project (MP) [3, 4], the
Open Quantum Materials Database (OQMD) [5, 6], and AFLOW [7].

However, the reliability of the calculated $\Delta$Hf using PBE depends on the
specific material class. For instance, PBE-calculated $\Delta$Hf values are
reasonably accurate for intermetallic alloys, such as Al-based and transition
metal-based alloys [8, 9]. On the other hand, PBE is unreliable in predicting
$\Delta$Hf for systems that combine metals and non-metals, such as oxides [10,
11] and nitrides [12].

The primary source of these increased errors is related to the significant
self-interaction error (SIE) inherent in semilocal functionals such as PBE
[13, 14]. The meta GGA functionals such as SCAN [15], have been demonstrated
to enhance the accuracy of $\Delta$Hf predictions by as much as 50% relative
to PBE [9]. However, issues remain in overestimating $\Delta$Hf for compounds
that are considered to be weakly bound (typically indicated by $\Delta$Hf
values between -1.0 eV/atom to -0.5 eV/atom) [9], such as in the case of
intermetallics. A recent advancement of SCAN, known as r2SCAN [16], has been
shown to further improve the MAE of SCAN by approximately 15 meV/atom for a
dataset of over 1000 solids [17]. However, the persistent problem of
overestimation for intermetallics still remains, as pointed out in Ref [18].
Alternatively, methods like random phase approximations (RPA) and hybrid
functionals, such as HSE06 [19], address some of the limitations of GGA
functionals by, for example, incorporating nonlocality. However, these methods
are computationally demanding and may not always lead to higher accuracies in
the predicted $\Delta$Hf values [20, 21].

A cheaper approach to improving PBE-calculated $\Delta$Hf
($\Delta$H${}_{f}^{PBE}$), which is widely used in databases such as MP, OQMD,
and AFLOW, is to apply an on-site Hubbard U to $d$\- or $f$-orbitals to reduce
the effect of SIE [22, 10]. However, determining the appropriate +U value
remains an open question, and different strategies have been proposed, such as
performing self-consistent calculations using linear response [23, 24, 25] and
tuning +U to recover higher accuracies for specific properties [26, 27].
Additionally, since +U corrections are commonly applied only to states with
$d$ or $f$ characters (eg. strongly correlated materials), SIE associated to
states with $s$ and $p$ characters may still remain [28, 29].

It is also possible to correct much of the errors made by GGA (or GGA+U)
functionals by fitting corrections to the elemental reference energies. This
approach was first shown for oxides [10] and generalized to compounds spanning
the periodic table [30, 31, 3, 6]. Fitted corrections are most effective when
errors in $\Delta$Hf are systematic with respect to certain elements (GGAs)
and were shown to yield less of an improvement for SCAN than for PBE [32].

Alternatively, machine learning (ML) has conventionally been applied to
predict the DFT-calculated $\Delta$Hf of solids at a large scale (> 1000
materials) [33, 34, 35, 36, 37, 38, 39, 40, 41, 42]. However, these models
inherit the underlying errors of the DFT functional (typically PBE). A more
recent study [43] instead uses DFT-calculated $\Delta$Hfs reported in the MP
database (computed using PBE with elemental corrections) to predict
experimental $\Delta$Hfs ($\Delta$H${}_{f}^{expt}$) with 30% more accuracy
compared to the MP database using compositional [33] features (i.e., just
considering the chemical formulas of each compound). Although this study
reports a substantial correction to the DFT $\Delta$Hf errors, it does not
explain the evolution or origin of the error. An understanding of $\Delta$Hf
errors due to the selection of a specific DFT functional, necessitates using
properties as features beyond compositional ones. These properties should be
able to characterize the electron density distribution computed by the chosen
functional.

In this work, we use ML to both understand where the DFT-calculated $\Delta$Hf
values are reliable and correct them to align with the corresponding
corresponding $\Delta$H${}_{f}^{expt}$ values. To achieve this, we initially
formulate a set of numerical electronic-structure-based features that are
sensitive to the choice of functional, for example, charge transfer and the
density of states. We subsequently analyze the influence of these features
using an interpretable ML approach to explain the varying degrees of
improvement in $\Delta$H${}_{f}^{DFT}$ for two DFT functionals, PBE and SCAN.

## Dataset and features

This work utilizes the dataset of $\Delta$H${}_{f}^{expt}$s previously
examined by Bartel et al. [32], which consists of 714 binary, 270 ternary, and
28 quaternary compounds. The dataset spans over 62 elements across the
periodic table and a diverse set of chemical families such as oxides,
sulfides, nitrides, phosphides, halides, and intermetallics.

To generate the electronic features, we use LOBSTER (version 4.1.0) [44] and
Density Derived Electrostatic and Chemical Methods (DDEC6, version 3.5) [45,
46, 47], which post-processes the outputs from VASP calculations (see METHODS
section) to generate a set of features (described below in more detail). The
initial dataset by Bartel et al. [32] contained 1012 compounds; however, we
were only able to generate all the features summarized in Table 1 for 1011
compounds using PBE and 984 compounds using SCAN because of issues with post-
processing via LOBSTER and DDEC6.

Table 1: Features used as inputs into our ML model

Features | Label | Definition | Unit
---|---|---|---
LOBSTER |  $\mu$HP, MAXHP,$~{}\sigma$HP |  averaged, maximum, and standard deviation value of the integral of pCOHP | eV
$\mu$OP, MAXOP,$~{}\sigma$OP |  averaged, maximum, and standard deviation value of the integral of pCOOP | eV
$\mu$BI, MAXBI,$~{}\sigma$BI |  averaged, maximum, and standard deviation value of the integral of COBI | eV
s, p |  relative contribution of s-, p-, d-, and f-orbitals |
d, f |  within $\pm$ 3 eV relative to Fermi level | -
Structure | N |  number of atoms per unit volume of the unit cell | at./Å-3
Eg |  band gap computed as the minimum difference of kohn-sham eigenvalues | eV
CN |  average coordination number of atoms in the cell using Brunners Algorithm [48] | -
$\rho$ |  density of the compound | gm/cm3
$\eta$ |  product of $N$ and $CN$ | at./Å-3
DDEC6 and Bader based-features |  $\zeta$B, $\zeta$D |  average charge transfer using Bader and DDEC6 analysis | e-/at.
PB, PD |  dipole moment per unit volume from Bader and DDEC6 analysis | e-/Å2
p |  net dipole moment per unit cell from DDEC6 analysis | e-/Å2
I |  ratio of net charge transfer to the summed bond order [49] | e-/at.
Atomic |  Z |  a compound index based on modified Pettifor values [50] defined for constituent elements | –

### Structure-based features:

Table 1 lists several structure-based features that are generated from the PBE
and SCAN optimized geometries from Ref. [32]. These features includes the
total number of atoms per unit volume ($N$) and the average coordination
number ($CN$), which is defined as the nearest neighbors for each atom within
a cutoff radius of 8 Å. The $CN$ is calculated using the Brunners algorithm
[48] within pymatgen [51]. We have also included the product of CN and N,
which is labeled as packing ($\eta$).

### LOBSTER-based features:

LOBSTER [44] was used to calculate the integrated projected crystal orbital
overlap population (IpCOOP) [52], projected crystal orbital Hamilton
population (IpCOHP) [53] and crystal orbital bond index (ICOBI) [54] for all
atom pairs separated up to 5 Å. IpCOOP, IpCOHP, and, ICOBI quantify the number
of electrons, contribution to the band-structure energy, and bond index
(degree of covalency, ionicity), respectively, associated with the given bond.
As a result, a set of IpCOOP, IpCOHP, and ICOBI values are calculated for all
unique pairs of atoms within the cutoff radius of 5 Å. For example, for the
ternary compound ABC, each of the three quantities will be calculated for all
A-B, B-C, and A-C atom pairs within 5 Å. To ensure that a consistent set of
nine numbers are computed for each compound, the maximum, average, and the
standard deviation were calculated across each unique pair of atoms for each
compound (e.g., maximum, average, and standard deviation across the values
generated for A-B, B-C, and A-C). If multiple IpCOOP, IpCOHP, and ICOBI values
were produced for each unique pair of atoms in a given compound, potentially
due to minor variations in the bond distance within the structure, we
consolidated this information by taking the average of these values for each
pair.

The normalized contribution of s-, p-, d- and f- orbitals to the bands within
an energy window of $\pm$3 eV around the Fermi level for each compound was
calculated using LOBSTER.

### DDEC6 features:

Several features have been previously used to quantify the charge-transfer
character in compounds such as the net dipole moment per unit cell ($p$) [45,
46] and ionicity ($I$) [49] by evaluating the contribution of each atom in the
system to the total dipole moment and the extent of ionic bonding,
respectively. For illustration, $I$ for a quaternary compound AαBβCθDγ is
computed as:

$I=\frac{1}{\alpha+\beta+\theta+\gamma}({\alpha\frac{\delta_{A}}{s_{A}}+\beta\frac{\delta_{B}}{s_{B}}+\theta\frac{\delta_{C}}{s_{C}}+\gamma\frac{\delta_{D}}{s_{D}}})$
(1)

where $\delta_{A}$ and $s_{A}$, respectively denote the net charges assigned
and summed bond orders for the element A in the compound [49] and so on.
Building off of this previous work, the dipole moment per unit volume ($P$)
was calculated according to the equation:

$P=\frac{1}{2}\frac{l}{V}{\sum_{i=1}^{N}|Q_{i}|},$ (2)

where, $N$ is the total number of atoms, Qi is the net charges acquired
(positive) or lost (negative) by the ith atom in the system, and $l$ is the
distance between the center of positive and negative charges. The average
charge transfer ($\zeta$) in the system is also calculated based on the
equation:

$\zeta=\sum_{i=1}^{N}\frac{|Q_{i}|}{2N},$ (3)

In the expression of both $P$ and $\zeta$, the factor of 2 in the denominator
compensates for the double-counting of the charge transfer (charge acquired
vs. charge lost) and for all samples, both $P$ and $\zeta$ are greater than
zero.
The net charges and dipole moment contribution of each atom computed with
DDEC6 [45, 46, 47] was used to compute $p$ [45, 46], $\it{I}$ [49], $\zeta$D
and PD, where the subscript, $D$, indicates computation using the DDEC6
methods. Separate from DDEC6, we also computed the net charges on each atoms
using Bader analysis [55] and calculated $\zeta$B and PB starting from these
charges (hence, the subscript, $B$, for Bader). Both features are incorporated
as the calculation of atomic charges differs when using DDEC6 in comparison to
using Bader atomic charges, which are both widely used. Bader charges are
calculated by dividing the electron density at zero flux surfaces, whereas
DDEC6-based charges are calculated as a functional of the electron density.
More importantly the computed net atomic charges in both cases are derived
from the total electron density, which naturally incorporates the effects of
SIE.

### Atomic features:

In addition to these electronic features, the compositionally averaged
modified Pettifor index [50] ($Z$) was also used to distinguish between
compounds based on the chemical elements. Modified Pettifor’s index is a
unique value assigned to each element in the periodic table that gives the
measure of the extent of its replaceability in the crystal structure.

## Methods

DFT calculations using the PBE [2] and SCAN [15] functionals were performed
using the projector augmented wave (PAW) formalism as implemented in the
Vienna ab initio simulation package (VASP) code version 5.4.4. The structure
files required for the calculations were obtained from Ref [32]. We used the
PAW pseudopotentials as recommended in the VASP manual, a plane-wave cutoff of
520 eV, energy convergence criterion of 10-6 eV, the smearing parameter
($k_{B}T$) of 0.01 eV (following the first order Methfessel-Paxton scheme),
and gamma-centered Monkhorst-Pack k-point grid with 20|bi| discretizations
along each reciprocal lattice vector, bi, for the Brillouin zone sampling in
all the calculations.

The ML task here adopted the so-called “$\Delta$” learning scheme [56] where
the electronic structure properties from DFT calculations (Table 1) are used
as model inputs to predict the difference between $\Delta$H${}_{f}^{expt}$ and
$\Delta$H${}_{f}^{DFT}$, where the $DFT$ functionals are PBE and SCAN. This
difference is then added back to $\Delta$H${}_{f}^{DFT}$ to predict
$\Delta$H${}_{f}^{expt}$ (with $\Delta$-learned predictions denoted
$\Delta$H${}_{f}^{ML}$). Linear ridge regression (Linear), random forest
regressor (RFR), and kernel ridge regression (KRR) using the Laplacian kernel
(KRR+Lap) and the Gaussian rbf kernel (KRR+rbf) were performed with scikit-
learn [57] version 1.1.1. All hyperparameters were tuned via grid search using
$5$-fold cross-validation. The GAM models were trained using the pyGAM package
[58]. A linear response function (LinearGAM) was used in training the
$\Delta$-learning model. We note that to evaluate the variability of each of
the ML models, 51 distinct models were trained using 51 unique random seeds,
leading to different random 80 % / 20 % training/test splits for each model
(see Table 2 for the mean accuracy and standard deviation of the test error).
$\Delta$H${}_{f}^{DFT/ML}$s were utilized to calculate the enthalpy of
decomposition ($\Delta$H${}_{d}^{DFT/ML}$) for stability analysis using the
approach described in Ref [59].

For consistency of the notation, we use ${\delta}\Delta$H${}_{f}^{DFT}$ and
${\delta}\Delta$H${}_{f}^{ML}$ to represent the difference between $\Delta$Hf
for the DFT functional and ML model, respectively, with experiment (i.e.,
${\delta}\Delta$H${}_{f}^{DFT/ML}~{}$=$~{}\Delta$H${}_{f}^{expt}~{}$ \-
$~{}\Delta$H${}_{f}^{DFT/ML}$). To evaluate model performance, the mean error
(ME) and mean absolute error (MAE) of $\Delta$H${}_{f}^{DFT/ML}$ are computed
relative to $\Delta$H${}_{f}^{expt}$. The ME (MAE) is computed by averaging
${\delta}\Delta$H${}_{f}^{DFT}$ (|${\delta}\Delta$H${}_{f}^{DFT}$|). Figures
S1 - S25 display the distribution of ${\delta}\Delta$H${}_{f}^{DFT}$s versus
all 25 features listed in Table 1.

Finally, for our analysis of the errors for the functional with experiment,
the partial dependence plots (PDPs), generated from the gamma response
function (GammaGAM) were used to predict |${\delta}\Delta$H${}_{f}^{DFT}$|.
Both LinearGAM and GammaGAM are referred to as ‘GAM’ in the results section.
To specifically obtain regions where a given set of $\Delta$H${}_{f}^{DFT}$
has a decreased error relative to experiment (i.e., reliable regions), we used
subgroup discovery (SGD) as described in Ref. [60]. The SGD target variable
for the identification of the reliable regions was the absolute errors
|${\delta}\Delta$H${}_{f}^{DFT}$|. All computations were performed using the
SGD implementation in realKD 0.7.2.

## Results and Discussion

###  $\Delta$-learning $\Delta$Hfs

Table 2 summarizes the average $\Delta$Hf errors calculated from several ML
methods.

Table 2:  Mean absolute errors (meV/atom) predicted by different machine
learning (ML) models trained using the features from PBE and SCAN listed in
Table 1 for the 20% left-out set. The MAEs for each ML model correspond to the
median model across 51 runs. The standard deviations for the 51 runs
(train/test splits) are also provided alongside.

| MAE (meV/atom)
---|---
ML models | PBE features | SCAN features
Linear | 102 $\pm$ 6 | 89 $\pm$ 5
GAM | 91 $\pm$ 6 | 86 $\pm$ 7
KRR+rbf | 92 $\pm$ 6 | 82 $\pm$ 5
KRR+Lap | 84 $\pm$ 6 | 81 $\pm$ 5
RFR | 80 $\pm$ 5 | 76 $\pm$ 5

Figure 1: Parity plot of $\Delta$H${}_{f}^{expt}$ and $\Delta$H${}_{f}^{ML}$
and corresponding histograms of ${\delta}\Delta$H${}_{f}^{DFT}$ (black) and
${\delta}\Delta$H${}_{f}^{ML}$ (gray; for train-set and blue for test-set) for
the RFR model based on (a) PBE and (b) SCAN (Table 2). The mean absolute
error, MAE, of ${\delta}\Delta$H${}_{f}^{ML}$ is 80 meV/atom, and 76 meV/atom
when used with PBE, and SCAN, respectively.

The highest performing $\Delta$-learned ML model for PBE was RFR, which has an
MAE of 80 meV/atom for the 20% test set which is substantially smaller than
PBE (MAE = 195 meV/atom) for the same 20% test set. Fig. 1 shows the
distribution of ${\delta}\Delta$H${}_{f}^{ML}$s, which is equally distributed
above and below zero indicating that the ML model has a mean near zero and
hence no systematic bias which is in stark contrast to PBE. The significant
reduction in ${\delta}\Delta$H${}_{f}^{ML}$ is attributed to the compounds
with a large ${\delta}\Delta$H${}_{f}^{PBE}$ (Fig. 1), which trends in the
number of elements in each compounds. More specifically, the set of compounds
for which |${\delta}\Delta$H${}_{f}^{PBE}$| > 200 meV includes 70% of the
ternary compounds (184 out of 270) and all the 28 quaternary compounds. The
gain in accuracies with $\Delta$H${}_{f}^{ML}$ compared with
$\Delta$H${}_{f}^{PBE}$ is factors of 7.5 ( 43 meV/atom vs. 322 meV/atom) for
quaternary and 3.7 (69 meV/atom vs. 252 meV/atom) for ternary both of which
are much higher than a factor of 2 for binary (86 meV/atom vs. 168 meV/atom).

Figure 2: Heatmap of ${\delta}\Delta$H${}_{f}^{ML}$s vs (a)
${\delta}\Delta$H${}_{f}^{PBE}$s and ${\delta}\Delta$H${}_{f}^{PBE+}$s and (b)
${\delta}\Delta$H${}_{f}^{SCAN}$s and ${\delta}\Delta$H${}_{f}^{SCAN+}$s on
the corresponding 20% left-out test set. The color used in the heatmap is
based on the $\Delta$H${}_{f}^{expt}$ values (in eV) as shown in the colorbar.

Elemental corrections to PBE, so-called PBE+, as fit by Bartel et al. [32] for
the same set of experimental compounds have an MAE = 103 meV/atom, which is
approximately 50% lower than PBE (MAE = 195 meV/atom). Moreover, the MAEs in
$\Delta$H${}_{f}^{PBE+}$ are even for all binary (103, meV/atom), ternary (106
meV/atom), and quaternary (82 meV/atom) compounds, which indicates that the
elemental reference energy corrections address most of the systematic errors
present with PBE. However, the errors are still higher with PBE+ than with the
$\Delta$-learned ML models as can be observed from Fig. 2(a) which shows a
parity plot of ${\delta}\Delta$H${}_{f}^{ML}$ and
${\delta}\Delta$H${}_{f}^{PBE/PBE+}$ indicating that the ML model outperforms
both PBE and PBE+ for the prediction of materials’ $\Delta$Hfs. The overall
reduced errors in $\Delta$H${}_{f}^{ML}$ suggest that the feature set is
contributing additional, material-specific knowledge, which aids in the
increases the model’s accuracy.

Moreover, as mentioned in the Introduction, databases such as OQMD, MP, and,
AFLOW often use +U corrections for specific elements to mitigate the effect of
the SIE and where generally ${\delta}\Delta$H${}_{f}^{PBE}$ is larger.
Analysis of the MAEs of PBE, PBE+, and ML (Table 3) shows that, PBE is most
problematic for oxides of few transition metals and actinides where OQMD would
use +U (365 meV/atom), followed by oxides and fluorides of transition metals
(348 meV/atom) where MP would use +U. The errors are much higher for ternary
and quaternary compounds (432 meV/atom (OQMD), and 418 meV/atom (MP)) compared
to binary compounds. Across each category in Table 3, ML models show
consistently lower MAE of $\sim$ 20 meV/atom compared with PBE+.

Table 3: MAE (in meV/atom) of PBE, PBE+, and ML, for the compounds in the test
set (203) where OQMD, MP, and AFLOW databases would (Y) and would not (N) use
GGA+U. The number inside the parenthesis represents the count of compounds.

| (+U) OQMD | (+U) MP | (+U) AFLOW
---|---|---|---
| Y (14) | N (189) | Y (17) | N (186) | Y (98) | N (105)
PBE | 365 | 182 | 348 | 181 | 198 | 191
PBE+ | 213 | 95 | 209 | 94 | 124 | 84
ML | 200 | 71 | 178 | 71 | 107 | 55

Figure 3: Violin plots showing the distribution of absolute values of
${\delta}\Delta$H${}_{f}^{PBE}$, ${\delta}\Delta$H${}_{f}^{PBE+}$,
${\delta}\Delta$H${}_{f}^{SCAN}$, and ${\delta}\Delta$H${}_{f}^{SCAN+}$ (in
meV/atom) across binary (B), ternary (T), and quaternary (Q) compositions in
our dataset. Out of 1012 compounds, our dataset has 714 binary, 270 ternary,
and 28 quaternary compositions.

For SCAN, ${\delta}\Delta$H${}_{f}^{SCAN}$s has a test-set MAE of 76 meV/atom,
which is comparable to the MAE of SCAN and SCAN fitted with elemental
corrections (SCAN+) of 85 and 65 meV/atom. The mean of
${\delta}\Delta$H${}_{f}^{SCAN}$ is equal to zero indicating that it has no
systematic bias. In contrast to what was observed for $\Delta$H${}_{f}^{PBE}$,
$\Delta$H${}_{f}^{SCAN}$s are relatively more accurate for the ternary (MAE =
59 meV/atom) and quaternary compounds (MAE = 45 meV/atom), but less so for the
binary compounds (MAE = 99 meV/atom), see Fig 3. This supports the findings in
the literature that SCAN is not as accurate as PBE for intermetallics (which
are primarily binary compounds in this dataset) [9].

In order to gain an understanding of the ability of ML to substantially
increase in the accuracy of the predictions of $\Delta$H${}_{f}^{PBE}$, while
providing little improvement for $\Delta$H${}_{f}^{SCAN}$, we have trained an
interpretable GAM using the same $\Delta$-learned procedure previously
described for the RFR model. Despite the fact that the $\Delta$-learned ML
model within the GAM framework exhibits a slightly higher error for
$\Delta$Hfs compared to the RFR model (91 vs. 80 for SCAN; 86 meV/atom vs. 76
meV/atom for PBE, respectively, see Table 2), the advantage of GAMs is the
simple, additive structure which facilitates a more straightforward analysis
of the feature set through partial dependence plots (PDP). The
interpretability of GAMs is preferred over RFR, which does not consider the
effect of each feature separately but instead by design incorporate
interactions between features. As a result, the PDPs generated from RFR-based
models may potentially show unreliable trends, and therefore, are avoided in
this work.

Using the GAM model trained to the set of |${\delta}\Delta$H${}_{f}^{PBE}$|
and |${\delta}\Delta$H${}_{f}^{SCAN}$| values, PDPs were generated (see Figs
S26-S50) for all 25 features inputted into our model (see Table 1). Out of
those 25 different input features, PDP analysis shows that the features
$\zeta$B (in e-/atom), $Z$ and I (in e-/atom) have the highest impact on
|${\delta}\Delta$H${}_{f}^{PBE/SCAN}$|. Although $\zeta$B, Z PDPs displayed
qualitatively similar trends for |${\delta}\Delta$H${}_{f}^{PBE}$| and
|${\delta}\Delta$H${}_{f}^{SCAN}$| values, the PDP of I indicates a stark
difference. Thus, in this text, we focus our discussion on the PDP for I (see
Fig LABEL:pdp:_ionicity).

Figure 5: (a) shows the relation of the $I$ vs ${\delta}\Delta$H${}_{f}^{PBE}$
(left; across whole dataset) and ${\delta}\Delta$H${}_{f}^{ML}$ (right; across
the test set). Blue, yellow, and red dots, respectively, indicate the errors
for binary, ternary, and quaternary compounds. Similarly, (b) shows
${\delta}\Delta$H${}_{f}^{PBE}$ as a function of binary, ternary, and
quaternary compositions where PBE is reliable (RR, blue) vs. unreliable (UR,
gray). The sub-group with $\zeta$D<=0.43 e-/atom, $P$D<=0.02 e-/Å2, and
I<=0.52 e-/atom (Table 1) is reliable (RR), the rest is unreliable (UR).

PDP plot of $I$ shows a positive impact for I > 0.22 on
|${\delta}\Delta$H${}_{f}^{PBE}$|, which indicates a higher error in the
calculation of $\Delta$H${}_{f}^{PBE}$ for compounds with more ionic character
(see Fig LABEL:pdp:_ionicity).

There are 636 compounds that have an I > 0.22. Of which, nearly all (301 out
of 314 compounds in the overall dataset) of the compounds that contain either
alkali or alkaline earth elements are included in this subset. These 301
compounds have larger ${\delta}\Delta$H${}_{f}^{PBE}$ values (MAE of 261
meV/atom), which is comparable to the MAE (246 meV/atom) of all 636 compounds.
The other 13 of 314 alkali or alkaline earth elements compounds, which have an
I < 0.22, the MAE = 71 meV/atom for ${\delta}\Delta$H${}_{f}^{PBE}$.

Although the current literature indicates compounds containing diatomic
elements (O2, Cl2, N2, F2, H2) have a very high error
(${\delta}\Delta$H${}_{f}^{PBE}$), for the 41 (out of a total of 528)
compounds within I < 0.22, the MAE is relatively low (116 meV/atom) despite
containing a few transition metal chlorides (NiCl2 (390 meV/atom), TiCl2 (360
meV/atom), FeCl3 (300 meV/atom), FeCl2 (350 meV/atom)), in this subset. In
comparison, the 487/528 compounds that contain diatomic elements and are in
the range I > 0.22 have an MAE=250 meV/atom.

A similar trend is observed for binary compounds with metallic bonding
(intermetallics), such as CaMg2, AlNi, AlTi, where 13/25 with I < 0.22, have
an MAE of 61 meV/atom. In comparison, the remaining 12 intermetallics with an
I > 0.22, which include Ba2Pb, Ca2Sn, FeTi, etc., have significantly higher
MAE (233 meV/atom).

These results indicate that the evolution of $\Delta$H${}_{f}^{PBE}$ errors
and hence the corrections required for such compounds can vary based on the
bonding environments. This could potentially be influenced by the number of
elements (e.g., binary, ternary, or quaternary) present in the compound, which
could lead to more variation in the bonding of some material.

In contrast, for |${\delta}\Delta$H${}_{f}^{SCAN}$|, a highly uncertain trend
with I is observed as indicated by the flat black line in Figure
LABEL:pdp:_ionicity). The explicit values of the I are also provided as
subfigures in the bottom panels confirming that the qualitative trend with I
changes between the PBE and SCAN. This observation is in agreement with the
previous literature that suggests SCAN improves thermochemical calculations
compared with PBE by better treatment of diversely bonded systems [61].

An examination of the ${\delta}\Delta$H${}_{f}^{ML}$ as a function of $I$,
compound class (binary, ternary, quaternary), and the presence of diatomic
elements within the compound, indicates it is exactly these compounds (e.g.,
compounds with diatomic elements with I>0.22) where a considerable decrease in
${\delta}\Delta$H${f}^{ML}$ compared to ${\delta}\Delta$H${f}^{PBE}$ is
observed (see Fig 5 (a)). This observation indicates there exists a specific
domain formed by conditions on various properties where
${\delta}\Delta$H${}_{f}^{PBE}$ is relatively low, thus constituting a
reliable region. To identify such a reliable region, we employ the approach
introduced in the Ref. [60] and identified a set of 368 compounds (out of
1011) with an MAE = 107 meV/atom for $\Delta$H${}_{f}^{PBE}$ which are
selected by this combination of properties: I<=0.52, $\zeta$D<=0.43, and
PD<=0.02 (Fig 5(b)). In comparison, the remaining 643 compounds have an MAE =
247 meV/atom. The constraint I<=0.52 identified by SGD is in agreement with
the analysis from the PDPs ($I$<0.22) discussed above, but expands the range
of compounds in this reliable region.

Finally, because a material competes with all possible phases in its
compositional space (enthalpy of decomposition; $\Delta$Hd) for stability
rather than just the elemental phases ($\Delta$Hf),[32, 59] we assess how the
lower errors in ${\delta}\Delta$H${}_{f}^{ML}$ translates to $\Delta$Hds
(i.e., $\Delta$H${}_{d}^{ML}$) using the leave-one-chemical-space-out scheme
and compared them with the $\Delta$Hds obtained using $\Delta$H${}_{f}^{expt}$
($\Delta$H${}_{d}^{expt}$) and $\Delta$H${}_{f}^{DFT}$
($\Delta$H${}_{d}^{DFT}$).

As analyzed in Ref [32], the majority of compounds in the MP database [4]
compete for stability with compound phases (Type 2) or a mixture of compounds
and elemental phases (Type 3) rather than elemental phases only (Type 1).
Although $\Delta$H${}_{d}^{DFT}$s (Type 2 and Type 3) are typically 1-2 orders
of magnitude smaller compared to $\Delta$H${}_{f}^{DFT}$s, previous literature
suggests that the ML models trained to predict accurate
$\Delta$H${}_{f}^{DFT}$ values perform poorly for $\Delta$H${}_{d}^{DFT}$s.
[59] This is because errors made by ML models for $\Delta$H${}_{f}^{DFT}$ are
not as systematic as the DFT errors with respect to the chemical composition
[59, 62]. In contrast, our results show that (Table S1) the
$\Delta$H${}_{d}^{ML}$ error is comparable to $\Delta$H${}_{d}^{PBE}$ error
for both type 2 and type 3 compounds suggesting that ML corrections to
$\Delta$H${}_{f}^{PBE}$ are mostly canceling out for stability prediction
($\Delta$Hds). We observe similar findings for $\Delta$H${}_{d}^{ML}$s in the
case of SCAN too (Table S2).

## Conclusion

In this work, we used ML to correct the DFT-calculated enthalpy of formation
($\Delta$H${}_{f}^{DFT}$) to predict experimental enthalpy of formation
($\Delta$H${}_{f}^{expt}$) of 1011 solids. The ML model inputs 25 numerical
properties that capture electronic structure information computed with either
the PBE or SCAN density functional. Among various ML methods, random forest
regressor (RFR) performs best and corrects $\Delta$H${}_{f}^{PBE}$ to within
an MAE = 80 meV/atom of $\Delta$H${}_{f}^{expt}$, significantly improving upon
the MAE of 195 meV/atom of $\Delta$H${}_{f}^{PBE}$s. In contrast, the ML model
correcting $\Delta$H${}_{f}^{SCAN}$ (MAE=76 meV/atom) has comparable accuracy
to SCAN (85 meV/atom).

To analyze why a significant correction is observed for PBE while not so for
SCAN, we studied the impact and trends of all 25 features on the absolute
difference between $\Delta$H${}_{f}^{expt}$ and $\Delta$H${}_{f}^{DFT}$
(errors), using the partial dependence plots (PDPs) within generalized
additive models. PDPs showed that for a number of features, there is a clear
trend between feature values and the $\Delta$H${}_{f}^{PBE}$ errors indicating
that the errors made by PBE are more systematic and straightforward to
correct. More specifically, the interpretable PDP analysis identified the
ionicity (I), which quantifies the ratio of average charge transferred and
summed bond orders, as a strong indicator for high $\Delta$H${}_{f}^{PBE}$
errors. Based on the PDP analysis, we explicitly showed that the MAE for
compounds with I < 0.22 (113 meV/atom) is a factor 2 smaller compared to the
ones with I > 0.22 (246 meV/atom) indicating different evolution of errors
across different regions. This trend from the PDPs is further supported by the
application of subgroup discovery (SGD) to the PBE errors using the same set
of 25 features to identify the reliable regions (i.e., lowest errors). Out of
these 25 features, SGD selects a subgroup with the lowest error (108
meV/atom), containing 368 compounds out of 1011 total, and is predominantly
comprised of compounds containing low charge transfer between atoms based on
the selector: I<=0.52, $\zeta_{D}$<=0.43, and $P_{D}$<=0.02, containing about
one-third of compounds (368/1011) with $\Delta$H${}_{f}^{PBE}$ errors less
than 108 meV/atom. Here $\zeta$ and $P$ quantify the average charge transfer
and dipole moment per volume from DDEC6 analysis. Interestingly, although the
literature suggests PBE is reliable for intermetallics and less so for oxides
and halides, our analysis reveals intermetallics pose a challenge for PBE only
when the charge transfer is significant (I>0.22) and oxides and halides may be
more amenable to PBE if the charge transfer is relatively low (I<0.22), for
example, compounds such as MoCl5 (98 meV/atom), PtCl3 (2 meV/atom), Ag2O (29
meV/atom), PdO (116 meV/atom), etc. Overall, our work leads to a better
understanding of the strengths and limitations of the semilocal density
functional predictions of materials thermochemistry and can potentially guide
the future development of more accurate methods.

###### Acknowledgements.

We acknowledge support by departmental start-up funds at the University of
South Carolina and the NSF-funded MadeinSC (Award number OIA-1655740).

## References

  * Bartel [2022] C. J. Bartel, Review of computational approaches to predict the thermodynamic stability of inorganic solids, Journal of Materials Science 57, 10475 (2022).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Physical Review Letters 77, 3865 (1996).
  * Jain  _et al._ [2011] A. Jain, G. Hautier, C. J. Moore, S. P. Ong, C. C. Fischer, T. Mueller, K. A. Persson, and G. Ceder, A high-throughput infrastructure for density functional theory calculations, Computational Materials Science 50, 2295 (2011).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, _et al._ , Commentary: The Materials Project: A materials genome approach to accelerating materials innovation, APL materials 1, 011002 (2013).
  * Saal  _et al._ [2013] J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, and C. Wolverton, Materials design and discovery with high-throughput density functional theory: the open quantum materials database (OQMD), Jom 65, 1501 (2013).
  * Kirklin  _et al._ [2015] S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl, and C. Wolverton, The Open Quantum Materials Database (OQMD): assessing the accuracy of DFT formation energies, npj Computational Materials 1, 1 (2015).
  * Curtarolo  _et al._ [2012] S. Curtarolo, W. Setyawan, G. L. Hart, M. Jahnatek, R. V. Chepulskii, R. H. Taylor, S. Wang, J. Xue, K. Yang, O. Levy, _et al._ , AFLOW: An automatic framework for high-throughput materials discovery, Computational Materials Science 58, 218 (2012).
  * Wolverton and Ozoliņš [2006] C. Wolverton and V. Ozoliņš, First-principles aluminum database: Energetics of binary Al alloys and compounds, Physical Review B 73, 144104 (2006).
  * Isaacs and Wolverton [2018] E. B. Isaacs and C. Wolverton, Performance of the strongly constrained and appropriately normed density functional for solid-state materials, Physical Review Materials 2, 063801 (2018).
  * Wang  _et al._ [2006] L. Wang, T. Maxisch, and G. Ceder, Oxidation energies of transition metal oxides within the GGA+ U framework, Physical Review B 73, 195107 (2006).
  * Hinuma  _et al._ [2017] Y. Hinuma, H. Hayashi, Y. Kumagai, I. Tanaka, and F. Oba, Comparison of approximations in density functional theory calculations: Energetics and structure of binary oxides, Physical Review B 96, 094102 (2017).
  * Fuchs  _et al._ [2002] M. Fuchs, J. L. F. Da Silva, C. Stampfl, J. Neugebauer, and M. Scheffler, Cohesive properties of group-III nitrides: A comparative study of all-electron and pseudopotential calculations using the generalized gradient approximation, Physical Review B 65, 245212 (2002).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Physical Review B 23, 5048 (1981).
  * Mori-Sánchez  _et al._ [2008] P. Mori-Sánchez, A. J. Cohen, and W. Yang, Localization and delocalization errors in density functional theory and implications for band-gap prediction, Physical Review Letters 100, 146401 (2008).
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly Constrained and Appropriately Normed Semilocal Density Functional, Physical Review Letters 115, 036402 (2015).
  * Furness  _et al._ [2020] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, Accurate and numerically efficient r2SCAN meta-generalized gradient approximation, The Journal of Physical Chemistry Letters 11, 8208 (2020).
  * Kingsbury  _et al._ [2022] R. Kingsbury, A. S. Gupta, C. J. Bartel, J. M. Munro, S. Dwaraknath, M. Horton, and K. A. Persson, Performance comparison of r2SCAN and SCAN metaGGA density functionals for solid materials via an automated, high-throughput computational workflow, Physical Review Materials 6, 013801 (2022).
  * Kothakonda  _et al._ [2022] M. Kothakonda, A. D. Kaplan, E. B. Isaacs, C. J. Bartel, J. W. Furness, J. Ning, C. Wolverton, J. P. Perdew, and J. Sun, Testing the r2SCAN Density Functional for the Thermodynamic Stability of Solids with and without a van der Waals Correction, ACS Materials Au (2022).
  * Krukau  _et al._ [2006] A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, Influence of the exchange screening parameter on the performance of screened hybrid functionals, The Journal of Chemical Physics 125, 224106 (2006).
  * Ong  _et al._ [2011] S. P. Ong, V. L. Chevrier, and G. Ceder, Comparison of small polaron migration and phase separation in olivine LiMnPO 4 and LiFePO 4 using hybrid density functional theory, Physical Review B 83, 075112 (2011).
  * Nepal  _et al._ [2020] N. K. Nepal, S. Adhikari, B. Neupane, and A. Ruzsinszky, Formation energy puzzle in intermetallic alloys: Random phase approximation fails to predict accurate formation energies, Physical Review B 102, 205121 (2020).
  * Dudarev  _et al._ [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. Humphreys, and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+ U study, Physical Review B 57, 1505 (1998).
  * Cococcioni and De Gironcoli [2005] M. Cococcioni and S. De Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the LDA+ U method, Physical Review B 71, 035105 (2005).
  * Kulik  _et al._ [2006] H. J. Kulik, M. Cococcioni, D. A. Scherlis, and N. Marzari, Density functional theory in transition-metal chemistry: A self-consistent Hubbard U approach, Physical Review Letters 97, 103001 (2006).
  * Shishkin and Sato [2016] M. Shishkin and H. Sato, Self-consistent parametrization of DFT+ U framework using linear response approach: Application to evaluation of redox potentials of battery cathodes, Physical Review B 93, 085135 (2016).
  * Loschen  _et al._ [2007] C. Loschen, J. Carrasco, K. M. Neyman, and F. Illas, First-principles LDA+ U and GGA+ U study of cerium oxides: Dependence on the effective U parameter, Physical Review B 75, 035115 (2007).
  * Singh  _et al._ [2015] V. Singh, M. Kosa, K. Majhi, and D. T. Major, Putting DFT to the test: a first-principles study of electronic, magnetic, and optical properties of Co3O4, Journal of Chemical Theory and Computation 11, 64 (2015).
  * Bajaj  _et al._ [2017] A. Bajaj, J. P. Janet, and H. J. Kulik, Communication: Recovering the flat-plane condition in electronic structure theory at semi-local DFT cost, The Journal of Chemical Physics 147, 191101 (2017).
  * Moore  _et al._ [2022] G. C. Moore, M. K. Horton, A. M. Ganose, M. Siron, and K. A. Persson, High-throughput determination of Hubbard U and Hund J values for transition metal oxides via linear response formalism, arXiv preprint arXiv:2201.04213 (2022).
  * Lany [2008] S. Lany, Semiconductor thermochemistry in density functional calculations, Physical Review B 78, 245207 (2008).
  * Stevanović  _et al._ [2012] V. Stevanović, S. Lany, X. Zhang, and A. Zunger, Correcting density functional theory for accurate predictions of compound enthalpies of formation: Fitted elemental-phase reference energies, Physical Review B 85, 115104 (2012).
  * Bartel  _et al._ [2019] C. J. Bartel, A. W. Weimer, S. Lany, C. B. Musgrave, and A. M. Holder, The role of decomposition reactions in assessing first-principles predictions of solid stability, npj Computational Materials 5, 1 (2019).
  * Ward  _et al._ [2016] L. Ward, A. Agrawal, A. Choudhary, and C. Wolverton, A general-purpose machine learning framework for predicting properties of inorganic materials, npj Computational Materials 2, 1 (2016).
  * Deml  _et al._ [2016] A. M. Deml, R. O’Hayre, C. Wolverton, and V. Stevanović, Predicting density functional theory total energies and enthalpies of formation of metal-nonmetal compounds by linear regression, Physical Review B 93, 085142 (2016).
  * Zhang  _et al._ [2020] Z. Zhang, M. Li, K. Flores, and R. Mishra, Machine learning formation enthalpies of intermetallics, Journal of Applied Physics 128, 105103 (2020).
  * Meredig  _et al._ [2014] B. Meredig, A. Agrawal, S. Kirklin, J. E. Saal, J. W. Doak, A. Thompson, K. Zhang, A. Choudhary, and C. Wolverton, Combinatorial screening for new materials in unconstrained composition space with machine learning, Physical Review B 89, 094104 (2014).
  * Jha  _et al._ [2018] D. Jha, L. Ward, A. Paul, W.-k. Liao, A. Choudhary, C. Wolverton, and A. Agrawal, Elemnet: Deep learning the chemistry of materials from only elemental composition, Scientific reports 8, 1 (2018).
  * Goodall and Lee [2020] R. E. Goodall and A. A. Lee, Predicting materials properties without crystal structure: Deep representation learning from stoichiometry, Nature communications 11, 6280 (2020).
  * Huang and Ling [2020] L. Huang and C. Ling, Practicing deep learning in materials science: An evaluation for predicting the formation energies, Journal of Applied Physics 128, 124901 (2020).
  * Schütt  _et al._ [2018] K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko, and K.-R. Müller, Schnet–a deep learning architecture for molecules and materials, The Journal of Chemical Physics 148, 241722 (2018).
  * Xie and Grossman [2018] T. Xie and J. C. Grossman, Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties, Physical Review Letters 120, 145301 (2018).
  * Chen  _et al._ [2019] C. Chen, W. Ye, Y. Zuo, C. Zheng, and S. P. Ong, Graph networks as a universal machine learning framework for molecules and crystals, Chemistry of Materials 31, 3564 (2019).
  * Gong  _et al._ [2022] S. Gong, S. Wang, T. Xie, W. H. Chae, R. Liu, Y. Shao-Horn, and J. C. Grossman, Calibrating DFT formation enthalpy calculations by multifidelity machine learning, JACS Au 2, 1964 (2022).
  * Maintz  _et al._ [2016] S. Maintz, V. L. Deringer, A. L. Tchougréeff, and R. Dronskowski, LOBSTER: A tool to extract chemical bonding from plane-wave based DFT (2016).
  * Manz and Limas [2016] T. A. Manz and N. G. Limas, Introducing DDEC6 atomic population analysis: part 1. Charge partitioning theory and methodology, RSC Advances 6, 47771 (2016).
  * Limas and Manz [2016] N. G. Limas and T. A. Manz, Introducing DDEC6 atomic population analysis: part 2. Computed results for a wide range of periodic and nonperiodic materials, RSC Advances 6, 45727 (2016).
  * Manz [2017] T. A. Manz, Introducing DDEC6 atomic population analysis: part 3. Comprehensive method to compute bond orders, RSC Advances 7, 45552 (2017).
  * Brunner [1977] G. Brunner, A definition of coordination and its relevance in the structure types AlB2 and NiAs, Acta Crystallographica Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography 33, 226 (1977).
  * Sun  _et al._ [2019] W. Sun, C. J. Bartel, E. Arca, S. R. Bauers, B. Matthews, B. Orvañanos, B.-R. Chen, M. F. Toney, L. T. Schelhas, W. Tumas, _et al._ , A map of the inorganic ternary metal nitrides, Nature materials 18, 732 (2019).
  * Glawe  _et al._ [2016] H. Glawe, A. Sanna, E. Gross, and M. A. Marques, The optimal one dimensional periodic table: a modified Pettifor chemical scale from data mining, New Journal of Physics 18, 093011 (2016).
  * Ong  _et al._ [2013] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, and G. Ceder, Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis (2013).
  * Hughbanks and Hoffmann [1983] T. Hughbanks and R. Hoffmann, Chains of trans-edge-sharing molybdenum octahedra: metal-metal bonding in extended systems, Journal of the American Chemical Society 105, 3528 (1983).
  * Dronskowski and Bloechl [1993] R. Dronskowski and P. E. Bloechl, Crystal orbital Hamilton populations (COHP): energy-resolved visualization of chemical bonding in solids based on density-functional calculations, The Journal of Physical Chemistry 97, 8617 (1993).
  * Müller  _et al._ [2021] P. C. Müller, C. Ertural, J. Hempelmann, and R. Dronskowski, Crystal Orbital Bond Index: Covalent Bond Orders in Solids, The Journal of Physical Chemistry C 125, 7959 (2021).
  * Henkelman  _et al._ [2006] G. Henkelman, A. Arnaldsson, and H. Jónsson, A fast and robust algorithm for Bader decomposition of charge density, Computational Materials Science 36, 354 (2006).
  * Ramakrishnan  _et al._ [2015] R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. Von Lilienfeld, Big data meets quantum chemistry approximations: the $\delta$-machine learning approach, Journal of Chemical Theory and Computation 11, 2087 (2015).
  * Pedregosa  _et al._ [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Scikit-learn: Machine Learning in Python, Journal of Machine Learning Research 12, 2825 (2011).
  * [58] https://github.com/dswah/pyGAM.
  * Bartel  _et al._ [2020] C. J. Bartel, A. Trewartha, Q. Wang, A. Dunn, A. Jain, and G. Ceder, A critical examination of compound stability predictions from machine-learned formation energies, npj Computational Materials 6, 97 (2020).
  * Sutton  _et al._ [2020] C. Sutton, M. Boley, L. M. Ghiringhelli, M. Rupp, J. Vreeken, and M. Scheffler, Identifying domains of applicability of machine learning models for materials science, Nature communications 11, 4428 (2020).
  * Sun  _et al._ [2016] J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, _et al._ , Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional, Nature chemistry 8, 831 (2016).
  * Wang  _et al._ [2021] A. Wang, R. Kingsbury, M. McDermott, M. Horton, A. Jain, S. P. Ong, S. Dwaraknath, and K. A. Persson, A framework for quantifying uncertainty in DFT energy corrections, Scientific reports 11, 1 (2021).


$-{{\frac{1}{2}}}e^{mL[C]}\oint_{C}\frac{\delta}{\delta\sigma_{\mu\nu}(s)}\left(ds\frac{e^{-mL[C]}}{L[C]}\frac{\delta\psi[C]}{\delta\sigma^{\mu\nu}(s)}\right)+r\psi[C]=0~{}.$
(D.3)

In (D.3) the action of the area derivative on most of the parts of the
expression is straightforward to understand. There is a subtlety however
associated with the action of the area derivative on $ds$, the integration
measure over the curve itself. To understand this, it is helpful to switch to
a general parameter $\lambda$, in terms of which we have
$ds=d\lambda\sqrt{\dot{X^{2}}}$. This part of the expression can now be
written as

$\frac{\delta}{\delta\sigma_{\mu\nu}(\lambda)}\sqrt{\dot{X^{2}}(\lambda)}=\lim_{\lambda^{\prime}\to\lambda}\frac{\delta}{\delta\sigma_{\mu\nu}(\lambda^{\prime})}\sqrt{\dot{X^{2}}(\lambda)}=-\lim_{\lambda^{\prime}\to\lambda}\frac{2}{\sqrt{\dot{X^{2}}}}\dot{t}_{[\mu}t_{{\alpha}]}\delta(\lambda-\lambda^{\prime})$
(D.4)

where we have used (A.18). The short distance divergence seems to appear from
the fact that we integrate the point where we compute the area derivative over
the whole curve; as we will see, this results in an effective renormalization
of the worldline tension $m$. Returning to the proper length parametrization
$s$, we denote $\delta(s-s^{\prime}=0)=\Lambda$.

The remainder of the expression is more straightforward. Let us now compare
these with the equations of motion previously obtained in Section 3.2. We
consider again the restricted form $\psi[C]=f(A[C])$. Plugging in and using
the expressions for area derivative and length derivatives in (A.17), we find

$-\left(\frac{1}{L[C]}\left(L[C]^{-1}+m+\Lambda\right)\oint
ds\;\dot{t}^{\nu}n_{\nu}\right)f^{\prime}(a)-f^{\prime\prime}(a)+rf(a)=0~{}.$
(D.5)

We now see an issue with this ansatz; the first term in $f^{\prime}(a)$
depends the integral of the tangent and normal vectors, which is not in
general a function of the area of the curve. Thus it is not possible to
satisfy this equation; however as argued in the bulk of the paper it is
possible to find a solution for asymptotically large areas. It would be very
interesting to solve this equation more generally.

##  Appendix E Computation of kinetic term for scalar phase modulation

Here we continue to compute the kinetic term for the field $t(x)$ arising from
evaluating the string field action on the scalar phase modulation ansatz
(4.26). As details have been explained in the bulk of the text, we will here
be somewhat brief. The area derivative is

$\frac{\delta\psi_{t}[C]}{\delta\sigma^{\mu{\alpha}}(s)}=2i\left(\frac{\dot{X}_{[\mu}\partial_{{\alpha}]}t(X(s))}{m\epsilon}-\frac{t(X(s))\ddot{X}_{[\mu}\dot{X}_{{\alpha}]}}{(m\epsilon)^{3}}\right)\psi_{t}[C]\
.$ (E.1)

This is then squared and path-integrated over in the action. The contribution
arising from the square of the second term has already been worked out in the
bulk of the paper; we must thus compute the square of the first term. The
cross term involves an odd number of $X$’s and so will vanish.

The square of the first term is

$\frac{4}{(m\epsilon)^{2}}\left(\dot{X}_{[\mu}\partial_{{\alpha}]}t\dot{X}^{[\mu}\partial^{{\alpha}]}t\right)=(\partial_{{\alpha}}t\partial^{{\alpha}}t)-\frac{2}{(m\epsilon)^{2}}(\dot{X}^{{\alpha}}\partial_{{\alpha}}t)^{2}$
(E.2)

The path integral over the first term involves no new kinematic dependence on
the worldline and so is essentially exactly the same as the computation of the
gauge field kinetic term in (4.25). To determine the second term we need to
work slightly harder. Define the Fourier transform of the function
$T_{{\alpha}\beta}(x)=v^{2}\partial_{{\alpha}}t(x)\partial_{\beta}t(x)$:

$T_{{\alpha}\beta}(x)=\int\frac{d^{d}k}{(2\pi)^{d}}e^{+ik\cdot
x}\tilde{T}_{{\alpha}\beta}(k)$ (E.3)

and then the answer for the contribution to the action is
$\int\frac{d^{d}k}{(2\pi)^{d}}K_{{\alpha}\beta}(k)\tilde{T}^{{\alpha}\beta}(k)$,
where the form factor is:

$K^{{\alpha}\beta}(k)\equiv-\frac{2{{\mathcal{N}}}}{N(m\epsilon)^{2}}(2\pi)^{d}\delta^{(d)}(k)\int\frac{dL}{L^{1+\frac{d}{2}}}e^{-\frac{m^{2}L^{2}}{2}}\oint\frac{ds^{\prime}}{2L}\langle\dot{X}^{{\alpha}}(s^{\prime})\dot{X}^{\beta}(s^{\prime})\rangle$
(E.4)

We now explicitly compute from the definition of the worldline propagator

$\langle\dot{X}^{{\alpha}}(s^{\prime})\dot{X}^{\beta}(s^{\prime})\rangle=-\lim_{s\to
s^{\prime}}\epsilon\left(\delta(s-s^{\prime})-\frac{1}{L}\right)\delta^{{\alpha}\beta}$
(E.5)

As before, we need to regulate the delta function; using the same
representation (4.32) we have $\delta(0)\to\frac{1}{\sqrt{2\pi}\Delta}$. As
$\Delta$ is the cutoff we always have $\Delta\ll L$ and can ignore the $L$
dependence in the integral. Absorbing the $L$ integral into the definition of
${{\mathcal{N}}}$, we find then that the form factor is

$K^{{\alpha}\beta}(k)=(2\pi)^{d}\delta^{(d)}(k)\frac{1}{\sqrt{2\pi}\Delta\epsilon
m^{2}}\delta^{{\alpha}\beta}$ (E.6)

Assembling the pieces, we find the following for the kinetic term:

$S_{\mathrm{kinetic}}[t]=v^{2}\int d^{d}x\left(\partial
t\right)^{2}\left({{\frac{1}{2}}}+\frac{1}{\sqrt{2\pi}\Delta\epsilon
m^{2}}\right)\ .$ (E.7)

## References

  * [1] L. D. Landau and E. M. Lifshitz, Statistical Physics, Part 1, Vol. 5, 3rd Edition.  Elsevier, 1980.
  * [2] F. J. Wegner, “Duality in generalized ising models and phase transitions without local order parameters,” Journal of Mathematical Physics 12 no. 10, (1971) 2259–2272. http://scitation.aip.org/content/aip/journal/jmp/12/10/10.1063/1.1665530.
  * [3] X.-G. Wen, “Topological orders in rigid states,” International Journal of Modern Physics B 4 no. 02, (1990) 239–271.
  * [4] X.-G. Wen, Quantum Field Theory of Many-Body Systems.  Oxford Univ. Press, Oxford, 2004.
  * [5] D. Gaiotto, A. Kapustin, N. Seiberg, and B. Willett, “Generalized Global Symmetries,” JHEP 02 (2015) 172, arXiv:1412.5148 [hep-th].
  * [6] Z. Nussinov and G. Ortiz, “Sufficient symmetry conditions for Topological Quantum Order,” Proc. Nat. Acad. Sci. 106 (2009) 16944–16949, arXiv:cond-mat/0605316.
  * [7] Z. Nussinov and G. Ortiz, “A symmetry principle for topological quantum order,” Annals Phys. 324 (2009) 977–1057, arXiv:cond-mat/0702377.
  * [8] E. Sharpe, “Notes on generalized global symmetries in QFT,” Fortsch. Phys. 63 (2015) 659–682, arXiv:1508.04770 [hep-th].
  * [9] D. M. Hofman and N. Iqbal, “Goldstone modes and photonization for higher form symmetries,” arXiv:1802.09512 [hep-th].
  * [10] A. A. Migdal, “Loop equations and $1/n$ expansion,” Physics Reports 102 no. 4, (1983) 199–290.
  * [11] Y. Makeenko and A. A. Migdal, “Quantum Chromodynamics as Dynamics of Loops,” Sov. J. Nucl. Phys. 32 (1980) 431.
  * [12] A. M. Polyakov, “Gauge Fields as Rings of Glue,” Nucl. Phys. B 164 (1980) 171–188.
  * [13] E. Witten, “Noncommutative Geometry and String Field Theory,” Nucl. Phys. B 268 (1986) 253–294.
  * [14] E. Witten, “On background independent open string field theory,” Phys. Rev. D 46 (1992) 5467–5473, arXiv:hep-th/9208027.
  * [15] S.-J. Rey, “The Higgs Mechanism for Kalb-ramond Gauge Field,” Phys. Rev. D 40 (1989) 3396.
  * [16] M. Franz, “Vortex-boson duality in four space-time dimensions,” EPL 77 no. 4, (2007) 47005, arXiv:cond-mat/0607310.
  * [17] A. J. Beekman, D. Sadri, and J. Zaanen, “Condensing nielsen-olesen strings and the vortex boson duality in 3+1 and higher dimensions,” New Journal of Physics 13 no. 3, (Mar, 2011) 033004. https://doi.org/10.1088/1367-2630/13/3/033004.
  * [18] A. A. Migdal, “Loop equation and area law in turbulence,” Int. J. Mod. Phys. A 9 (1994) 1197–1238, arXiv:hep-th/9310088.
  * [19] S. Mandelstam, “Quantum electrodynamics without potentials,” Annals Phys. 19 (1962) 1–24.
  * [20] A. M. Polyakov, “Gauge Fields and Strings,” Contemp. Concepts Phys. 3 (1987) 1–301.
  * [21] J. Distler, “A Note on the 3-D Ising model as a string theory,” Nucl. Phys. B388 (1992) 648–670, arXiv:hep-th/9205100 [hep-th].
  * [22] A. Kovner and B. Rosenstein, “New look at QED in four-dimensions: The Photon as a Goldstone boson and the topological interpretation of electric charge,” Phys. Rev. D 49 (1994) 5571–5581, arXiv:hep-th/9210154.
  * [23] E. Lake, “Higher-form symmetries and spontaneous symmetry breaking,” arXiv:1802.07747 [hep-th].
  * [24] N. Iqbal and J. McGreevy , to appear .
  * [25] M. J. Strassler, “Field theory without Feynman diagrams: One loop effective actions,” Nucl. Phys. B 385 (1992) 145–184, arXiv:hep-ph/9205205.
  * [26] N. D. Mermin, “The topological theory of defects in ordered media,” Rev. Mod. Phys. 51 (1979) 591–648.
  * [27] J. Polchinski, “Monopoles, duality, and string theory,” Int. J. Mod. Phys. A 19S1 (2004) 145–156, arXiv:hep-th/0304042.
  * [28] J. M. Maldacena, G. W. Moore, and N. Seiberg, “D-brane charges in five-brane backgrounds,” JHEP 10 (2001) 005, arXiv:hep-th/0108152 [hep-th].
  * [29] T. Banks and N. Seiberg, “Symmetries and Strings in Field Theory and Gravity,” Phys. Rev. D83 (2011) 084019, arXiv:1011.5120 [hep-th].
  * [30] M. Kardar, Statistical physics of fields.  Cambridge University Press, 2007.
  * [31] A. M. Somoza, P. Serna, and A. Nahum, “Self-dual criticality in three-dimensional $\mathbb{Z}_{2}$ gauge theory with matter,” arXiv:2012.15845 [cond-mat.stat-mech].
  * [32] S. Grozdanov, D. M. Hofman, and N. Iqbal, “Generalized global symmetries and dissipative magnetohydrodynamics,” Phys. Rev. D95 no. 9, (2017) 096003, arXiv:1610.07392 [hep-th].
  * [33] J. Armas and A. Jain, “Magnetohydrodynamics as superfluidity,” Phys. Rev. Lett. 122 no. 14, (2019) 141603, arXiv:1808.01939 [hep-th].
  * [34] J. Armas and A. Jain, “One-form superfluids & magnetohydrodynamics,” arXiv:1811.04913 [hep-th].
  * [35] N. Iqbal and N. Poovuttikul, “2-group global symmetries, hydrodynamics and holography,” arXiv:2010.00320 [hep-th].
  * [36] A. M. Polyakov, “Compact Gauge Fields and the Infrared Catastrophe,” Phys. Lett. B 59 (1975) 82–84.
  * [37] A. M. Polyakov, “Quark Confinement and Topology of Gauge Groups,” Nucl. Phys. B 120 (1977) 429–458.
  * [38] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, “Topological quantum memory,” J. Math. Phys. 43 (2002) 4452–4505, arXiv:quant-ph/0110143.
  * [39] M. B. Hastings, “Topological order at nonzero temperature,” Physical Review Letters 107 no. 21, (Nov, 2011) . http://dx.doi.org/10.1103/PhysRevLett.107.210501.
  * [40] M. B. Hastings, G. H. Watson, and R. G. Melko, “Self-Correcting Quantum Memories Beyond the Percolation Threshold,” Phys. Rev. Lett. 112 (2014) 070501, arXiv:1309.2680 [cond-mat.str-el].
  * [41] B. Swingle and J. McGreevy, “Mixed s-sourcery: Building many-body states using bubbles of Nothing,” Phys. Rev. B 94 no. 15, (2016) 155125, arXiv:1607.05753 [cond-mat.str-el].
  * [42] M. Creutz, L. Jacobs, and C. Rebbi, “Monte Carlo Computations in Lattice Gauge Theories,” Phys. Rept. 95 (1983) 201–282.
  * [43] Y. Blum, P. Coyle, S. Elitzur, E. Rabinovici, H. Rubinstein, and S. Solomon, “Investigation of the critical behavior of the critical point of the z2 gauge lattice,” Nuclear Physics B 535 no. 3, (1998) 731–738. https://www.sciencedirect.com/science/article/pii/S0550321398006361.
  * [44] T. Filk, M. Marcu, and K. Fredenhagen, “Line of second-order phase transitions in the four-dimensional z2 gauge theory with matter fields,” Physics Letters B 169 no. 4, (1986) 405–412. https://www.sciencedirect.com/science/article/pii/0370269386903813.
  * [45] J. Drouffe, “Large dimensional lattice gauge theories:(i). z2 pure gauge system,” Nuclear Physics B 170 no. 2, (1980) 211–227.
  * [46] O. Akerlund and P. de Forcrand, “U(1) lattice gauge theory with a topological action,” JHEP 06 (2015) 183, arXiv:1505.02666 [hep-lat].
  * [47] H. Kawai, M. Nio, and Y. Okamoto, “On existence of nonrenormalizable field theory: Pure SU(2) lattice gauge theory in five-dimensions,” Prog. Theor. Phys. 88 (1992) 341–350.
  * [48] A. Florio, J. a. M. V. P. L. J. Matos, and J. a. Penedones, “Searching for continuous phase transitions in 5D SU(2) lattice gauge theory,” arXiv:2103.15242 [hep-lat].
  * [49] S. Trebst, P. Werner, M. Troyer, K. Shtengel, and C. Nayak, “Breakdown of a topological phase: Quantum phase transition in a loop gas model with tension,” Phys. Rev. Lett. 98 (2007) 070602, arXiv:cond-mat/0609048.
  * [50] I. S. Tupitsyn, A. Kitaev, N. V. Prokof’ev, and P. C. E. Stamp, “Topological multicritical point in the phase diagram of the toric code model and three-dimensional lattice gauge higgs model,” Physical Review B 82 no. 8, (Aug, 2010) . http://dx.doi.org/10.1103/PhysRevB.82.085114.
  * [51] F. Wu, Y. Deng, and N. Prokof’ev, “Phase diagram of the toric code model in a parallel magnetic field,” Phys. Rev. B 85 (2012) 195104, arXiv:1201.6409 [cond-mat.stat-mech].
  * [52] G. Parisi, “Hausdorff Dimensions and Gauge Theories,” Phys. Lett. B 81 (1979) 357–360.
  * [53] K. Symanzik, Local quantum theory, R. Jost (ed), ch. Euclidean quantum field theory.  Academic Press, New York, London, 1969.
  * [54] D. J. Gross, “The Size of Random Surfaces,” Phys. Lett. B 138 (1984) 185–190.
  * [55] B. Widom, “Surface tension and molecular correlations near the critical point,” The Journal of Chemical Physics 43 no. 11, (1965) 3892–3897, https://doi.org/10.1063/1.1696617. https://doi.org/10.1063/1.1696617.
  * [56] M. Hasenbusch, “Direct Monte Carlo measurement of the surface tension in Ising models,” J. Phys. I(France) 3 no. 3, (1993) 753–765, arXiv:hep-lat/9209016.
  * [57] N. Iqbal and J. McGreevy, “Toward a 3d Ising model with a weakly-coupled string theory dual,” SciPost Phys. 9 no. 2, (2020) 019, arXiv:2003.04349 [hep-th].
  * [58] B. Le Floch, “A slow review of the AGT correspondence,” arXiv:2006.14025 [hep-th].
  * [59] S. S. Razamat and B. Willett, “Down the rabbit hole with theories of class $\mathcal{S}$,” JHEP 10 (2014) 099, arXiv:1403.6107 [hep-th].
  * [60] J. J. Heckman and T. Rudelius, “Top Down Approach to 6D SCFTs,” J. Phys. A 52 no. 9, (2019) 093001, arXiv:1805.06467 [hep-th].
  * [61] L. Zou, Y.-C. He, and C. Wang, “Stiefel liquids: possible non-Lagrangian quantum criticality from intertwined orders,” arXiv:2101.07805 [cond-mat.str-el].
  * [62] T. R. Morris, “Renormalizable extra-dimensional models,” Journal of High Energy Physics 2005 no. 01, (2005) 002.
  * [63] H. Gies, “Renormalizability of gauge theories in extra dimensions,” Phys. Rev. D 68 (2003) 085015, arXiv:hep-th/0305208.