+m +m +m

#1
#2 #3

+m +m +m #1#2#3

# Nonlocal Machine-Learned Exchange Functional for Molecules and Solids

Kyle Bystrom<EMAIL_ADDRESS>Harvard John A. Paulson School of
Engineering and Applied Sciences  Boris Kozinsky<EMAIL_ADDRESS>Harvard
John A. Paulson School of Engineering and Applied Sciences  Robert Bosch LLC
Research and Technology Center, Cambridge, MA, USA

###### Abstract

The design of better exchange-correlation functionals for Density Functional
Theory (DFT) is a central challenge of modern electronic structure theory.
However, current developments are limited by the mathematical form of the
functional, with efficient semilocal functionals being inaccurate for many
technologically important systems and the more accurate hybrid functionals
being too expensive for large solid-state systems due to the use of the exact
exchange operator. In this work, we use machine learning combined with exact
physical constraints to design an exchange functional that is both orbital-
dependent and nonlocal, but which can be evaluated at roughly the cost of
semilocal functionals and is significantly faster than hybrid DFT in plane-
wave codes. By training functionals with several different feature sets, we
elucidate the roles of orbital-dependent and nonlocal features in learning the
exchange energy and determine that both types of features provide vital and
independently important information to the model. Having trained our new
exchange functional with an expressive, nonlocal feature set, we substitute it
into existing hybrid functionals to achieve hybrid-DFT accuracy on
thermochemical benchmark sets and improve the accuracy of band gap predictions
over semilocal DFT. To demonstrate the scalability of our approach as well as
the practical benefits of improved band gap prediction, we compute charged
defect transition levels in silicon using large supercells. Due to its
transferability and computational efficiency for both molecular and extended
systems, our model overcomes the cost-accuracy trade-off between semilocal and
hybrid DFT, and our general approach provides a feasible path toward a
universal exchange-correlation functional with post-hybrid DFT accuracy and
semilocal DFT cost.

††preprint: APS/123-QED

##  I Introduction

One of the central challenges of density functional theory (DFT) research is
to design accurate approximations to the exchange-correlation (XC) functional
[1]. As the one unknown term in Kohn-Sham DFT [2], errors in the XC functional
limit the predictive power of DFT for a wide variety of technological
problems, including catalysis [3], battery materials [4], and semiconductor
physics [5, 6]. Increasing the complexity of the XC functional can improve its
accuracy, but often at drastically increased computational cost, which can be
prohibitive for large condensed matter simulations. In this work, we use
machine learning (ML) to overcome this cost-accuracy trade-off by training a
functional that obeys exact physical constraints and is accurate, transferable
across a broad range of chemistries, and scalable to hundreds of atoms.

The key cost-accuracy trade-off we address is that between semilocal and
hybrid XC functionals. Semilocal functionals, which include the local density
(LDA), generalized gradient (GGA), and meta-generalized gradient (meta-GGA)
approximations [7], are computationally efficient but insufficiently accurate
for many applications. For example, LDA and GGA functionals are incapable of
correctly predicting band gaps due to their lack of derivative discontinuity
[8, 9], and they are of limited accuracy for many molecular systems [10, 1].
Meta-GGAs introduce derivative discontinuities into the XC functional and can
therefore improve band gap predictions, but they typically account for only
20-50% of the GGA band gap error [11]. Meta-GGAs can be designed to further
improve band gap prediction [12, 13], but this typically makes them less
accurate for other chemical and material properties [14, 13]. On the other
hand, hybrid functionals, which mix a fraction of the nonlocal, orbital-
dependent exact exchange energy into an XC functional [15], provide
significantly improved predictions of band gaps [11, 16] and thermochemical
properties [10, 1] over semilocal DFT, but this accuracy comes at a steep
additional computational cost, especially within plane-wave DFT. In spite of
the remarkable recent progress of efficient implementations for solids [17,
18, 19, 20, 21, 22], hybrid DFT is not as efficient or practical as semilocal
DFT for these calculations.

Several approaches to design improved XC functionals have been pursued in the
literature. Many popular functionals, like PBE [23] and SCAN [24], are
designed primarily or entirely based on mathematical constraints obeyed by the
exact functional. Constraint-based functionals are often considered
trustworthy because they can predict molecular properties without being
explicitly tuned to match them, reducing the risk of bias toward a particular
property or set of systems [25, 26]. However, some of the most popular and
accurate XC functionals for main-group chemistry are empirically fit to
chemical data [27, 28]. Recently, significant progress has been made toward
combining constraint and data-based design approaches using machine learning
[29, 30, 31, 32, 33, 34, 35], culminating in the development of the DM21
functional [36], which achieves state-of-the-art accuracy on main-group
chemistry benchmarks and also improves the description of static correlation
by fitting fractional charge constraints. However, these developments have a
few key shortcomings, the most notable of which is their computational cost.
DM21 uses a local range-separated hybrid formalism, which is both
computationally expensive (even more so than global hybrid DFT) and not widely
available for the plane-wave DFT codes typically used for investigating large
extended systems. As a result, DM21 has not been applied to solids. To the
best of our knowledge, no nonlocal machine-learned functional has yet been
implemented for both molecular (Gaussian-type orbital) and extended (plane-
wave) DFT calculations. All ML XC functionals that have been applied to solids
have been limited to the semilocal form [37, 38, 39, 33, 13, 32] (with the
exception of minimally parameterized nonlocal van der Waals terms). This
limitation restricts the application space of ML functionals by preventing
their use in applications where hybrid DFT accuracy is needed. For example,
hybrid DFT has proven useful for studying problems in catalysis within both
Gaussian-type orbital and plane-wave codes [40, 41] and for addressing the
“band-gap problem” in DFT [42]. In addition to this technical limitation,
applications of ML functionals to complex and heterogeneous systems are
limited by their lack of transferability and universality, as they have
previously only been trained on molecular systems.

These limitations raise the question of whether machine learning models can
leverage more efficient features to balance accuracy and computational cost.
One possible approach in this direction is to generate nonlocal features by
taking integrals of the density over some set of kernel functions [43, 30].
This way, an ML model of the XC functional can capture nonlocal information
about the density without evaluating the computationally demanding exact
exchange operator required for hybrid DFT. In line with this strategy, we
recently introduced the Compressed Scale-Invariant Density Representation
(CIDER) formalism for learning the exchange functional [35], in which the
nonlocal features were designed to enforce the uniform scaling rule [44] on
the learned exchange functional. We demonstrated that the exchange functional
can be accurately learned with the CIDER approach. However, the high
computational cost of computing the features, lack of implementation for
plane-wave DFT, and need to compute and fit the exchange energy density
(rather than the total exchange energy) limited the practicality of our model
and restricted its application to small molecules.

In this work, we overcome all of these limitations by revising the form of the
features, implementing efficient algorithms to evaluate the new features for
both molecular and extended systems, and redesigning our Gaussian process
model to fit the total exchange energy. These improvements enable the
construction of the first nonlocal ML exchange functional for solids, which
can be applied to large systems with similar computational cost to semilocal
DFT. We use this exchange functional to implement efficient and accurate
hybrid DFT surrogates that can be used in a variety of applications. Because
our new functional form is both transferable and numerically stable, we
believe that our approach will also be useful for learning the full XC
functional, allowing us to address problems where DFT accuracy is limited by
both exchange and correlation effects.

By training and testing exchange functionals with different feature sets, we
find that to learn the exchange functional well enough to reproduce hybrid DFT
accuracy, it is necessary to include both a semilocal orbital-dependent
feature (i.e. the kinetic energy density) and nonlocal density-dependent
features in the model input. Having developed this “nonlocal meta-GGA,” we
show that it accurately predicts molecular properties and improves solid-state
band gap predictions over standard meta-GGAs, which is a central problem in
DFT. The performance for solid-state band gaps is particularly promising
because only 9 systems in the training set are solids, with the others being
isolated molecules and complexes. We also encounter only a few convergence
problems in our self-consistent field (SCF) calculations of 2462 molecular and
453 periodic systems, indicating that the functional has good numerical
stability. To demonstrate the transferability and efficiency of the model, we
simulate systems with over 500 atoms to predict charged defect transition
levels in silicon.

The paper is organized as follows. Section II describes the features and
models used in this work. Section III presents our results and discusses the
accuracy of the CIDER functionals as well as their computational efficiency
within plane-wave DFT. Section IV contains our conclusion. Appendix A provides
a detailed description of the computational methods used in this study, and
Appendix B provides additional mathematical details for the plane-wave
implementation of CIDER.

##  II Theory

In Kohn-Sham DFT [2], it is often convenient to separate the XC energy into
exchange and correlation parts, i.e.
$E_{\text{xc}}[n]=E_{\text{x}}[n]+E_{\text{c}}[n]$. The exchange energy
$E_{\text{x}}[n]$ accounts for the antisymmetric nature of the fermionic wave
function, while the correlation energy $E_{\text{c}}[n]$ accounts for
additional many-body correlation effects between the electrons. The exchange
energy is the main focus of this work and can be trivially generalized from
non-spin-polarized DFT to spin-polarized DFT as [45]

$E_{\text{x}}[n_{\uparrow},n_{\downarrow}]=\frac{1}{2}E_{\text{x}}[2n_{\uparrow}]+\frac{1}{2}E_{\text{x}}[2n_{\downarrow}],$
(1)

so all formulas in the theory section are written for the non-spin-polarized
case for simplicity.

The exact exchange energy can be expressed in terms of the occupied Kohn-Sham
orbitals $\\{\phi_{i}(\mathbf{r})\\}$ as

$E_{\text{x}}^{\text{exact}}[n]=-\sum_{ij}\int\differential[3]\mathbf{r}\int\differential[3]\mathbf{r}^{\prime}\frac{\phi_{i}^{*}(\mathbf{r})\phi_{j}^{*}(\mathbf{r})\phi_{i}(\mathbf{r}^{\prime})\phi_{j}(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}.$
(2)

Evaluating Eq. 2 is computationally expensive. However, the improved accuracy
of hybrid DFT arises from mixing a fraction of
$E_{\text{x}}^{\text{exact}}[n]$ into an otherwise semilocal XC functional.
These considerations explain the cost-accuracy trade-off between semilocal and
hybrid DFT. Our goal is to fit Eq. 2 as accurately as possible with a
computationally efficient surrogate model.

The remainder of this section sets up our machine learning-based approach to
this problem. Specifically, Sections II.1-II.3 cover the design of the
exchange functional form, nonlocal features, and Gaussian process models, and
Sections II.4-II.6 cover the computationally efficient implementation of the
nonlocal features. The feature implementation requires different algorithms
for Gaussian-type orbital and plane-wave DFT codes due to the use of different
basis sets and integration grids within these codes.

###  II.1 Uniform Scaling and Semilocal Exchange Functionals

The simplest approximate exchange functional is the local density
approximation (LDA), in which the exchange energy is written as an integral
over a function of the density $n(\mathbf{r})$ at each point in space

$E_{\text{x}}[n]=\int\differential[3]\mathbf{r}\,e_{\text{x}}^{\text{LDA}}\left(n(\mathbf{r})\right),$
(3)

where in the Kohn-Sham formalism [2], the density is expressed in terms of the
occupied orbitals as

$n(\mathbf{r})=\sum_{i}\left|\phi_{i}(\mathbf{r})\right|^{2}.$ (4)

Typically, $e_{\text{x}}^{\text{LDA}}\left(n\right)$ is defined as the
exchange energy density of the uniform electron gas [46]

$e_{\text{x}}^{\text{LDA}}(n)=-\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}n^{4/3}.$
(5)

In addition to making Eq. 3 exact for the uniform electron gas, Eq. 5 also
satisfies an important property of the exchange functional known as uniform
scaling [44, 47], which states that for a scaled density distribution
$n_{\lambda}(\mathbf{r})=\lambda^{3}n(\lambda\mathbf{r})$,

$E_{\text{x}}[n_{\lambda}]=\lambda E_{\text{x}}[n],$ (6)

where $\lambda$ is a positive scalar.

More sophisticated functionals that obey Eq. 6 can be designed by writing the
exchange energy density as the product of Eq. 5 and an exchange enhancement
factor $F_{\text{x}}$:

$E_{\text{x}}[n]=\int\differential[3]\mathbf{r}\,e_{\text{x}}^{\text{LDA}}\left(n(\mathbf{r})\right)F_{\text{x}}(\mathbf{x}(\mathbf{r})).$
(7)

In the above equation, $\mathbf{x}(\mathbf{r})$ is a vector of scale-invariant
features $x_{i}$, where scale-invariance is defined as

$x_{i}[n_{\lambda}](\mathbf{r})=x_{i}[n](\lambda\mathbf{r}).$ (8)

For example, a generalized gradient approximation (GGA) exchange functional is
constructed by making $F_{\text{x}}$ a function of the reduced gradient $s$ of
the density

$s=\frac{|\nabla n|}{2(3\pi^{2})^{1/3}n^{4/3}}.$ (9)

Because $s$ satisfies Eq. 8, the GGA form

$E_{\text{x}}[n]=\int\differential[3]\mathbf{r}\,e_{\text{x}}^{\text{LDA}}\left(n(\mathbf{r})\right)F_{\text{x}}(s(\mathbf{r}))$
(10)

satisfies Eq. 6. Likewise, a meta-GGA can be constructed by making the
exchange enhancement factor a function of two ingredients
$F_{\text{x}}(s,\alpha)$, where

$\alpha=\frac{\tau-\tau_{W}}{\tau_{0}}$ (11)

is a scale-invariant quantity that depends on the kinetic energy density

$\tau(\mathbf{r})=\frac{1}{2}\sum_{i}\left|\nabla\phi_{i}(\mathbf{r})\right|^{2},$
(12)

and $\tau_{W}=\frac{|\nabla n|^{2}}{8n}$ and
$\tau_{0}=\frac{3}{10}(3\pi^{2})^{2/3}n^{5/3}$ are the single-orbital and
uniform electron gas kinetic energies, respectively.

###  II.2 Updated Nonlocal Features for CIDER

In our initial construction of the CIDER formalism [35], the exchange
functional is defined as

$\displaystyle E_{\text{x}}[n]=\int\differential[3]\mathbf{r}\,$
$\displaystyle e_{\text{x}}^{\text{LDA}}\left(n(\mathbf{r})\right)$
$\displaystyle\times
F_{\text{x}}\left(s(\mathbf{r}),\alpha(\mathbf{r}),\\{G_{nlm}(\mathbf{r})\\}\right).$
(13)

In Eq. 13, the functions $\\{G_{nlm}(\mathbf{r})\\}$ are _nonlocal_ features
of the density defined by

$\displaystyle G_{nlm}(\mathbf{r})$
$\displaystyle=\int\differential[3]\mathbf{r}^{\prime}\,\,g_{nlm}(\mathbf{r}^{\prime};\mathbf{r})\,n(\mathbf{r}+\mathbf{r}\,^{\prime}),$
(14)

where $g_{nlm}(\mathbf{r}^{\prime};\mathbf{r})$ is the kernel function

$\displaystyle g_{nlm}(\mathbf{r}^{\prime};\mathbf{r})$
$\displaystyle=M_{l}Y_{lm}(\mathbf{\hat{r}}^{\prime})\left(a(\mathbf{r})|\mathbf{r}^{\prime}|^{2}\right)^{n+l/2}\mathrm{e}^{-a(\mathbf{r})|\mathbf{r}^{\prime}|^{2}},$
(15)

with
$M_{l}=B_{0}^{3/2}\sqrt{4\pi^{l-1}}\left(\frac{8\pi}{3}\right)^{\frac{l}{3}}$.
The position-dependent exponent $a(\mathbf{r})$ is defined as a semilocal
quantity

$\displaystyle a[n](\mathbf{r})$
$\displaystyle=\pi\left(\frac{n(\mathbf{r})}{2}\right)^{2/3}\left[B_{0}+C_{0}\left(\frac{\tau(\mathbf{r})}{\tau_{0}(\mathbf{r})}-1\right)\right],$
(16)

where $B_{0}$ and $C_{0}$ are adjustable constants that control how quickly
the width of the Gaussian-type function in Eq. 15 decreases as the density and
kinetic energy density increase. Equation 16 behaves quadratically under
uniform scaling, i.e.

$a[n_{\lambda}](\mathbf{r})=\lambda^{2}a[n](\lambda\mathbf{r}).$ (17)

This property ensures that $G_{nlm}(\mathbf{r})$ is scale-invariant (Eq. 8)
and therefore that Eq. 13 satisfies the uniform scaling rule (Eq. 6).

Using these scale-invariant nonlocal features, we were able to train an ML
exchange functional that accurately reproduced hybrid DFT atomization energies
[35], but there is room for improvement in making the features more physically
intuitive and easier to implement with production-level performance. In this
work, we simplify the features by removing the spherical harmonics and instead
adding a second position-dependent exponent to the kernel function, resulting
in the new features

$\displaystyle G_{i}(\mathbf{r}_{1})=$
$\displaystyle\left(\frac{B_{i}+B_{0}}{2}\right)^{3/2}$
$\displaystyle\times\int\differential[3]\mathbf{r}_{2}\,\Phi\left(a(\mathbf{r}_{2}),b_{i}(\mathbf{r}_{1}),r_{12}\right)n(\mathbf{r}_{2})$
(18) $\displaystyle\Phi\left(a,b,r\right)=$
$\displaystyle\exp\left[-\left(a+b\right)r^{2}\right],$ (19)

where $r_{12}=|\mathbf{r_{1}}-\mathbf{r_{2}}|$. Multiple features
$G_{1},G_{2},...$ are generated by evaluating Eq. 18 with different choices of
exponent $b_{i}(\mathbf{r})$, which in turn are tuned by constants $B_{i}$ and
$C_{i}$ as discussed below. Because the position-dependent exponents increase
as the density increases, the $a$ exponent effectively damps the value of the
integrand when $\mathbf{r}_{2}$ is near the atomic core, which we intend to
limit the interaction between the valence region and core region. The
prefactor $\left(\frac{B_{i}+B_{0}}{2}\right)^{3/2}$ is chosen so that
$G_{i}=2$ for the non-spin-polarized uniform electron gas.

To retain the scale-invariance of the features, any choice of
$a[n](\mathbf{r})$ and $b_{i}[n](\mathbf{r})$ that satisfies Eq. 17 is
sufficient. For example, the exponents can be an orbital-independent function
of the density and its gradient

$\displaystyle a[n](\mathbf{r})$
$\displaystyle=\pi\left(\frac{n}{2}\right)^{2/3}\left[B_{0}+C_{0}\left(\frac{|\nabla
n|^{2}}{8n\tau_{0}}\right)\right]$ (20) $\displaystyle b_{i}[n](\mathbf{r})$
$\displaystyle=\pi\left(\frac{n}{2}\right)^{2/3}\left[B_{i}+C_{i}\left(\frac{|\nabla
n|^{2}}{8n\tau_{0}}\right)\right]$ (21)

or can be made kinetic-energy dependent, as in our previous work:

$\displaystyle a[n](\mathbf{r})$
$\displaystyle=\pi\left(\frac{n}{2}\right)^{2/3}\left[B_{0}+C_{0}\left(\frac{\tau}{\tau_{0}}-1\right)\right]$
(22) $\displaystyle b_{i}[n](\mathbf{r})$
$\displaystyle=\pi\left(\frac{n}{2}\right)^{2/3}\left[B_{i}+C_{i}\left(\frac{\tau}{\tau_{0}}-1\right)\right].$
(23)

Equations 20 and 21 are used for the nonlocal GGA presented in this work (c.f.
Section II.3), while Eqs. 22 and 23 are used for the nonlocal meta-GGA.
Throughout Sections II.4–II.6, the $\left(\frac{B_{i}+B_{0}}{2}\right)^{3/2}$
factor and the feature index $i$ will be dropped for simplicity, as the
inclusion of the scaling factor and the introduction of additional features by
introducing new $B_{i}$ and $C_{i}$ is trivial.

Before continuing, we briefly discuss the construction of the XC potential
when nonlocal features are employed. A typical semilocal functional has a
multiplicative XC potential

$\functionalderivative{E_{\text{xc}}}{n(\mathbf{r})}=\partialderivative{e_{\text{xc}}(\mathbf{r})}{n(\mathbf{r})}.$
(24)

A meta-GGA also has a non-multiplicative potential

$\functionalderivative{E_{\text{xc}}}{\tau(\mathbf{r})}=\partialderivative{e_{\text{xc}}(\mathbf{r})}{\tau(\mathbf{r})}.$
(25)

This above term is non-multiplicative because it depends on the orbitals
through $\tau(\mathbf{r})$, rather than depending exclusively on the density
$n(\mathbf{r})$. For a given basis set
$\left\\{\chi_{\mu}(\mathbf{r})\right\\}$, the contributions of the XC
potential to the effective Kohn-Sham Hamiltonian matrix are

$\displaystyle V_{\mu\nu}^{\text{xc}}=$
$\displaystyle\int\differential[3]\mathbf{r}\,\chi_{\mu}^{*}(\mathbf{r})\chi_{\nu}(\mathbf{r})\functionalderivative{E_{\text{xc}}}{n(\mathbf{r})}$
$\displaystyle+\frac{1}{2}\int\differential[3]\mathbf{r}\,\nabla\chi_{\mu}^{*}(\mathbf{r})\cdot\nabla\chi_{\nu}(\mathbf{r})\functionalderivative{E_{\text{xc}}}{\tau(\mathbf{r})}.$
(26)

When nonlocal features are included in the XC energy density, e.g.

$E_{\text{xc}}=\int\differential[3]\mathbf{r}\,e_{\text{xc}}\left(n(\mathbf{r}),|\nabla
n(\mathbf{r})|,\tau(\mathbf{r}),G(\mathbf{r})\right),$ (27)

both $\functionalderivative{E_{\text{xc}}}{n(\mathbf{r})}$ and
$\functionalderivative{E_{\text{xc}}}{\tau(\mathbf{r})}$ include terms that
depend on other coordinates $\mathbf{r}^{\prime}$:

$\displaystyle\functionalderivative{E_{\text{xc}}}{n(\mathbf{r})}=$
$\displaystyle\partialderivative{e_{\text{xc}}(\mathbf{r})}{n(\mathbf{r})}+\partialderivative{e_{\text{xc}}(\mathbf{r})}{G(\mathbf{r})}\partialderivative{G(\mathbf{r})}{n(\mathbf{r})}$
$\displaystyle+\int\differential[3]\mathbf{r}^{\prime}\,\partialderivative{e_{\text{xc}}(\mathbf{r}^{\prime})}{G(\mathbf{r}^{\prime})}\functionalderivative{G(\mathbf{r}^{\prime})}{n(\mathbf{r})}$
(28) $\displaystyle\functionalderivative{E_{\text{xc}}}{\tau(\mathbf{r})}=$
$\displaystyle\partialderivative{e_{\text{xc}}(\mathbf{r})}{\tau(\mathbf{r})}+\partialderivative{e_{\text{xc}}(\mathbf{r})}{G(\mathbf{r})}\partialderivative{G(\mathbf{r})}{\tau(\mathbf{r})}$
$\displaystyle+\int\differential[3]\mathbf{r}^{\prime}\,\partialderivative{e_{\text{xc}}(\mathbf{r}^{\prime})}{G(\mathbf{r}^{\prime})}\functionalderivative{G(\mathbf{r}^{\prime})}{\tau(\mathbf{r})}$
(29)

In Eq. 28, $\partialderivative{G(\mathbf{r})}{n(\mathbf{r})}$ describes the
local dependence of $G(\mathbf{r})$ on the density, and
$\functionalderivative{G(\mathbf{r}^{\prime})}{n(\mathbf{r})}$ describes the
nonlocal dependence on the density. While the expressions for
$\functionalderivative{E_{\text{xc}}}{n(\mathbf{r})}$ and
$\functionalderivative{E_{\text{xc}}}{\tau(\mathbf{r})}$ are more complicated
than for a meta-GGA, they still contribute to the XC potential matrix
exclusively via Eq. 26. Therefore, any DFT code supporting meta-GGAs can
evaluate this type of functional if a routine is added to evaluate the
nonlocal contributions to Eqs. 28 and 29.

###  II.3 Gaussian Process Models for the Exchange Functional

####  II.3.1 Overview

To study the impact of different features on the accuracy of an ML exchange
functional, we implement four functional types that use different sets of
features.

SL-GGA is a standard, semilocal GGA like PBE [23], except trained via machine
learning. The exchange enhancement factor ($F_{\text{x}}$ in Eq. 7) is a
function of $s$ only.

NL-GGA is a _nonlocal_ GGA. The exchange enhancement factor is a function of
$s$ and also three nonlocal features of the form in Eq. 18, with Eqs. 20 and
21 for the exponents.

SL-MGGA is a standard, semilocal meta-GGA like SCAN [24], except trained via
machine learning. The exchange enhancement factor is a function of $s$ and
$\alpha$.

NL-MGGA is a _nonlocal_ meta-GGA. The exchange enhancement factor is a
function of $s$, $\alpha$, and three nonlocal features of the form in Eq. 18,
with Eqs. 22 and 23 for the exponents.

In short, the SL-GGA and SL-MGGA are semilocal functionals, while the NL-GGA
and NL-MGGA augment the semilocal functional form with nonlocal features of
the density. As described below, Gaussian process models for the exchange
energy are trained using each of these four feature sets. Section II.3.2
explains how a Gaussian process can be fit to total exchange energies of
systems, Section II.3.3 lists the specific mathematical form of each feature
vector that is used as input to the Gaussian processes, and Section II.3.4
lists the covariance kernels for each feature set.

####  II.3.2 Fitting the Total Exchange Energy

In this section, we introduce a method to train a Gaussian process to the
total exchange energy, as opposed to the nonunique exchange energy density
used in our previous work [35]. While we only learn the exchange energy in
this work, the approach in this section also applies to the XC energy, so we
write the formulas for $E_{\text{xc}}$ for the sake of generality. We index
the chemical systems in the training set by $i$ and $j$; while the feature
implementations are different for Gaussian-type orbital and plane-wave DFT,
the role of systems computed within these two frameworks is identical with
regard to model training. The XC energy of a system $i$ is given by a
numerical integral of the energy density

$E_{\text{xc}}^{i}=\sum_{g\in i}w_{g}^{i}e_{\text{xc}}(\mathbf{x}_{g}^{i}),$
(30)

where $g$ indexes grid points, $w_{g}^{i}$ are the quadrature weights, and
$\mathbf{x}_{g}^{i}$ is the feature vector of molecule $i$ at point $g$.

Constructing a Gaussian process for $E_{\text{xc}}$ requires defining the
covariance between the XC energies for two different systems, i.e.
$\text{Cov}(E_{\text{xc}}^{i},E_{\text{xc}}^{j})$ (c.f. Chapter 2.2 of Ref.
[48]). The covariance of sums of random variables $a,b,c$, and $d$ is

$\displaystyle\text{Cov}(a+b,c+d)=$
$\displaystyle\text{Cov}(a,c)+\text{Cov}(b,c)$
$\displaystyle+\text{Cov}(a,d)+\text{Cov}(b,d).$ (31)

Using Eq. 31, $\text{Cov}(E_{\text{xc}}^{i},E_{\text{xc}}^{j})$ can be
expressed in terms of the covariances between the exchange energy densities as

$\displaystyle\text{Cov}(E_{\text{xc}}^{i},E_{\text{xc}}^{j})=$
$\displaystyle\sum_{g\in i}\sum_{h\in j}w_{g}^{i}w_{h}^{j}$
$\displaystyle\times\text{Cov}(e_{\text{xc}}(\mathbf{x}_{g}^{i}),e_{\text{xc}}(\mathbf{x}_{h}^{j}))$
(32) $\displaystyle=$ $\displaystyle\sum_{g\in i}\sum_{h\in
j}w_{g}^{i}w_{h}^{j}k_{e_{\text{xc}}}(\mathbf{x}_{g}^{i},\mathbf{x}_{h}^{j}),$
(33)

where $k_{e_{\text{xc}}}(\mathbf{x},\mathbf{x}^{\prime})$ is the kernel
function for the covariance between the exchange energy densities for feature
vectors $\mathbf{x}$ and $\mathbf{x}^{\prime}$. Unfortunately, this expression
can be quite expensive to compute because there are tens of thousands of grid
points for each small molecule in the training database, and there are several
hundred training molecules. Therefore, even for a database of this size, at
least $10^{12}$ evaluations of the kernel are required. To circumvent this
problem, a set of control points $\mathbf{\tilde{x}}_{a}$ can be sampled from
the training database to use as a basis set for the feature space (as
described in Appendix A.6), effectively creating a sparse Gaussian process
[48]. Equation 33 is then approximated as

$\displaystyle\text{Cov}(E_{\text{xc}}^{i},E_{\text{xc}}^{j})\approx$
$\displaystyle\left(\mathbf{\tilde{k}}_{e_{\text{xc}}}^{i}\right)^{T}\mathbf{\widetilde{K}}^{-1}\mathbf{\tilde{k}}_{e_{\text{xc}}}^{j}$
(34) $\displaystyle\left(\mathbf{\widetilde{K}}\right)_{ab}=$ $\displaystyle
k_{e_{\text{xc}}}(\mathbf{\tilde{x}}_{a},\mathbf{\tilde{x}}_{b})$ (35)
$\displaystyle\left(\mathbf{\tilde{k}}_{e_{\text{xc}}}^{i}\right)_{a}=$
$\displaystyle\sum_{g\in
i}w_{g}^{i}k_{e_{\text{xc}}}(\mathbf{x}_{g}^{i},\mathbf{\tilde{x}}_{a}).$ (36)

In this work, the exchange functional alone is trained, and the kernel
function is expressed in terms of the exchange enhancement factor
$F_{\text{x}}=e_{\text{x}}/e_{\text{x}}^{\text{LDA}}$. This leads to the
adjusted expressions

$\displaystyle\text{Cov}(E_{\text{x}}^{i},E_{\text{x}}^{j})$
$\displaystyle\approx\left(\mathbf{\tilde{k}}^{i}\right)^{T}\mathbf{\widetilde{K}}^{-1}\mathbf{\tilde{k}}^{j}=K_{ij}$
(37) $\displaystyle\left(\mathbf{\widetilde{K}}\right)_{ab}$
$\displaystyle=k_{F_{\text{x}}}(\mathbf{\tilde{x}}_{a},\mathbf{\tilde{x}}_{b})$
(38) $\displaystyle\left(\mathbf{\tilde{k}}^{i}\right)_{a}$
$\displaystyle=\sum_{g\in
i}w_{g}^{i}(e_{\text{x}}^{\text{LDA}})_{g}^{i}k_{F_{\text{x}}}(\mathbf{x}_{g}^{i},\mathbf{\tilde{x}}_{a}),$
(39)

where $k_{F_{\text{x}}}(\mathbf{x},\mathbf{x}^{\prime})$ is the covariance
kernel for the exchange enhancement factor and
$(e_{\text{x}}^{\text{LDA}})_{g}^{i}=e_{\text{x}}^{\text{LDA}}\left(n_{i}(\mathbf{r}_{g})\right)$
is the LDA exchange energy density at coordinate $g$ for molecule $i$ (see Eq.
5). Having obtained these covariance expressions, they can be plugged into the
predictive function for a Gaussian process to obtain the predicted exchange
energy density [48]

$\displaystyle e_{\text{x}}^{\text{CIDER}}(\mathbf{x}_{*})$
$\displaystyle=e_{\text{x}}^{\text{LDA}}(n_{*})\sum_{a}k_{F_{\text{x}}}(\mathbf{x}_{*},\mathbf{\tilde{x}}_{a})\alpha_{a}$
(40) $\displaystyle\bm{\alpha}$
$\displaystyle=\sum_{i}\mathbf{\tilde{k}}^{i}\left\\{\left[\mathbf{K}+\Sigma_{\text{noise}}\right]^{-1}\mathbf{y}\right\\}_{i},$
(41)

where $n_{*}$ and $\mathbf{x}_{*}$ are the density and feature vector at a
test point and $\Sigma_{\text{noise}}$ is the noise matrix for the training
points.

####  II.3.3 Model Types and Feature vectors

This section lists the feature vectors for each model type in terms of $s$
(Eq. 9), $\alpha$ (Eq. 11), and the nonlocal features $G_{i}$ (Eq. 18). All of
the feature vectors below are designed such that for the uniform electron gas,
$\mathbf{x}=\mathbf{0}$. Since the exchange enhancement factor is
$F_{\text{x}}=1$ for the uniform electron gas [46], the uniform electron gas
constraint is enforced with a single noiseless training point
$F_{\text{x}}^{\text{CIDER}}(\mathbf{0})=1$. All functionals in this work
satisfy the uniform electron gas constraint.

For the SL-GGA, the feature vector $\mathbf{x}$ (Eq. 7) is of dimension 1:

$x_{1}=\frac{\gamma s^{2}}{1+\gamma s^{2}},$ (42)

where $\gamma=0.243$. For the SL-MGGA, the feature vector is of dimension 2:

$\displaystyle x_{1}$ $\displaystyle=\frac{\gamma s^{2}}{1+\gamma s^{2}}$ (43)
$\displaystyle x_{2}$ $\displaystyle=\frac{2}{1+\alpha^{2}}-1.$ (44)

The NL-GGA feature vector is of dimension 4:

$\displaystyle x_{1}$ $\displaystyle=\frac{\gamma s^{2}}{1+\gamma s^{2}}$ (45)
$\displaystyle x_{2}$ $\displaystyle=\frac{G_{1}}{2+G_{1}}-\frac{1}{2}$ (46)
$\displaystyle x_{3}$ $\displaystyle=\frac{G_{2}}{2+G_{2}}-\frac{1}{2}$ (47)
$\displaystyle x_{4}$ $\displaystyle=\frac{G_{3}}{2+G_{3}}-\frac{1}{2}.$ (48)

Lastly, the NL-MGGA feature vector is of dimension 5:

$\displaystyle x_{1}$ $\displaystyle=\frac{\gamma s^{2}}{1+\gamma s^{2}}$ (49)
$\displaystyle x_{2}$ $\displaystyle=\frac{2}{1+\alpha^{2}}-1$ (50)
$\displaystyle x_{3}$ $\displaystyle=\frac{G_{1}}{2+G_{1}}-\frac{1}{2}$ (51)
$\displaystyle x_{4}$ $\displaystyle=\frac{G_{2}}{2+G_{2}}-\frac{1}{2}$ (52)
$\displaystyle x_{5}$ $\displaystyle=\frac{G_{3}}{2+G_{3}}-\frac{1}{2}.$ (53)

The three nonlocal features $G_{1},G_{2},G_{3}$ are given by Eq. 18, and the
parameterization of the kernel exponents (i.e. the constants $B_{i}$ and
$C_{i}$ in Eqs. 20-23) is described in detail in Appendix A.8.

####  II.3.4 Kernels

To fully define the Gaussian process model for the exchange energy, we need to
specify the kernel function $k_{F_{\text{x}}}(\mathbf{x},\mathbf{x}^{\prime})$
and the noise matrix $\Sigma_{\text{noise}}$. The noise matrix is constructed
to heuristically approximate the model uncertainty for each data point as
described in Appendix A.7. To reflect the different number of features in each
functional type, a different kernel is used for
$k_{F_{\text{x}}}(\mathbf{x},\mathbf{x}^{\prime})$ for each, as specified
below:

$\displaystyle k_{\text{SL-GGA}}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=\Sigma_{\text{Cov}}k_{1}(x_{1},x_{1}^{\prime})$ (54)
$\displaystyle k_{\text{SL-MGGA}}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=\Sigma_{\text{Cov}}k_{1}(x_{1},x_{1}^{\prime})k_{2}(x_{2},x_{2}^{\prime})$
(55) $\displaystyle k_{\text{NL-GGA}}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=\Sigma_{\text{Cov}}k_{1}(x_{1},x_{1}^{\prime})$
$\displaystyle\times\sum_{i=2}^{4}\sum_{j=i+1}^{4}k_{i}(x_{i},x_{i}^{\prime})k_{j}(x_{j},x_{j}^{\prime})$
(56) $\displaystyle k_{\text{NL-MGGA}}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=\Sigma_{\text{Cov}}k_{1}(x_{1},x_{1}^{\prime})$
$\displaystyle\times\sum_{i=2}^{5}\sum_{j=i+1}^{5}k_{i}(x_{i},x_{i}^{\prime})k_{j}(x_{j},x_{j}^{\prime}),$
(57)

where $\Sigma_{\text{Cov}}$ is a covariance hyperparamter that can be tuned
for each model, and the base kernel $k_{i}(x,x^{\prime})$ for all of the above
is the squared-exponential kernel

$k_{i}(x,x^{\prime})=\exp[-\frac{1}{2}\left(\frac{x-x^{\prime}}{l_{i}}\right)^{2}].$
(58)

Equations 56 and 57 use the additive Gaussian process approach of Duvenaud _et
al._ [49].

Because the cost of evaluating Gaussian processes scales linearly with the
number of training points, we map the models to cubic splines after training
so they can be evaluated efficiently. The semilocal kernels have 1 and 2
features, so they can be mapped to 1 and 2-dimensional cubic splines,
respectively. The nonlocal kernels are a sum of 3-dimensional kernels, so they
can be mapped to a sum of 3-dimensional cubic splines [35, 50, 51, 52]. For a
description of the hyperparameter selection process for these models, see
Appendix A.8.

###  II.4 Auxiliary Expansion of the CIDER Features

The rest of this section focuses on the efficient implementation of nonlocal
features within Gaussian-type orbital and plane-wave DFT, which is vital for
the practicality of CIDER as well as any nonlocal density functional.
Evaluating Eq. 18 via simple numerical integration is computationally
expensive due to quadratic scaling with the number of grid points. The key
insight to designing a faster approach comes from Román-Pérez and Soler [53],
who developed a method to expand van der Waals density functionals as a sum of
convolutions that could be computed efficiently using the Fast Fourier
Transform (FFT). They observed that a function of the form in Eq. 19, which
depends exclusively on the distance between two points and on a semilocal
functional of the density at each point, can be approximately expanded as [53]

$\Phi\left(a(\mathbf{r}_{2}),b(\mathbf{r}_{1}),r_{12}\right)\approx\sum_{\alpha\beta}p^{a}_{\alpha}(a(\mathbf{r}_{2}))p^{b}_{\beta}(b(\mathbf{r}_{1}))\Phi_{\alpha\beta}(r_{12}).$
(59)

In Eq. 59, $\alpha$ and $\beta$ are indices which each correspond to
interpolation control points $q_{\alpha}$ and $q_{\beta}$, respectively, and
$p^{a}_{\alpha}(a(\mathbf{r}))$ and $p^{b}_{\beta}(b(\mathbf{r}))$ are
interpolating functions which take the values
$p^{a}_{\alpha}(q_{\alpha^{\prime}})=\delta_{\alpha\alpha^{\prime}}$ and
$p^{b}_{\beta}(q_{\beta^{\prime}})=\delta_{\beta\beta^{\prime}}$ at the
interpolating points. The term
$\Phi_{\alpha\beta}(r)\equiv\Phi(q_{\alpha},q_{\beta},r)$.

The Gaussian form of Eq. 19 is particularly conducive to the use of an even-
tempered basis for interpolation, in which $q_{\alpha}$ takes the form

$q_{\alpha}=q_{0}\lambda^{\alpha}$ (60)

for some constant $\lambda>1$. Román-Pérez and Soler [53] used cubic splines
for interpolation, and for the solid-state implementation of CIDER we follow
this approach, with $p^{a}_{\alpha}(a(\mathbf{r}))$ being a cubic spline in
the transformed coordinate
$\gamma(\mathbf{r})=\frac{\ln(a(\mathbf{r})/q_{0})}{\ln(\lambda)}$. However,
for the molecular implementation, we instead treat the control points as a
basis set of Gaussians $e^{-q_{\alpha}r_{12}^{2}}$ onto which the function
$e^{-a(\mathbf{r_{1}})r_{12}^{2}}$ is projected:

$\displaystyle p^{a}_{\alpha}(a(\mathbf{r}))$
$\displaystyle\equiv\sum_{\beta}(S_{a}^{-1})_{\alpha\beta}s_{\alpha}(a(\mathbf{r}))$
(61) $\displaystyle(S_{a})_{\alpha\beta}$
$\displaystyle=\int\differential[3]\mathbf{r}\exp\left[-(q_{\alpha}+q_{\beta})|\mathbf{r}|^{2}\right]$
(62) $\displaystyle s_{\alpha}(a(\mathbf{r}))$
$\displaystyle=\int\differential[3]\mathbf{r^{\prime}}\exp\left\\{-\left(a(\mathbf{r})+q_{\alpha}\right)|\mathbf{r}^{\prime}|^{2}\right\\}.$
(63)

As $\lambda\rightarrow 1$, the number of control points (for the interpolation
scheme) and number of basis functions (for the basis expansion scheme) both
become infinite, and the kernel expansion becomes exact for both approaches.
We determined that the cubic interpolation yields better numerical stability
for the plane-wave implementation, while the basis expansion form is
sufficiently stable for the molecular implementation.

The benefit of expanding the kernel in this manner is that it allows the
nonlocal features (Eq. 18) to be expressed as a sum of convolutions

$G(\mathbf{r}_{1})=\sum_{\beta}p^{b}_{\beta}(b(\mathbf{r}_{1}))\sum_{\alpha}\int\differential[3]\mathbf{r}_{2}\Phi_{\alpha\beta}(r_{12})\theta_{\alpha}(\mathbf{r}_{2}),$
(64)

with $\theta_{\alpha}(\mathbf{r})=p^{a}_{\alpha}(a(\mathbf{r}))n(\mathbf{r})$.
Within plane-wave DFT, the convolutions can be performed with quasi-linear
computational complexity by using FFTs [53]. In addition, because CIDER uses a
squared exponential for $\Phi_{\alpha\beta}(r_{12})$ (Eq. 19), the integrals
of Eq. 64 can be evaluated analytically if $\theta_{\alpha}(\mathbf{r})$ is
projected onto a Gaussian-type orbital basis. We use this approach to design
an efficient implementation of nonlocal features for molecular DFT (Section
II.5).

The Román-Pérez and Soler technique effectively solves the problem of
implementing CIDER for an all-electron plane-wave DFT calculation. However,
such calculations are infeasible for all but the smallest atoms because a
prohibitively large plane-wave basis would be required to describe the core
electrons. In addition, pseudopotentials cannot be reliably used with nonlocal
features because $G(\mathbf{r})$ calculated on the pseudo-density is not
sufficiently similar to the exact, all-electron feature to get accurate
results. To solve this problem, in Section II.6 a method is introduced to
implement nonlocal features within the Projector-Augmented Wave (PAW)
formalism, which requires careful attention to the contributions to
$G(\mathbf{r})$ from the core regions around the atoms.

###  II.5 Efficient Implementation of Nonlocal CIDER Features for Molecular
DFT

Computing Eq. 18 using simple numerical integration is prohibitively expensive
because the number of operations is proportional to the square of the number
of grid points. While van der Waals functionals like VV10 [54] are usually
computed this way within Gaussian-type orbital codes, the contribution of such
functionals is sufficiently small and smooth that a very sparse grid can be
used. The nonlocal features require a grid of the same density as the
semilocal XC grid, and the number of operations per pair of grid points is
higher for evaluating Eq. 18 than, for example, the VV10 kernel. It is
therefore necessary to reduce the computational complexity of evaluating the
nonlocal features for molecular systems.

In this section, we present a charge-partitioning and density fitting-based
approach to compute the nonlocal features within Gaussian-type orbital DFT
codes for molecular systems, which we have implemented for use in the PySCF
code [55, 56]. Much of this approach is based on the work of Franchini et al.
[57], who developed a similar method for computing the Coulomb energy. To
start, the density is partitioned into atomic contributions using Becke
partitioning [58],

$\displaystyle n_{A}(\mathbf{r})$
$\displaystyle=W_{A}(\mathbf{r})n(\mathbf{r})$ (65)
$\displaystyle\sum_{A}W_{A}(\mathbf{r})$ $\displaystyle=1,$ (66)

where $A$ is the atom index and $W_{A}(\mathbf{r})$ are the atomic Becke
weights. The first step in evaluating $G(\mathbf{r})$ is to project
$\theta_{\alpha}(\mathbf{r})$ onto spherical harmonics for each atom (with $L$
being the combined $l$ and $m$ indices and $\mathbf{R}_{A}$ being the position
of atom $A$):

$\displaystyle\theta_{\alpha,L}^{A}(r)$
$\displaystyle=\int\differential\Omega_{A}\theta_{\alpha}(\mathbf{r})W_{A}(\mathbf{r})Y_{L}(\Omega_{A})$
(67) $\displaystyle\mathbf{r}$
$\displaystyle=\mathbf{R}_{A}+(r\cos\phi\sin\theta,r\sin\phi\sin\theta,r\cos\theta),$
(68)

where $\Omega_{A}=(\theta,\phi)$ is the solid angle and $Y_{L}(\Omega)$ are
the spherical harmonics. This projection is done numerically using Lebedev
grids [59]. Then each $L$ channel is projected onto an even-tempered Gaussian
basis $\\{\chi_{A,L,\mu}(\mathbf{r})\\}$ as follows:

$\displaystyle\chi_{A,L,\mu}(\mathbf{r})$
$\displaystyle=Y_{L}(\hat{\mathbf{r}})N_{l,\mu}r^{l}\exp(-\mu r^{2})$ (69)
$\displaystyle\tilde{\theta}_{\alpha,L,\mu}^{A}$
$\displaystyle=\int\differential r\,r^{2+l}N_{l,\mu}\exp(-\mu
r^{2})\theta_{\alpha,L}^{A}(r),$ (70)

where $N_{l,\mu}$ is a normalization constant. The exponents $\mu$ are
constructed using the even-tempered Gaussian approach with tunable parameter
$\beta^{\text{aux}}$

$\mu_{n}=\mu_{\text{min}}\left(\beta^{\text{aux}}\right)^{n}\,\,\,\,\,\,\,n=\\{0,1,...,N\\}$
(71)

A transformation to $\theta_{\alpha,L,\mu}^{A}$ is then performed by
multiplying $\tilde{\theta}_{\alpha,L,\mu}^{A}$ by the inverse overlap matrix
of the $\chi$ basis for each $l$ channel:

$\displaystyle\theta_{\alpha,L,\mu}^{A}$
$\displaystyle=\sum_{\nu}\left[\left(\mathbf{S}_{A,L}^{\chi}\right)^{-1}\right]_{\mu\nu}\tilde{\theta}_{\alpha,L,\nu}^{A}$
(72) $\displaystyle\left(\mathbf{S}_{A,L}^{\chi}\right)_{\mu\nu}$
$\displaystyle=\innerproduct{\chi_{A,L,\mu}}{\chi_{A,L,\nu}}.$ (73)

For each atom, a second even-tempered basis $\\{\xi_{A,L,\mu}(\mathbf{r})\\}$
is introduced to serve as a basis for the convolutions of
$\chi_{A,L,\mu}(\mathbf{r})$. The following projections are computed using
analytical Gaussian integrals:

$\displaystyle\tilde{F}_{\beta,L,\nu}^{A}=$
$\displaystyle\sum_{\alpha,\mu}\theta_{\alpha,L,\mu}^{A}$
$\displaystyle\times\int\differential[3]\mathbf{r}_{1}\differential[3]\mathbf{r}_{2}\,\Phi_{\alpha\beta}(r_{12})\chi_{A,L,\mu}(\mathbf{r}_{1})\xi_{A,L,\nu}(\mathbf{r}_{2}).$
(74)

The transformation to $F_{\beta,L,\nu}^{A}$ is then performed by multiplying
$\tilde{F}_{\beta,L,\nu}^{A}$ by the inverse overlap matrix of the $\xi$
basis, similarly to Eq. 72.

The second to last step is to transform $F_{\beta,L,\nu}^{A}$ to
$F_{\beta}(\mathbf{r}_{g})$, where $\mathbf{r}_{g}$ are the DFT integration
grid points. Because evaluating every $\xi$ orbital of every atom at every
$\mathbf{r}_{g}$ would be computationally expensive, $F_{\beta,L,\nu}^{A}$ is
first transformed to a cubic spline $F_{\beta,L,G,p}^{A}$, where $G$ is a
radial coordinate index and $p$ indexes the polynomial coefficients of the
spline in the interval between points $G$ and $G+1$. The transformation from
$F_{\beta,L,G,p}^{A}$ to $F_{\beta}(\mathbf{r}_{g})$ is then performed as

$F_{\beta}(\mathbf{r}_{g})=\sum_{A}\sum_{L}Y_{L}(\hat{\mathbf{r}}_{gA})\sum_{p=0}^{3}(r_{gA}-r_{G})^{p}F^{A}_{\beta,L,G,p},$
(75)

where $\mathbf{r}_{gA}=\mathbf{r}_{g}-\mathbf{R}_{A}$. The spline index $G$ is
a step function of $|\mathbf{r}_{gA}|$. The computational cost of evaluating
the splines depends heavily on the efficiency of memory access patterns on the
processor, so the coordinates are sorted by $G$ before evaluating the
contributions to $F_{\beta}(\mathbf{r}_{g})$ from each atom.

The last and simplest step to compute $G(\mathbf{r})$ is to evaluate the sum

$G(\mathbf{r})=\sum_{\beta}p^{b}_{\beta}(\mathbf{r}_{g})F_{\beta}(\mathbf{r}_{g}).$
(76)

After computing the XC energy from $G(\mathbf{r})$ and the semilocal features,
the contributions to the XC potential from the $\theta_{\alpha}(\mathbf{r})$
terms can then be computed by passing backward through the matrix
multiplications and spline evaluations detailed above.

The gradient of the total energy with respect to nuclear positions is required
for the evaluation of forces. For this implementation of the CIDER functional
with nonlocal features, the energy gradient with respect to the nuclear
position of atom $A$ is given by

$\displaystyle\nabla_{A}E_{\text{xc}}=$
$\displaystyle\nabla_{A}^{\text{sl}}E_{\text{xc}}+\sum_{g}\left(\sum_{\alpha}\functionalderivative{E_{\text{xc}}}{\theta_{\alpha}(\mathbf{r}_{g})}\theta_{\alpha}(\mathbf{r}_{g})\right)\nabla_{A}w_{g}$
$\displaystyle+\sum_{g}\left(\sum_{\beta}\functionalderivative{E_{\text{xc}}}{F_{\beta}(\mathbf{r}_{g})}\nabla_{A}F_{\beta}(\mathbf{r}_{g})\right)w_{g},$
(77)

where $w_{g}$ is the numerical integration quadrature weight corresponding to
the coordinate $\mathbf{r}_{g}$. The analytical energy gradient for a
semilocal functional $\nabla_{A}^{\text{sl}}E_{\text{xc}}$ is given by Johnson
_et al._ [60] Eq. 11.

The computational bottleneck of this algorithm is the spline evaluation in Eq.
75, whose cost scales quadratically with system size. Specifically, the cost
scales as
$\Theta(N_{\text{grid}}\,N_{\text{atom}}\,N_{\beta}\,l_{\text{max}}^{2})$,
where $N_{\text{grid}}$ is the number of integration grid points,
$N_{\text{atom}}$ is the number of atoms, $N_{\beta}$ is the number of kernel
interpolation points in Eq. 59, and $l_{\text{max}}$ is the maximum spherical
harmonic order used in the $\xi$ basis. All other components of the algorithm
are performed on individual atoms and therefore scale linearly with system
size.

###  II.6 Quasi-linear-scaling Implementation of Nonlocal CIDER Features for
Plane-Wave DFT with PAW

This section covers our algorithm for evaluating the nonlocal CIDER features
within the Projector-Augmented Wave (PAW) method [61], a generalization of the
pseudopotential method that can recover the full wave function and electron
density. We have implemented this algorithm in the GPAW code [62, 63]. As
discussed above, Eqs. 59 and 64 enable the evaluation of $G(\mathbf{r})$ in
$O(N\log N)$ complexity if the density $n(\mathbf{r})$ is represented on an
FFT grid. However, in the PAW formalism of Blöchl [61], the all-electron
density $n(\mathbf{r})$ is represented as

$n(\mathbf{r})=\tilde{n}(\mathbf{r})+\sum_{A}(n_{A}^{1}(\mathbf{r})-\tilde{n}_{A}^{1}(\mathbf{r})),$
(78)

where $A$ indexes the atoms, $\tilde{n}(\mathbf{r})$ is the pseudo-density on
the FFT grid, and the terms with the 1 superscript are represented on radial
support grids centered around each atom and extending to some cutoff radius
$r_{c}$. Within the cutoff radius of atom $A$,
$\tilde{n}(\mathbf{r})=\tilde{n}_{A}^{1}(\mathbf{r})$ and
$n_{A}^{1}(\mathbf{r})$ is the all-electron density. Outside the cutoff radius
of atom $A$, $n_{A}^{1}(\mathbf{r})=\tilde{n}_{A}^{1}(\mathbf{r})$.

For a semilocal functional, Eq. 78 yields a simple expression for the XC
energy $E_{\text{xc}}[n]$:

$E_{\text{xc}}[n]=E_{\text{xc}}[\tilde{n}]+\sum_{A}(E_{\text{xc}}[n_{A}^{1}]-E_{\text{xc}}[\tilde{n}_{A}^{1}]).$
(79)

However, the case is not so simple for a nonlocal functional such as CIDER. If
Eq. 18 is evaluated on the pseudo-density $\tilde{n}(\mathbf{r})$, the
resulting $G(\mathbf{r})$ is not equal to the all-electron feature even
outside the cutoff radius because the convolutions used to compute the feature
are nonlocal. For the same reason, evaluating Eq. 18 on the atomic density
$n_{A}^{1}(\mathbf{r})$ does not yield the all-electron feature inside the
cutoff radius. Therefore, an implementation of CIDER (and other nonlocal
functionals with similar forms, such as some kinetic energy and van der Waals
functionals) for PAW must construct a “pseudo-feature” $\tilde{G}(\mathbf{r})$
that is smooth but yields the correct feature outside the cutoff radius, as
well as construct the “all-electron feature” $G(\mathbf{r})$ on the radial
support grids inside the cutoff radius. Because of this difficulty, current
plane-wave implementations of nonlocal functionals like van der Waals
correlation [64] use the pseudo-density rather than the full PAW density. When
this approach is used, one may need to use PAW datasets with more valence
electrons, increasing computational cost [64]. In addition, because CIDER
computes the exchange energy, which is of much larger magnitude than the van
der Waals correlation energy, higher precision is required. A similar problem
exists for kinetic energy density functionals. There is an implementation of
semilocal kinetic energy functionals with PAW [65], and of nonlocal kinetic
energy functionals with ultrasoft pseudopotentials [66], but no nonlocal
kinetic energy functional implementation within PAW. Therefore, by designing
an algorithm for computing nonlocal features for CIDER within the PAW
formalism, we also enable the full integration of van der Waals functionals
and orbital-free DFT [67] into the PAW method.

Starting from Eq. 78, the all-electron value of $G(\mathbf{r})$ _outside the
core region_ is given by

$\displaystyle\tilde{G}(\mathbf{r})$
$\displaystyle=\sum_{\beta}p^{b}_{\beta}\left(\tilde{b}(\mathbf{r})\right)\left[\tilde{F}_{\beta}(\mathbf{r})+\sum_{A}\Delta
F_{A,\beta}(\mathbf{r})\right]$ (80)
$\displaystyle\tilde{F}_{\beta}(\mathbf{r}_{1})$
$\displaystyle=\sum_{\alpha}\int\differential[3]\mathbf{r}_{2}\,\Phi_{\alpha\beta}(r_{12})\tilde{\theta}_{\alpha}(\mathbf{r}_{2})$
(81) $\displaystyle\Delta F_{A,\beta}(\mathbf{r}_{1})$
$\displaystyle=\sum_{\alpha}\int\differential[3]\mathbf{r}_{2}\,\Phi_{\alpha\beta}(r_{12})\left(\theta^{1}_{A,\alpha}(\mathbf{r}_{2})-\tilde{\theta}^{1}_{A,\alpha}(\mathbf{r}_{2})\right).$
(82)

The quantities $\tilde{b}(\mathbf{r})$ and
$\tilde{\theta}_{\alpha}(\mathbf{r})$ are evaluated on the pseudo-density
$\tilde{n}(\mathbf{r})$; $\theta^{1}_{A,\alpha}(\mathbf{r})$ and
$\tilde{\theta}^{1}_{A,\alpha}(\mathbf{r})$ are evaluated on
$n_{A}^{1}(\mathbf{r})$ and $\tilde{n}_{A}^{1}(\mathbf{r})$, respectively.
$\tilde{F}_{\beta}(\mathbf{r})$ can be evaluated simply using FFTs. To compute
$\Delta F_{A,\beta}(\mathbf{r})$, the term
$\Delta\theta^{1}_{A,\alpha}(\mathbf{r})=\theta^{1}_{A,\alpha}(\mathbf{r})-\tilde{\theta}^{1}_{A,\alpha}(\mathbf{r})$
is first projected onto spherical harmonics and Gaussians as described in the
previous section. Because $\Delta\theta^{1}_{A,\alpha}(\mathbf{r})$ is
localized within the augmentation region and mostly spherical, far fewer
spherical harmonics and Gaussians are required per atom than for the Gaussian-
type orbital implementation. This projection yields
$\Delta\theta_{A,\alpha,L,\mu}^{1}$. The Gaussians are then analytically
projected into reciprocal space to yield $\Delta\theta_{A,\alpha,L}^{1}(k)$.
In reciprocal space,

$\Delta
F_{A,L,\beta}(k)=\sum_{\alpha}\Phi_{\alpha\beta}(k)\Delta\theta_{A,\alpha,L}^{1}(k).$
(83)

The above equation applies because the convolution $\Phi_{\alpha\beta}(k)$ has
angular momentum number $l=0$. Note that in this section, denoting a symbol as
function of $k$ implies that it is the Fourier transform of the real-space
function corresponding to that symbol. For example, $\Phi_{\alpha\beta}(k)$ is
the Fourier transform of $\Phi_{\alpha\beta}(r)$.

Inspired by the PAW method itself, we perform the following trick to
approximate $\Delta F_{A,\beta}(\mathbf{r})$ on the FFT grid. We introduce two
sets of functions, $\\{g_{i}^{A}(\mathbf{r})\\}$ and
$\\{h_{j}^{A}(\mathbf{r})\\}$. The $h_{j}^{A}(\mathbf{r})$ functions form a
localized basis (i.e. each function vanishes for $r>r_{c}$) and are meant to
fit the high-frequency components of $\Delta F_{A,\beta}(\mathbf{r})$
introduced by the all-electron density inside the core region. For the current
implementation of these functions, we use differences between Gaussians and
smooth polynomials that match their derivatives at $r_{c}$. The other
functions $g_{i}^{A}(\mathbf{r})$ are both localized and smooth, and are meant
to augment $\Delta\theta^{1}_{A,\alpha}(\mathbf{r})$ to yield the correct
values of $\Delta F_{A,\beta}(\mathbf{r})$ outside of the augmentation region.
To do this, least-squares linear regression is used to fit the following
function to $\Delta F_{A,\beta}(\mathbf{r})$:

$\displaystyle\hat{y}_{A,\beta}(\mathbf{r})$
$\displaystyle=\hat{y}^{1}_{A,\beta}(\mathbf{r})+\hat{y}^{2}_{A,\beta}(\mathbf{r})$
(84) $\displaystyle\hat{y}^{1}_{A,\beta}(\mathbf{r})$
$\displaystyle=\sum_{j}D_{j\beta}h_{j}^{A}(\mathbf{r})$ (85)
$\displaystyle\hat{y}^{2}_{A,\beta}(\mathbf{r}_{1})$
$\displaystyle=\sum_{i,\alpha}C_{i\alpha}\int\differential[3]\mathbf{r}_{2}\phi_{\alpha\beta}(r_{12})g_{i}^{A}(\mathbf{r}_{2}),$
(86)

where the coefficients $C_{i\alpha}$ and $D_{j\beta}$ are determined by the
regression. The necessary integrals are actually performed in reciprocal space
since Eq. 83 yields $\Delta F_{A,L,\beta}(k)$ rather than $\Delta
F_{A,\beta}(\mathbf{r})$. At first glance, this appears to be a complicated
and expensive fitting problem, but by separating the basis functions into
their individual angular momentum channels and doing a separate regression for
each channel, it becomes tractable. The details for this procedure are
described in Appendix B.1.

Because the functions $h_{j}^{A}(\mathbf{r})$ are localized, in the complete
basis set limit $\Delta
F_{A,\beta}(\mathbf{r})=\hat{y}^{2}_{A,\beta}(\mathbf{r})$ outside the
augmentation region. Therefore, by defining

$\tilde{\theta}^{\text{aug}}_{\alpha}(\mathbf{r})\equiv\tilde{\theta}_{\alpha}(\mathbf{r})+\sum_{A}\sum_{i}C_{i\alpha}^{A}g_{i}^{A}(\mathbf{r}),$
(87)

one can compute

$\displaystyle\tilde{G}(\mathbf{r})$
$\displaystyle=\sum_{\beta}p^{b}_{\beta}\left(\tilde{b}(\mathbf{r})\right)\tilde{F}^{\text{aug}}_{\beta}(\mathbf{r})$
(88) $\displaystyle\tilde{F}^{\text{aug}}_{\beta}(\mathbf{r})$
$\displaystyle=\sum_{\alpha}\int\differential[3]\mathbf{r}_{2}\Phi_{\alpha\beta}(r_{12})\tilde{\theta}^{\text{aug}}_{\alpha}(\mathbf{r}),$
(89)

which satisfies the condition that the pseudo-feature $\tilde{G}(\mathbf{r})$
is equal to the all-electron feature $G(\mathbf{r})$ outside the augmentation
region (within the limit of complete basis sets).

Within the core regions, $\tilde{F}^{\text{aug}}_{\beta}(\mathbf{r})$ is
projected onto the radial support grids using a set of projectors in a very
similar formalism to PAW itself. This quantity will be called
$\tilde{F}^{1,\text{aug}}_{A,\beta}(\mathbf{r})$, in keeping with the PAW
notation, and it takes the form

$\displaystyle\tilde{F}^{1,\text{aug}}_{A,\beta}(\mathbf{r})=$
$\displaystyle\sum_{j}B_{j\beta}^{A}\tilde{f}^{A}_{j}(\mathbf{r})+\tilde{F}_{A,\beta}^{1,0}(\mathbf{r})$
$\displaystyle-\sum_{j}\tilde{f}_{j}^{A}(\mathbf{r})\left[\int_{1}\differential[3]\mathbf{r}^{\prime}\tilde{p}^{A}_{j}(\mathbf{r}^{\prime})\tilde{F}_{A,\beta}^{1,0}(\mathbf{r}^{\prime})\right]$
(90) $\displaystyle B_{j\beta}^{A}=$ $\displaystyle\differential v\sum_{g\in
A}\tilde{F}^{\text{aug}}_{\beta}(\mathbf{r}_{g})\tilde{p}_{j}^{A}(\mathbf{r}_{g}).$
(91)

In the above equations, $\tilde{f}_{i}^{A}(\mathbf{r})$ and
$\tilde{p}_{j}^{A}(\mathbf{r})$ are sets of atom-centered functions satisfying
$\innerproduct{\tilde{f}_{i}^{A}}{\tilde{p}_{j}^{A}}=\delta_{ij}$, and
$\tilde{p}_{j}^{A}(\mathbf{r})$ vanish for $r>r_{c}$ so that the numerical
integral of Eq. 91 can be performed efficiently on the FFT grid points inside
the cutoff sphere of atom $A$. The symbol $\differential v$ represents the FFT
grid volume element. This is directly analagous to the projector-partial wave
pairs in PAW [61]. The construction of $\tilde{f}_{i}^{A}(\mathbf{r})$ and
$\tilde{p}_{j}^{A}(\mathbf{r})$ is described in Appendix B.2.
$\tilde{F}_{A,\beta}^{1,0}(\mathbf{r})$ is an approximate convolution of
$\tilde{\theta}^{1}_{A,\alpha}(\mathbf{r})$, and while arbitrary, it helps
make the $\tilde{f}_{j}^{A}(\mathbf{r})$ expansion more efficient. Using this
construction, the all-electron and pseudo-features on the radial support grids
are, respectively,

$\displaystyle G_{A}^{1}(\mathbf{r})$
$\displaystyle=\sum_{\beta}p^{b}_{\beta}\left(b_{A}^{1}(\mathbf{r})\right)\left(\tilde{F}^{1,\text{aug}}_{A,\beta}(\mathbf{r})+\hat{y}^{1}_{A,\beta}(\mathbf{r})\right)$
(92) $\displaystyle\tilde{G}_{A}^{1}(\mathbf{r})$
$\displaystyle=\sum_{\beta}p^{b}_{\beta}\left(\tilde{b}_{A}^{1}(\mathbf{r})\right)\tilde{F}^{1,\text{aug}}_{A,\beta}(\mathbf{r}).$
(93)

The exchange correlation energy is then

$E_{\text{xc}}[n]=E_{\text{xc}}[\tilde{n},\tilde{G}]+\sum_{A}\left(E_{\text{xc}}[n^{1}_{A},G^{1}_{A}]-E_{\text{xc}}[\tilde{n}^{1}_{A},\tilde{G}^{1}_{A}]\right),$
(94)

where all three terms are semilocal functionals of the density and of the
nonlocal features. The first term is evaluated on the FFT grid and the second
two on the radial support grids on the atoms.

The use of nonlocal features introduces a few additional terms to the forces
and stresses not present for semilocal DFT. The force terms are

$\displaystyle\mathbf{F}_{A}$
$\displaystyle=\mathbf{F}^{\text{sl}}_{A}+\mathbf{F}^{\text{PASDW}}_{A}$ (95)
$\displaystyle\mathbf{F}^{\text{PASDW}}_{A}$
$\displaystyle=\sum_{\alpha}\differential v\sum_{ig}C^{A}_{i\alpha}\nabla
g_{i}^{A}(\mathbf{r}_{g})\partialderivative{E_{\text{xc}}}{\theta_{\alpha}(\mathbf{r}_{g})}$
$\displaystyle+\sum_{\beta}\differential
v\sum_{jg}\partialderivative{E_{\text{xc}}}{B^{A}_{j\beta}}\nabla\tilde{p}_{j}^{A}(\mathbf{r}_{g})\tilde{F}_{\beta}(\mathbf{r}_{g}).$
(96)

The stress terms are

$\displaystyle\sigma_{\mu\nu}$
$\displaystyle=\sigma_{\mu\nu}^{\text{sl}}+\sigma_{\mu\nu}^{\text{PASDW}}$
(97) $\displaystyle\sigma_{\mu\nu}^{\text{PASDW}}$
$\displaystyle=\sum_{\alpha}\sum_{A}\differential
v\sum_{g}\left\\{\left[\left(\mathbf{r}_{g}-\mathbf{R}_{A}\right)_{\mu}\partialderivative{E_{\text{xc}}}{\theta_{\alpha}(\mathbf{r}_{g})}\right]\left[\sum_{i}C_{i\alpha}^{A}\nabla_{\nu}g_{i}^{A}(\mathbf{r}_{g})\right]\right\\}+\delta_{\mu\nu}\sum_{\beta}\sum_{A}\sum_{j}\left(\partialderivative{E_{\text{xc}}}{B_{j\beta}^{A}}B_{j\beta}^{A}\right)$
$\displaystyle+\sum_{\beta}\sum_{A}\differential
v\sum_{g}\left\\{\left[(\mathbf{r}_{g}-\mathbf{R}_{A})_{\mu}F_{\beta}(\mathbf{r}_{g})\right]\left[\sum_{j}\partialderivative{E_{\text{xc}}}{B_{j\beta}^{A}}\nabla_{\nu}\tilde{p}_{j}^{A}(\mathbf{r}_{g})\right]\right\\}$
$\displaystyle+\sum_{k}\frac{G_{k\mu}G_{k\nu}}{|\mathbf{G}_{k}|}\sum_{\beta}\left(\partialderivative{E_{\text{xc}}}{\tilde{F}_{\beta}(\mathbf{G}_{k})}\right)^{*}\sum_{\alpha}\Phi_{\alpha\beta}^{\prime}(|\mathbf{G}_{k}|)\theta_{\alpha}(\mathbf{G}_{k}).$
(98)

In the above equations, $\mathbf{F}^{\text{sl}}_{A}$ and
$\sigma_{\mu\nu}^{\text{sl}}$ refer to the force and stress equations for
semilocal functionals within the PAW formalism (see Section IIID of Ref.
[68]).

There are several significant contributions to the computational cost of the
above approach for evaluating $G(\mathbf{r})$, which all scale linearly or
quasi-linearly with system size. The first contribution comes from computing
the Fourier transforms of $\tilde{\theta}_{\alpha}^{\text{aug}}(\mathbf{r})$
and $\tilde{F}^{\text{aug}}_{\beta}(\mathbf{r})$ (as well as the functional
derivatives with respect to these quantities), which is necessary to compute
the convolutions like those in Eq. 89 in reciprocal space. The cost of these
operations scales as
$\Theta(N_{\alpha}\,N_{\text{FFT}}\,\log(N_{\text{FFT}}))$, where $N_{\alpha}$
is the number of kernel interpolation points in Eq. 59 and $N_{\text{FFT}}$ is
the size of the FFT grid used for integrating the XC energy. The second
contribution arises from computing the convolutions in Eq. 89 in reciprocal
space, and scales as $\Theta(N_{\alpha}^{2}\,N_{\text{FFT}})$. The third
contribution comes from computing the PAW correction to
$\tilde{\theta}_{\alpha}^{\text{aug}}(\mathbf{r})$ (the sum in Eq. 87). This
and related terms scale as
$\Theta(N_{\text{FFT,atom}}\,N_{\text{atom}}\,N_{\alpha})$, where
$N_{\text{FFT,atom}}$ is the number of FFT grid points inside the core region
of an atom and $N_{\text{atom}}$ is the number of atoms. The last contribution
consists of all the other PAW correction-related routines (e.g. computing
$C_{i\alpha}^{A}$), which depend only on the density on the atomic radial
grids and therefore scale linearly with $N_{\text{atom}}$. The overall scaling
with system size $N$ is then $\Theta(N\log(N))$. Notably, none of the
components of the algorithm scale with the number of k-points.

##  III Results and Discussion

Figure 1: Accuracy of PBE0/CIDER (Eq. 100) compared to PBE0 for each
functional type on the GMTKN55 database of molecular properties, using the MoM
error metric. The database subsets are [10]: “all”: Full database; “small”:
Basic properties and reaction energies for small systems; “big”: Reaction
energies for large systems and isomerization reactions; “barrier”: Barrier
heights; “inter NCI”: Intermolecular noncovalent interactions; “intra NCI”:
Intramolecular noncovalent interactions.

Using Gaussian process regression, we fit each of the model types in Section
II.3 (SL-GGA, NL-GGA, SL-MGGA, and NL-MGGA) to match the exact exchange energy
on subsets of the GMTKN55 [10] and SOL62 [69, 33] databases, as well as a few
dozen other training points. We provide a summary of our methods in this
section and explain them in more technical detail in Appendix A.

To train the ML functionals, ground-state DFT calculations were first
performed using the PBE functional [23]. Then, for each chemical system in the
training set, the density and orbitals obtained from the PBE calculation were
used to compute the features described in Section II.3.3, as well the exact
exchange energy of Eq. 2. The computed features were used as the inputs to the
Gaussian process to predict the exact exchange energy. See Appendices A.2 and
A.3 for details.

The training and validation systems were chosen as described in Appendix A.4.
The SOL62 database [69] contains cohesive energies for a mix of metals and
non-metals, for which 9 were used for training, 22 for validation, and 31 for
testing. The GMTKN55 database contains 55 sub-databases blocked into five
categories: small-molecule properties (“small”), larger molecule reactions and
isomerizations (“big”), barrier heights (“barrier”), intermolecular
noncovalent interactions (“inter NCI”), and intramolecular noncovalent
interactions (“intra NCI”). The training and validation data consist of twelve
sub-databases of GMTKN55 and contain some representative systems from each
category, but disproportionately taken from the “small” category. All data not
used for training or validation is considered the GMTKN55 test set. This
GMTKN55 partitioning results in 280 train data, 163 validation data, and 1062
test data. Further technical details for training are provided in Appendices
A.6 and A.7.

The SL-GGA, NL-GGA, SL-MGGA, and NL-MGGA functionals used for benchmarking
below are the functionals of each type with the best validation set
performance, as described in Appendix A.8. One important technical note is
that each functional was trained via delta-learning on top of some baseline
exchange functional. We tried PBE and the Chachiyo GGA [70] as baselines, and
Chachiyo proved more accurate for each functional class. However, as described
in Section III.3 below, slight numerical stability improvements were found
when using the NL-MGGA with a PBE baseline (dubbed NL-MGGA-PBE). As described
in more detail in Section III.3, we also found that retraining NL-MGGA-PBE
with a more diverse training set from GMTKN55 (which is described in Appendix
A.5) results in improved numerical stability. Therefore, the resulting
functional, dubbed NL-MGGA-DTR, was used for the applications in Sections
III.4-III.7.

When testing and benchmarking the functionals, all calculations (both
conventional and CIDER functionals) were performed self-consistently, with the
exception of the SOL62 cohesive energies. These were obtained non-self-
consistently from PBE orbitals, due to the difficulty of computing the
isolated atom energies in plane-wave DFT. This is a small technical barrier
that will be addressed in future work. All calculations in the main text used
static geometries, with the exception of the defect calculations in Section
III.7; in this section, the silicon lattice parameters and defect structures
were optimized with both PBE and CIDER.

In Section III.1, we demonstrate that the NL-MGGA is capable of learning the
exchange functional more accurately than the other feature sets, and in
Section III.2 we use 1 and 2-electron systems to provide intuition for why
this is. Next, we fine-tune the accuracy and numerical stability of the NL-
MGGA (Section III.3) and demonstrate that the resulting exchange functional
can be used to match the accuracy of hybrid DFT (Section III.4).

Having designed an accurate exchange functional, we demonstrate that it can be
used to improve band gap predictions over standard meta-GGAs in Section III.5.
In Section III.6, we show that for large systems in plane-wave DFT, our
implementation of the NL-MGGA form is similarly efficient to semilocal meta-
GGAs and orders of magnitude more efficient than hybrid DFT. Lastly, in
Section III.7, we take advantage of our new functional’s efficiency and
accurate band gap predictions to compute charge transition levels of point
defects in silicon, which is a notoriously challenging task for DFT due to the
underestimation of band gaps and the large supercell sizes required to
simulate isolated point defects [6].

###  III.1 Accuracy of Different Model Types for Molecular and Solid-state
Benchmarks

We assess the accuracy of our exchange functionals by substituting them for
exact exchange in hybrid functionals to create a “surrogate hybrid”
functional. For example, one of the simplest and most popular hybrid
functionals is PBE0 [71], which mixes a 1/4 fraction of exact (Hartree-Fock,
HF) exchange into the PBE functional [23]:

$E_{\text{xc}}[n]=\frac{3}{4}E_{\text{x}}^{\text{PBE}}[n]+\frac{1}{4}E_{\text{x}}^{\text{HF}}[n]+E_{\text{c}}^{\text{PBE}}[n].$
(99)

In this section we perform so-called PBE0/CIDER calculations, where the
exchange-correlation energy is

$E_{\text{xc}}[n]=\frac{3}{4}E_{\text{x}}^{\text{PBE}}[n]+\frac{1}{4}E_{\text{x}}^{\text{CIDER}}[n]+E_{\text{c}}^{\text{PBE}}[n].$
(100)

One can also generalize PBE0/CIDER with different fractions of CIDER exchange,
giving the PBE0($\alpha$)/CIDER form, which is used in later sections:

$E_{\text{xc}}[n]=(1-\alpha)E_{\text{x}}^{\text{PBE}}[n]+\alpha
E_{\text{x}}^{\text{CIDER}}[n]+E_{\text{c}}^{\text{PBE}}[n].$ (101)

Throughout the results, we will also use abbreviations like PBE0/SL-GGA and
PBE0/NL-MGGA, which indicate that the specific CIDER functionals used in
PBE0/CIDER are the SL-GGA and NL-MGGA, respectively.

All ML functional results in this section refer to the PBE0/CIDER surrogate
hybrid form. Because our exchange functionals are trained on exact exchange
energies, they are “perfect” if they exactly reproduce the predicted reaction
energies of PBE0 on the test set. To measure the error of the exchange
functionals, we calculate the average error between PBE0/CIDER and PBE0 for
each functional type on GMTKN55 and SOL62. All calculations on GMTKN55 are SCF
calculations, but the SOL62 cohesive energy CIDER and PBE0 calculations are
performed non-self-consistently with PBE orbitals. See Appendices A.2 and A.3
for details.

The average deviations of the machine-learned functionals from PBE0, using the
mean of means (MoM) metric (i.e. the mean of the mean absolute deviations, or
MADs, for the sub-databases of GMTKN55), are presented in Fig. 1. Results for
all functionals include D4 dispersion corrections [72], with the PBE0/CIDER
functionals using the PBE0 dispersion parameters. The nonempirical PBE and
r2SCAN [73] functionals, while not designed to reproduce PBE0, are included in
Fig. 1 to show typical deviations between GGAs, meta-GGAs, and hybrids.

Figure 2: Average deviation of PBE and the PBE0/CIDER ML functionals from PBE0
on the SOL62 cohesive energies database.

While the trends vary slightly in the different subsets of GMTKN55, the key
trend is that the quality of the fit to exact exchange improves going from SL-
GGA to NL-GGA to SL-MGGA to NL-MGGA. Notably, the SL-GGA and NL-GGA suffer
from out-of-distribution transferability issues, resulting in larger
deviations from PBE0 even compared to the PBE functional on the test set. On
the other hand, the ML meta-GGAs are more transferable outside the training
data than the ML GGAs, so they match PBE0 more closely than PBE and r2SCAN on
the test set. The NL-MGGA is more accurate than the SL-MGGA, achieving an MoM
of less than 1 kcal/mol on the test set. The fact that the NL-MGGA outperforms
both the NL-GGA and SL-MGGA suggests that the kinetic energy and nonlocal
features provide non-redundant information, meaning that nonlocal density
features cannot replace the kinetic energy density nor vice versa. As shown in
Fig. 2, the same general trends apply for SOL62; the NL-MGGA is the most
accurate and transferable ML functional for reproducing PBE0, followed by the
SL-MGGA and then the NL-GGA.

Because GMTKN55 contains a variety of properties and systems, the performance
of the functionals is somewhat dependent on the choice of error metric.
However, the general trends are the same as for the MoM metric. See
Supplemental Material Section S2A for a discussion of the impact of error
metric on the results.

###  III.2 Impact of Nonlocality and Derivative Discontinuity on Model
Accuracy

Figure 3: Geometries and density distributions of the 1-electron systems
dataset. a) Hydrogen atom. b) H2+. c) H3^2+-linear. d) H3^2+-triangular.

Overall, the SL-MGGA model is more accurate than the NL-GGA model across the
databases tested here. However, the fact that NL-GGA improves over SL-GGA and
that NL-MGGA improves over SL-MGGA suggests that the nonlocal features
contribute useful, non-redundant information about the density for exchange
energy prediction. This raises the question: What types of systems require the
nonlocal features, and which require the semilocal orbital dependence of the
kinetic energy density? We hypothesize that the derivative discontinuity of
the exchange functional [74] plays a key role. Because only the kinetic energy
density—and not the nonlocal features of the NL-GGA—includes a derivative
discontinuity passing through integer electron number, the NL-GGA (as well as
any smooth, pure density functional) cannot be expected to accurately capture
energy differences involving a change in electron number, such as ionization
energies and electron affinities. On the other hand, existing semilocal meta-
GGAs are fundamentally limited by their locality and cannot be expected to
improve over GGAs for energy differences between systems of like electron
number. This is especially true of one-electron systems, for which the kinetic
energy is simply determined by the density gradient.

To explore this hypothesis, we set up two simple sets of systems. The first
set is one-electron systems for different arrangements of protons: the
hydrogen atom, the H2+ molecular ion, the linear H3^2+ molecular ion, and the
triangular H3^2+ molecular ion, all with bond lengths of 0.7 Å, as illustrated
in Fig. 3. Because each system contains only one electron with different
degrees of delocalization, we expect that only the functionals with nonlocal
features can accurately predict energy differences between these systems.
Table 1 shows the MAD between the ML exchange functionals and Hartree-Fock for
the energy differences between the hydrogen atom and the molecular ions. As
expected, the SL-GGA and SL-MGGA are both less accurate than NL-GGA and NL-
MGGA. This confirms that the nonlocal features successfully characterize
single-electron delocalization in a way that is not feasible with semilocal
functionals.

Table 1: Mean absolute deviation (kcal/mol) of the four types of ML functional for simple 1 and 2-electron systems compared to Hartree-Fock. Model | 1e Systems111The 1-electron system dataset consists of the energies of H2+, H3^2+-linear, and H3^2+-triangular referenced to the energy of the isolated hydrogen atom. All nearest-neighbor bond lengths are 0.7 Å. | 2e Systems222The 2-electron system dataset consists of three energy differences between He+, He-singlet, and He-triplet, where He-triplet is the neutral helium atom excited state with both electrons having the same spin.
---|---|---
SL-GGA | $2.9$ | $8.7$
NL-GGA | $2.5$ | $5.5$
SL-MGGA | $4.8$ | $3.9$
NL-MGGA | $1.1$ | $2.7$

The second set of systems consists of the He+ ion, the He singlet (closed-
shell ground state), and the He triplet (open-shell excited spin state with
two same-spin electrons). Because the exchange energy has a discontinuous
derivative as the number of electrons of a given spin passes through an
integer [75, 74], the energy differences between He+ and He triplet and
between He singlet and He triplet can only be described accurately by
functionals that capture this derivative discontinuity. As shown in Table 1,
the absolute deviations from Hartree Fock for the energy differences between
these three systems are larger for SL-GGA and NL-GGA than for SL-MGGA and NL-
MGGA. This indicates that the derivative discontinuity included in the kinetic
energy density is essential for capturing energy differences related to
electron addition/removal. Notably, the NL-MGGA is the most accurate of the
four functional types for both datasets.

These observations regarding the role of the derivative discontinuity explain
why a previous attempt at machine learning a GGA resulted in more accurate
atomization energies but less accurate ionization potentials than PBE [31].
Likewise, the role of nonlocality explains why Kovács _et al._ [13] found that
existing meta-GGAs, as well as 25 meta-GGAs they empirically trained, suffer
from a trade-off between cohesive energy accuracy and band gap accuracy. Meta-
GGAs are limited in their ability to describe both these properties at once
due their semilocality. In addition, while we do not exhaustively try many
possible NLDFs to conclusively show that they can only yield accurate models
when paired with orbital-dependent features, the NLDFs we have implemented do
need to be paired with a meta-GGA formalism to improve their accuracy. These
findings provide evidence that to accurately and transferably fit a nonlocal,
orbital-dependent quantity like the exchange energy, an ML functional should
contain both orbital dependence and nonlocality in its feature set.

###  III.3 Fine-tuning for Numerical Stability

Table 2: The number of systems in the GMTKN55 database for which a given functional did not achieve SCF convergence with thresholds of $10^{-8}$ Eh for total energy and $10^{-4}$ Eh for the orbital gradient of the energy. All functionals use the PBE0/CIDER form (Eq. 100). Functional | No. Unconverged
---|---
SL-GGA | 10
NL-GGA | 16
SL-MGGA | 5
NL-MGGA | 3
NL-MGGA-PBE | 1
NL-MGGA-DTR | 0

Because XC functionals are used within SCF calculations, they must be
numerically stable and consistently achieve electronic convergence on real
systems to be useful for practical applications. For the GMTKN55 dataset, we
found that some systems do not quite converge to our default convergence
thresholds, which were $10^{-8}$ Hartree atomic unit (Eh) for the total energy
and $10^{-4}$ Eh for the orbital gradient of the energy. (See Appendix A.9 for
more details on how we remedied convergence problems.) We found that these
calculations could only be converged with a $10^{-7}$ Eh energy threshold and
$2\times 10^{-3}$ Eh gradient threshold. Interestingly, this problem is least
severe for the NL-MGGA, as illustrated in Table 2, which shows the number of
systems for each functional type which did not converge to the $10^{-8}$ Eh
energy and $10^{-4}$ Eh gradient thresholds. We believe this is because the
increased model dimensionality of the NL-MGGA allows it to fit the exchange
energy with a smoother function of the input features.

Figure 4: Accuracy of GMTKN55 predictions for the NL-MGGA and the NL-MGGA-PBE
and NL-MGGA-DTR variants described in the text (all using the PBE0/CIDER
form), compared to PBE0. Train/test partition results are not shown because
the train and test data are different for the NL-MGGA-DTR functional.

Because the NL-MGGA still sees some convergence problems, we decided to
explore two modifications of the NL-MGGA for improved numerical stability. The
first is the NL-MGGA-PBE, which has the same hyperparameters and feature shape
as the NL-MGGA but uses the PBE functional as a baseline rather than the
Chachiyo functional [70]. NL-MGGA-PBE was screened during the validation
procedure (Appendix A.8) but was found to be slightly less accurate than the
NL-MGGA with the Chachiyo baseline. However, only 1 system sees a convergence
issue in Table 2 with NL-MGGA-PBE, compared to 3 with NL-MGGA. We also tried
training an NL-MGGA model with the PBE baseline and a more chemically diverse
subset of GMTKN55, as described in Appendix A.5. We refer to this new
train/validation/test partition as DTR (for “diverse training”), and it yields
the NL-MGGA-DTR, for which all GMTKN55 systems achieve an orbital gradient
convergence of $10^{-4}$ Eh or lower. In fact, all but 10 of the 2462 GMTKN55
systems converge to within $10^{-9}$ Eh for energy and $10^{-4.5}$ Eh for
orbital gradient with NL-MGGA-DTR, and the rest to within $10^{-8}$ and
$10^{-4}$ Eh, respectively.

As shown in Fig. 4, which displays MoM errors of PBE0/CIDER versus PBE0, the
accuracy of these three NL-MGGA variants is similar across the GMTKN55
dataset. The three functional variants also have nearly identical performance
for the SOL62 cohesive energies, with MADs compared to PBE0 of 0.13 eV/atom
for all three. Therefore, for the remainder of this work, we adopt the NL-
MGGA-DTR model as our primary functional of choice for applications, and we
name it CIDER23X-NL-MGGA-DTR.

###  III.4 Obtaining Hybrid DFT Accuracy

Figure 5: MoM error for all of GMTKN55 (left) and MAD on the DTR test set
(right) relative to GMTKN55 reference values, for semilocal, hybrid, and CIDER
functionals.

In this section, we show that by substituting NL-MGGA-DTR for exact exchange
in hybrid functionals, accurate predictions of the XC energy can be obtained
for molecular systems without additional parameter tuning. In particular, we
look at PBE0/CIDER (Eq. 100) and another type of functional called
PW6B95/CIDER. The original PW6B95 functional [76] is

$E_{\text{xc}}^{\text{PW6B95}}[n]=0.28E_{\text{x}}^{\text{HF}}[n]+E_{\text{xc,sl}}^{\text{PW6B95}}[n],$
(102)

where $E_{\text{xc,sl}}^{\text{PW6B95}}$ is a meta-GGA term. Just like with
PBE0/CIDER, we substitute CIDER in for Hartree-Fock exchange to obtain
PW6B95/CIDER:

$E_{\text{xc}}^{\text{PW6B95}}[n]=0.28E_{\text{x}}^{\text{CIDER}}[n]+E_{\text{xc,sl}}^{\text{PW6B95}}[n].$
(103)

We chose to construct PW6B95/CIDER functionals because the original PW6B95
[76] is among the most accurate functionals on the GMTKN55 database [10, 36].

For PBE0/NL-MGGA-DTR, PW6B95/NL-MGGA-DTR, and several semilocal [23, 77, 73,
78] and hybrid [27, 71, 76, 28, 79] functionals, Fig. 5 shows the MoM error on
all of GMTKN55 as well as the mean absolute error (MAE) on the test partition
used for NL-MGGA-DTR. The errors are relative to the high-accuracy benchmark
data collected in the GMTKN55 database [10]. All functionals include D4
dispersion corrections [72] except for B97M-V [78], $\omega$B97X-V [28], and
$\omega$B97M-V [79], which all use VV10 dispersion [54]. PBE0/CIDER
calculations use PBE0 dispersion parameters, and PW6B95/CIDER calculations use
PW6B95 dispersion parameters. We show conventional semilocal functionals in
shades of blue, conventional hybrid functionals in orange, and CIDER
functionals in green.

As Fig. 5 shows, the PBE0/NL-MGGA-DTR functional is comparable in accuracy to
r2SCAN, but the combination of systematic errors in the NL-MGGA-DTR functional
and in PBE0 itself prevent it from outperforming the most accurate meta-GGA
(B97M-V) or any of the hybrid functionals. On the other hand, because PW6B95
is much more accurate than PBE0 on the GMTKN55 dataset, PW6B95/NL-MGGA-DTR is
more accurate than all semilocal functionals explored here. It is also more
accurate than the PBE0 and B3LYP [27] hybrid functionals and equivalent in
accuracy to the range-separated hybrid $\omega$B97X-V [28]. This shows that
the NL-MGGA-DTR exchange functional can be used to calculate molecular
reaction energies with an accuracy that was previously not feasible without
exact exchange mixing.

As with Section III.1, the analysis of the errors on GMTKN55 depends somewhat
on the error metric used. These nuances are discussed in further detail in
Supplemental Material Section S2B, which also provides a comparison between
the nonlocal and semilocal meta-GGA models.

###  III.5 Using CIDER for Improved Band Gap Prediction

Table 3: Error statistics for 255 band gaps from the database published by Borlido _et al._ [42], in eV. The ML functionals use the PBE0($\alpha$)/CIDER formalism of Eq. 101. | PBE333This work. | r2SCAN11footnotemark: 1 | SL-MGGA11footnotemark: 1 | NL-MGGA-DTR11footnotemark: 1 | HSE06444Borlido _et al._ [42] | PBE022footnotemark: 2
---|---|---|---|---|---|---
|  |  | $\alpha=0.25$ | $\alpha=0.25$ | $\alpha=0.35$ |  |
ME555Mean error (eV) | $-1.14$ | $-0.81$ | $-0.74$ | $-0.61$ | $-0.38$ | $-0.17$ | $0.40$
MAE666Mean absolute error (eV) | $1.16$ | $0.87$ | $0.82$ | $0.71$ | $0.58$ | $0.53$ | $0.73$
RMSE777Root mean square error (eV) | $1.50$ | $1.18$ | $1.10$ | $1.02$ | $0.86$ | $0.80$ | $0.90$
MARE888Mean absolute relative error | $0.53$ | $0.41$ | $0.41$ | $0.35$ | $0.33$ | $0.33$ | $0.67$

The band gap problem [80, 81], in which semilocal DFT tends to drastically
underestimate the band gaps of solids, stands in the way of a variety of
important applications like high-throughput screening for optical and
electronic properties and semiconductor point defect physics [5, 6, 82]. By
introducing a derivative discontinuity into the XC functional, meta-GGAs
partly mitigate this problem [83], but they still tend to significantly
underestimate band gaps [84]. Exceptions include mBJ [85] and TASK [12], which
each have important drawbacks. The mBJ functional provides the exchange
potential only, not the energy [85], and TASK does not predict atomization
energies and lattice constants as well as SCAN [12, 13, 14]. This reinforces
the discussion in Section III.2 regarding the accuracy trade-offs meta-GGAs
must make due to their semilocality.

By introducing nonlocality into the exchange energy, PBE0/NL-MGGA-DTR (Eq. 100
or Eq. 101 with $\alpha=0.25$) improves band gap predictions over PBE and
r2SCAN, as shown in Table 3. This improvement occurs even though the CIDER
functionals are not fit to solid-state band gaps. The PBE0/SL-MGGA also
improves band gap predictions over r2SCAN, but significantly less so than
PBE0/NL-MGGA-DTR. The database used for testing covers 255 of the 472 band
gaps in the database published by Borlido _et al._ [42], with La, Yb, and Th-
containing systems being excluded as well as systems with a greater than
$0.05$ eV difference between the spin-orbit-corrected and non-spin-orbit-
corrected band gaps (as computed with PBE by Borlido _et al._ [42]). See
Appendix A.11 for details, and see Supplemental Material Section S3 for error
metrics for all 453 solids in the database that do not contain La, Yb, or Th.

To provide a clearer picture of the improved band gap prediction provided by
the NL-MGGA-DTR functional, Table 4 shows the band gap predictions for a
selected subset of systems that are commonly studied [11]. In addition, Figs.
6 and 7 contain band structures for silicon (Si) and boron phosphide (BP),
respectively, which show that PBE0/NL-MGGA-DTR provides similar band
structures to PBE and r2SCAN but with a wider band gap.

Table 4: Band gaps for a small subset of solids (in eV). CIDER refers to PBE0/NL-MGGA-DTR. | PBE | r2SCAN | CIDER | HSE06 | Expt. [42]
---|---|---|---|---|---
Si | $0.57$ | $0.76$ | $1.05$ | $1.15$ | $1.17$
Ge | $0.12$ | $0.48$ | $0.71$ | $0.81$ | $0.74$
InP | $0.69$ | $1.12$ | $1.47$ | $1.52$ | $1.42$
GaAs | $0.59$ | $1.08$ | $1.36$ | $1.43$ | $1.52$
CdSe | $0.71$ | $1.19$ | $1.48$ | $1.69$ | $1.74$
BP | $1.25$ | $1.42$ | $1.71$ | $1.98$ | $2.10$
GaP | $1.61$ | $1.88$ | $2.09$ | $2.29$ | $2.35$
CdS | $1.20$ | $1.65$ | $2.02$ | $2.27$ | $2.48$
GaN | $1.88$ | $2.33$ | $2.65$ | $3.18$ | $3.50$
ZnS | $2.13$ | $2.69$ | $3.04$ | $3.37$ | $3.72$
C | $4.13$ | $4.34$ | $4.70$ | $5.33$ | $5.50$
BN | $4.21$ | $4.79$ | $5.08$ | $5.63$ | $5.96$
CaO | $3.68$ | $4.23$ | $4.53$ | $5.32$ | $6.88$
MgO | $4.73$ | $5.63$ | $6.11$ | $6.42$ | $7.67$
NaCl | $5.11$ | $5.90$ | $6.32$ | $6.45$ | $8.75$
LiF | $9.07$ | $10.05$ | $10.40$ | $11.39$ | $13.60$
Ar | $8.71$ | $9.61$ | $9.71$ | $10.36$ | $14.15$
ME | $-1.93$ | $-1.42$ | $-1.11$ | $-0.74$ | -
MAE | $1.93$ | $1.42$ | $1.11$ | $0.76$ | -
RMSE | $2.41$ | $1.88$ | $1.66$ | $1.31$ | -
MARE | $0.45$ | $0.29$ | $0.18$ | $0.11$ | -

(a) fig:pbe_si_bs (b) fig:r2scan_si_bs (c) fig:cider_si_bs

Figure 6: Band structures of silicon with three different functionals. The
black dotted line is positioned at the valence band maximum, and the orange
dotted line is located at the experimental band gap of 1.17 eV.

(a) fig:pbe_bp_bs (b) fig:r2scan_bp_bs (c) fig:cider_bp_bs

Figure 7: Band structures of boron phosphide with three different functionals.
The black dotted line is positioned at the valence band maximum, and the
orange dotted line is located at the experimental band gap of 2.10 eV.

There is room for improvement in these band gaps, since PBE0/NL-MGGA-DTR still
systematically underestimates band gaps compared to the state-of-the-art HSE06
[86, 87] range-separated hybrid functional, as well as the PBE0 functional to
which PBE0/CIDER is an approximation. One possible solution is to tune the
fraction of CIDER exchange, which is already done with exact exchange within
hybrid DFT to optimize functional accuracy for different applications [88]. As
shown in Table 3, PBE0($\alpha$)/NL-MGGA-DTR with $\alpha=0.35$ has increased
band gaps over $\alpha=0.25$ and similar error statistics to HSE06. However,
one drawback of this approach is that because NL-MGGA-DTR increases band gaps
less than exact exchange for a given mixing fraction, one might need a large
fraction of NL-MGGA-DTR to fit band gaps. This large fraction might not
provide accurate predictions of other material properties. Alternative
approaches to improve CIDER band gaps include using physical intuition and
exact constraints to tune the derivative discontinuity contribution from the
kinetic energy density (as is done in the TASK functional [12]), introducing
more nonlocal features to improve the accuracy of the functional in general,
and training on fundamental band gaps (i.e. differences between ionization
potentials and electron affinities). However, considering the drastic increase
in computational cost associated with using HSE06 or PBE0 compared to
semilocal DFT, the ability to mitigate the underestimation of band gaps
without explicitly fitting them—and while accurately predicting other
properties like bond energies—already serves as a major step toward solving
the band gap problem.

###  III.6 Benchmarking Performance on Large Condensed Matter Systems

(a)  (b)

Figure 8: Computational cost for plane-wave DFT SCF calculations of large systems with different functionals. All semilocal (blue) and CIDER (green) calculations were performed with GPAW on one 64-core Intel Ice Lake node for diamond and two 64-core nodes for silicon. The PBE0 CPU calculation (bright orange) was performed on one 64-core node, while the GPU calculation (light orange) was performed on one Nvidia A100 GPU. a) Diamond 216-atom bulk supercell with $\Gamma$-point sampling only. Note the logarithmic time axis. b) Diamond and silicon neutral vacancy formation energies with 216-atom supercell with $2\times 2\times 2$ Monkhorst-Pack k-point mesh [89]. Table 5: The formation energies of the neutral vacancies in diamond and silicon (in eV), computed with 216-atom supercells and a $2\times 2\times 2$ Monkhorst-Pack k-point mesh [89]. All CIDER functionals use the PBE0/CIDER surrogate hybrid form. Functional | C | Si
---|---|---
PBE | $6.47$ | $3.61$
r2SCAN | $6.36$ | $4.19$
SL-GGA | $6.11$ | $3.48$
NL-GGA | $6.73$ | $3.50$
SL-MGGA | $6.66$ | $4.16$
NL-MGGA | $6.77$ | $4.58$
NL-MGGA-PBE | $6.77$ | $4.58$
NL-MGGA-DTR | $6.79$ | $4.42$
HSE06 | $6.96$ | $4.54$
PBE0 | $6.97$ | $4.65$

The key benefit of the nonlocal features used to learn the exchange functional
is that the quasi-linear-scaling algorithm used to compute them within plane-
wave DFT (Sections II.4 and II.6) is over an order of magnitude faster than
the evaluation of exact exchange for large systems and does not depend on the
number of k-points. To illustrate this efficiency, Fig. 8a shows the wall time
to perform a ground-state SCF calculation for a 216-atom diamond supercell
with semilocal and PBE0/CIDER functionals in GPAW [62, 63]. These wall times
are compared to those of the hybrid DFT implementation in JDFTx [90], which
uses the adaptively compressed exchange algorithm [17] to compute the exact
exchange energy as efficiently as possible. The CIDER functionals, including
those with nonlocal features, all take less than 170 seconds (compared to 60
seconds for PBE and 79 seconds for r2SCAN). When run on the same number of
processors, PBE0 is roughly 700 times slower than the NL-MGGA-DTR functional.
Even the GPU-accelerated evaluation of PBE0 is about 40 times slower than NL-
MGGA-DTR.

The efficiency of CIDER functionals for large systems makes it practical to
perform more complicated calculations like defect formation energies. Figure
8b shows the wall time for computing the neutral vacancy formation energies in
diamond and silicon using 216-atom supercells for semilocal and PBE0/CIDER
functionals. All of the functionals are of comparable computational cost, with
the NL-MGGA-DTR model being only 25% slower than r2SCAN for diamond and 50%
slower for silicon. This small decrease in efficiency comes with a significant
benefit: As shown in Table 5, of the functionals tested, the NL-MGGA model and
its variants give the closest match to PBE0 for predicting these vacancy
formation energies, with PBE0/NL-MGGA-DTR giving deviations of 0.18 eV and
0.23 eV for diamond for silicon, respectively. The HSE06 and PBE0 formation
energies were computed using VASP [91, 92, 93]. See Appendix A.12 for further
calculation details.

###  III.7 Charged Defect Transition Levels in Silicon

(a) Freysoldt finite-size corrections [94].

(b) Potential alignment correction only.

Figure 9: Electronic transitions levels for the single vacancy and several
substitutional defects in silicon computed using different functionals and
finite-size correction schemes. The grey bars represent the band gap for the
given functional (or experiment). The experimental references are from the
following sources: vacancy [95], boron [96], phosphorus [97], copper [98], and
sulfur [99, 100]. CIDER refers to PBE0(0.3)/NL-MGGA-DTR. Figure 10: Change in
charge distribution in the $(001)$ plane (integrated over the $z$ axis) for
several charge transition levels in point defects of silicon. Blue shades
represent an increase in electron density, while red shades indicate a
decrease in electron density. “MaxDisp” is the maximum ionic displacement
between the charge states, and “Loc.” is the fraction of the charge difference
contained within a 5 Å sphere around the defect site. Loc. is greater than 1
for the vacancy ++/0 transition because more than 2 electrons are drawn in to
the defect region upon addition of an electron pair, depleting electron
density in the surrounding region.

One of the consequences of the band gap problem is that transition levels in
charged defects are not accurately described by semilocal DFT, so hybrid DFT
is typically used when accurate transition level predictions are needed. In
this section, we show that by improving band gap prediction, CIDER is able to
accurately predict defect transition levels in silicon. The transition levels
associated with the single vacancy [95], as well as the boron [96], phosphorus
[97], copper [98], and sulfur [99, 100] substitutionals, were computed as
described in Appendix A.13, using a 512-atom supercell, $2\times 2\times 2$
$\Gamma$-centered k-point mesh, and full structural relaxation for each defect
and functional.

Figure 9 shows the computed transition levels compared to experiment with two
different finite-size correction schemes: the Freysoldt finite-size
corrections [94] (Fig. 9a) and the potential alignment correction (Fig. 9b).
We computed the transition levels with PBE and PBE0(0.3)/NL-MGGA-DTR (Eq. 101
with $\alpha=0.3$). The latter functional was chosen because its band gap
matches the 1.17 eV experimental gap; however, PBE0/NL-MGGA-DTR yields similar
results (see Supplemental Material Section S7). Supplemental Table S13
provides a more detailed list of all computed transition levels. In the rest
of this section, PBE0(0.3)/NL-MGGA-DTR is referred to simply as CIDER for
brevity.

PBE and CIDER both give reasonable descriptions of the ++/0 transition level
in the silicon vacancy and the boron acceptor level. However, the phosphorus
donor level is not well-described by PBE due to its severe underestimation of
the band gap; CIDER remedies the band gap underestimation and therefore
predicts a transition level in good agreement with experiment. To correct for
band gap underestimation, one could predict defect levels from PBE by scaling
the gap (and with it, the transition levels) to the experimental gap, but this
makes it difficult to distinguish true shallow levels from deep levels that
happen to sit near the (too-low) PBE conduction band edge [5]. For example,
PBE predicts a transition level of 0.66 eV for the -/– transition of
$\text{Cu}_{\text{Si}}$ (Fig. 9a), 0.05 eV above the conduction band minimum.
Without further information, it is unclear whether this state is resonant in
the conduction band or localized in the gap. CIDER, on the other hand, places
the -/– level inside the gap at 0.97 eV, in good agreement with the 1.00 eV
level from experiment.

More generally, Fig. 9a shows that paired with the Freysoldt correction, CIDER
matches each level in the figure with an error of 0.10 eV or less compared to
experiment, with the exception of the +/0 level in $\text{Cu}_{\text{Si}}$.
For this transition, the experimental level is 0.20 eV [98], while the CIDER
level is 0.02 eV with the Freysoldt correction and 0.07 eV with the potential
alignment correction only. The HSE06 level was previously reported to be 0.21
eV by Sharan _et al._ [101], in good agreement with experiment. However,
Sharan _et al._ used a 64-atom supercell with a $2\times 2\times 2$ Monkhorst-
Pack k-point mesh [89]. This supercell size is relatively small, and the
k-point mesh is both small and fails to sample the band edges, which can
interfere with the description of relatively delocalized defect states. When
we recalculated the CIDER transition level using Sharan _et al._ ’s supercell
size and k-point mesh, we obtained 0.18 eV with Freysoldt corrections, in good
agreement with their result and with experiment. In the Supplemental Material
(Section S7A), we provide further evidence that the HSE06 level itself is not
converged with respect to k-point mesh. This finding illustrates the
importance of using well-converged supercell sizes and k-point meshes in point
defect calculations, a requirement that is made much easier by the efficiency
of the CIDER functionals compared to hybrid DFT. Why this particular
transition level is underestimated by both hybrid DFT and CIDER is an
interesting question warranting future investigation.

In addition to yielding more accurate transition levels, CIDER predicts
somewhat different charge distributions than PBE as well. Figure 10 shows the
charge density differences between charge states for several of the transition
levels in Fig. 9. There are a few notable takeaways from these plots. First,
as expected, the charge density differences are larger for transitions that
result in large ionic displacements (as indicated by the MaxDisp in the bottom
of each plot). Second, these ionic displacements and resulting charge
distortions are noticeable more than 5 Å away from the defect site for some
transitions, indicating the importance of large supercells. Lastly, a slightly
larger fraction of the charge associated with the transition is contained
within a 5 Å sphere around the defect site for CIDER than for PBE (as
indicated by the Loc. metric in the bottom of each plot), suggesting that
CIDER localizes the defect charge more effectively than PBE. This is
reassuring because unphysical charge delocalization is one of the key driving
factors of erroneous charge transition level predictions in semiconductor
point defects [5].

In summary, our results indicate that the CIDER approach could be a route
toward high-accuracy defect calculations at similar cost to semilocal DFT,
which would significantly broaden the potential scope of point defects
research. In particular, it could be used to improve the accuracy of high-
throughput point defects studies [102], which are currently only
computationally feasible with semilocal functionals like PBE.

##  IV Conclusion

We have demonstrated that by combining the kinetic energy density (a
semilocal, orbital-dependent quantity) with nonlocal, density-dependent
features, it is possible to design functionals with comparable accuracy to
hybrid DFT at a drastically reduced computational cost for plane-wave DFT
calculations. We can do this by learning exact exchange energies, without the
use of any high-accuracy reference data from experiments or wave function
theory. We have shown that our nonlocal meta-GGA model has higher accuracy
than semilocal functionals on both molecular benchmarks and solid-state band
gap prediction, and that these improvements can assist with applications such
as semiconductor point defect calculations.

The developments presented here serve a variety of purposes. As is, the
exchange functional we designed can be used as an exact exchange surrogate in
scenarios where hybrid DFT is impractical, such as large extended systems.
This expands the scope of calculations that can be performed at hybrid DFT
accuracy. In addition, the techniques we have introduced to learn smooth,
physically constrained functionals and implement efficient nonlocal features
could be used in other contexts, such as orbital-free DFT [67] and van der
Waals functional design [53, 64, 103].

Importantly, the ability to systematically vary the descriptor complexity in
our machine learning approach allows us to extract physically nontrivial
information about the nature of the exchange interactions. Specifically, we
find that both nonlocality and derivative discontinuity are independently
important for the structure of the exchange energy density functional. This
clarifies a key limitation of existing semilocal functionals.

Lastly, our approach can be generalized to learn the full XC functional,
providing a path to highly transferable and universal functionals capable of
describing challenging and heterogeneous chemistry. We anticipate that these
developments will contribute to significantly improving the cost-accuracy
trade-off of DFT, accelerating the pace and efficacy of materials and chemical
research by providing a stronger theoretical backing in various sub-fields.

##  Appendix A Computational Details

###  A.1 Overview

To explore the effectiveness of different features for learning the exchange
energy, four types of functional were trained, with the model and feature
vector details provided in Section II.3:

  * •

SL-GGA: A semilocal GGA model based on the reduced squared gradient $s^{2}$
(Eq. 9).

  * •

NL-GGA: A nonlocal model based on $s^{2}$ and three nonlocal features with
Eqs. 20 and 21 for the exponents.

  * •

SL-MGGA: A semilocal meta-GGA model based on $s^{2}$ and the iso-orbital
indicator $\alpha$ (Eq. 11).

  * •

NL-MGGA: A nonlocal model based on $s^{2}$, $\alpha$, and three nonlocal
features with Eqs. 22 and 23 for the exponents.

The purpose of training semilocal GGA and meta-GGA models is to test the
expressive power of the nonlocal descriptors. In our previous work (see Table
S1 of Ref. [35]), we found that a standard, semilocal meta-GGA could not
effectively learn the exchange energy density, but in this work, we determined
that a semilocal meta-GGA can learn the _total_ exchange energy much more
effectively than it can learn the exchange energy density.

The two primary datasets used in this work were the GMTKN55 [10] and SOL62
[69, 33] databases. The GMTKN55 database contains 1505 molecular reaction
energies based on 2462 single-point calculations, covering five categories:
small-molecule properties, larger molecule reactions and isomerization,
barrier heights, intermolecular noncovalent interactions, and intramolecular
noncovalent interactions. The SOL62 dataset consists of 42 non-metals and 20
metals as categorized by Zhang _et al._ [69]. We follow Trepte and Voss [33]
in excluding Pb and Th from Zhang _et al._ ’s initial set of 64 solids. Twelve
subsets of GMTKN55 and half of SOL62, along with a few other systems, were
used for training and validation, as detailed in Section A.4.

Throughout this work, we used pymatgen [104] and ase [105] for structure
manipulation and calculation setup, fireworks [106] for automated workflow
management, and Scikit-learn [107] for the Gaussian process models. The
computational details for the PySCF and GPAW calculations are provided in
Appendices A.2 and A.3, respectively. Appendices A.4 and A.5 describe the
training and validation set selection, and Appendix A.6 covers the selection
of control points for the sparse Gaussian processes. Appendix A.7 describes
the selection of noise hyperparameters and training/validation loss functions,
and Appendix A.8 describes the list of hyperparameters tested for each model
and the model selection procedure. Appendix A.9 describes how we handled
convergence issues for CIDER models on GMTKN55, and Appendix A.10 describes
the PySCF and GPAW settings that are specific to the CIDER functionals.
Appendix A.11 describes the band gap benchmark methods. Appendix A.12
describes the computational cost benchmarking of CIDER for plane-wave DFT
calculations. Finally, Appendix A.13 describes our methodology for computing
charged defect transition levels.

###  A.2 Molecular Calculations and Reference Data Generation

Molecular systems, including the nanoclusters discussed in Appendix A.4, were
computed using the PySCF code [55, 56]. All calculations were performed using
the def2-QZVPPD basis set (with effective core potentials for atoms larger
than Kr) [108, 109] for the atomic orbitals and the def2-universal-jkfit
auxiliary basis set [110] for the classical Coulomb energy. The seminumerical
exchange (SGX) module in PySCF, which is similar to the computationally
efficient chain-of-spheres exchange algorithm [111], was used to compute the
exchange energy and potential for hybrid functionals. For SGX, level 1 PySCF
grids were used with P-junction screening turned off. For a given sub-database
in GMTKN55, restricted DFT was used if all systems in that sub-database had
$S_{z}=0$; otherwise, unrestricted DFT was used. The exception was the PCl_3
transition state in the INV24 database, for which spin symmetry breaking
occurs for some functionals. This system was calculated with unrestricted DFT
as well. For SCF calculations, the nonlocal features and XC energy were
computed on level 3 PySCF grids.

For benchmarking comparison, the GMTKN55 database was evaluated with the above
settings using the PBE [23], revPBE [77], r2SCAN [73], PBE0 [71], B3LYP [27],
and PW6B95 [76] functionals, all with D4 dispersion corrections [72], as well
as with the B97M-V [78], $\omega$B97X-V [28], and $\omega$B97M-V [79]
functionals.

To collect the exact exchange energy training data, the ground-state orbitals
from the PBE calculations were used to compute the exact exchange energy
$E_{\text{x}}^{\text{exact}}$ via Eq. 2, and the feature vectors were computed
from the corresponding PBE densities. The features were computed on level 1
PySCF grids. For training set feature generation only (not the SCF
calculations), the integral in Eq. 18 was computed numerically with the
$\mathbf{r}_{2}$ coordinate evaluated on level 3 PySCF grids without pruning.
This dense grid is necessary for numerical stability and precision in the core
region.

###  A.3 Solid-State Calculations and Reference Data Generation

The electron densities and ground-state energies of the SOL62 dataset were
computed in GPAW [62, 63] using the PBE geometries obtained by Trepte and Voss
[33], with an energy cutoff of 1000 eV and the same k-point meshes as Trepte
and Voss. The $h$ parameter in GPAW was set to the minimum value required to
avoid aliasing of the density, i.e. 0.096 Å. Fermi-Dirac smearing with a width
of 0.01 eV was used for periodic systems. For isolated atoms, fixed
occupations were used unless numerical instabilities occurred as a result,
which was the case for six atoms (Pb-2, Pd-2, Pt-2, Sn-2, Ta-5, and V-5, where
the number is the ground-state magnetic moment $2S$). In these cases, Fermi-
Dirac smearing with a width of 0.002 eV was used. The default Davidson solver
was used for all solids except the alkaline earth metals, which were solved
with the conjugate gradient (CG) method. The CG method was used for all
isolated atomic systems. To compute cohesive energies for training, the
experimental ground state magnetic moments were used for the isolated atom
references. The convergence threshold was set to $10^{-6}$ eV per valence
electron for both isolated atoms and periodic systems. The density convergence
threshold was set to $10^{-4}$ per electron for periodic systems and $10^{-3}$
per electron for isolated atoms. The eigenstates convergence threshold was set
to $4\times 10^{-8}$ eV2/electron for periodic systems and $10^{-5}$
eV2/electron for isolated atoms. The cell size for the isolated atoms was
chosen to be about 11.7 Å across, with perturbations included to break
symmetry. Compared to a larger box size of 17.5 Å, the PBE atomic energies had
a root mean square deviation of 0.012 eV and maximum deviation of 0.080 eV
(for Zr). Because of the small average deviation, and the fact that the train
and test reference values were computed on the atom-in-a-box systems rather
than determined by experiment, the smaller boxes were used to lower the
computational cost and memory requirements of training and testing. Symmetry
constraints were also turned off for the isolated atom calculations.

Training the CIDER models to the solid-state data required computing the exact
exchange energy and the nonlocal features for each system. To obtain exact
exchange energies, the hybrid module of GPAW was used to compute the non-self-
consistent exact exchange energy (EXX) from the PBE orbitals obtained in the
above calculations. However, because the k-point meshes above are too dense
for evaluation of EXX to be computationally feasible, the periodic systems
were first recomputed with a coarser k-point mesh (for each system, the
smallest even, $\Gamma$-centered mesh with a linear k-point density of at
least 4.5 k-points per Å-1), and the EXX was computed on this mesh. All
nonlocal density features were evaluated on the density and orbitals obtained
from the denser k-point mesh. They were evaluated with the internal CIDER
settings $q_{\text{max}}=300$ and $\lambda=1.8$ (see Section A.10 for
details).

The PBE0 cohesive energies were computed using the same procedure above, with
the PBE0 energies being defined as

$E^{\text{PBE0}}=E^{\text{PBE,dense}}-E^{\text{PBE,coarse}}+E^{\text{PBE0,coarse}},$
(104)

with _dense_ denoting the calculation done with the dense k-point mesh and
_coarse_ denoting the coarser k-point mesh discussed above.

###  A.4 Construction of Training and Validation Sets

The total exchange energy of a molecule is highly basis-set and code-
dependent, and it is dominated by the core contributions for large atoms. In
practice, only relative exchange energies between systems is relevant for most
chemistry and materials science applications. Therefore, rather than training
the model to match total exchange energies, we trained to match the exchange
energy differences between systems. For example, the target and predicted
values for the atomization energy of water would be

$\displaystyle E_{\text{x}}^{\text{target}}$
$\displaystyle=E_{\text{x}}^{\text{ex}}[n_{\text{\cf{O}}}]+2E_{\text{x}}^{\text{ex}}[n_{\text{\cf{H}}}]-E_{\text{x}}^{\text{ex}}[n_{\text{\cf{H2O}}}]$
(105) $\displaystyle E_{\text{x}}^{\text{pred}}$
$\displaystyle=\hat{E}_{\text{x}}[n_{\text{\cf{O}}}]+2\hat{E}_{\text{x}}[n_{\text{\cf{H}}}]-\hat{E}_{\text{x}}[n_{\text{\cf{H2O}}}],$
(106)

where $E_{\text{x}}^{\text{ex}}[n]$ and $\hat{E}_{\text{x}}[n]$ are the exact
and model exchange energies for density distribution $n(\mathbf{r})$,
respectively.

The training and validation sets were constructed from four sources:

Molecules: The W4-11 (atomization energies), BH76 (barrier heights), BH76RC
(reaction energies for BH76 reactions), MB16-43 (“mindless benchmarking” of
reaction energies), S22 (noncovalently bound dimers), G21IP (ionization
potentials), PA26 (proton affinities), RG18 (noble gas dimers), ACONF (alkane
conformers), ALKBDE10 (alkali/alkaline earth diatomic bond energies),
HEAVYSB11 (dissociation energies for molecules with heavy atoms), and SIE4x4
(self-interaction problems) subsets from the GMTKN55 database were used for
training and validation. See the original GMTKN55 paper [10] and references
therein for more details on the datasets. Each set was randomly split into
two-thirds training and one-third validation, except for the BH76 and BH76RC
subsets. Because the data points in these subsets are based on the same
systems, and also because they include the forward and backward reactions for
the same chemical equation, a different splitting procedure was chosen to
maintain train-validation independence. First, the 34 unique transition states
in BH76 were identified and split half-and-half into the training and
validation sets. All 40 reactions with the training transition states were
used for training, and all 36 reactions with the validation transition states
were used for validation. All 17 BH76RC reactions with systems in the BH76
training partition were used for training, and the other 13 reactions were
used for validation. All data not in the training or validation set was
classified as the GMTKN55 test set. This data partitioning scheme resulted in
280 train data, 163 validation data, and 1062 test data in GMTKN55.

Bulk solids: Half of the SOL62 cohesive energy database (11 randomly selected
elemental metals and 20 randomly selected non-metals) was used for training
and validation. Nine systems (C-diamond, Li, GaAs, MgO, LiF, ZrC, TiC, Pd, and
Pt) were hand-selected for the training set to cover a range of chemical
space, while the other 22 systems were used for validation. All data not in
the training or validation set was classified as the SOL62 test set.

Nanoclusters: Six nanocluster atomization energies were added to the training
set as a “bridge” between the isolated and extended systems in the study.
These were Ca13 and Zn26, with the initial structure obtained based on the
Materials Project [112] bulk lattice constant and structure optimization at
the PBE/def2-TZVP level in Orca [113, 114, 115]; C35H36 (hydrogen-capped
diamond), obtained from Refs. [116, 117] and also structurally optimized in
Orca at the PBE/def2-TZVP level; and Na4Cl4, [Na13Cl14]-, and [Na14Cl13]+, all
using the Materials Project bulk lattice constant without further structure
optimization.

Atoms: Absolute and relative energies of the isolated atoms H-Kr and a few
excited spin states were included in the training set. Specifically, the
training points consisted of the absolute energies of the H atom and the noble
gases He, Ne, Ar, and Kr; the relative energies of the elements with nuclear
charge $Z$ and $Z+1$ for $Z=\\{1,...,35\\}$, and the relative energies of the
following ground state/excited state spin pairs: Sc 1/3, Ti 2/4, V 3/5, Cr
6/4, and Ni 2/0. The preceding numbers correspond to the magnetic moment $2S$.

In addition to fitting to the above training data, the uniform electron gas
constraint was enforced as described in Section II.3.3, and the uniform
scaling constraint was enforced by the scale-invariance of the input features.
The ground-state total energies of the SOL62 solids and isolated atoms were
computed with GPAW as described in Section A.3. The ground-state energies for
all other training data were evaluated with PySCF as described in Section A.2.
The selection of train and validation databases from GMTKN55 includes a wide
variety of chemistries and properties. However, because training data is
sampled disproportionately from small-molecule properties, there is a large
distribution shift between the train/validation and test sets, with the test
sets containing more large molecules, isomerization reactions, and barrier
heights. The use of this train/test partition therefore serves as a strong
test of transferability. The accuracy on SOL62 also strongly reflects the
transferability of a model because only 9 solid-state systems were included in
the training set, and SOL62 covers a wide variety of chemistries.

###  A.5 A More Diverse Training Set for Numerical Stability

To design a more diverse training set for the NL-MGGA-DTR functional, a kernel
function was defined for the “covariance” between two training points in
GMTKN55:

$k_{\text{GMTKN}}(\mathbf{x},\mathbf{x}^{\prime})=\exp\left(-\frac{\mathbf{x}\cdot\mathbf{x}^{\prime}}{328}\right),$
(107)

where $\mathbf{x}$ is the vector of errors for 82 functionals benchmarked on
the GMTKN55 database by Goerigk _et al._ [10] with D3(BJ) dispersion
corrections [118, 119, 120], normalized so that the mean square error across
all functionals for that individual data point is 1. This kernel essentially
provides a larger similarity score between data points for which the
(normalized) errors tend to be similar for a given functional. The kernel
matrix was constructed for the entire GMTKN55 database and then factorized
using a method called Cholesky decomposition with complete pivoting [121,
122], which factors a given matrix $\mathbf{A}$ as follows:

$\mathbf{P}^{\top}\mathbf{A}\mathbf{P}=\mathbf{L}\mathbf{L}^{\top},$ (108)

In Eq. 108, $\mathbf{P}$ is a permutation matrix that orders the vectors
$\mathbf{x}$ in Equation 107 in decreasing order of their variance,
conditioned on the preceding feature vectors, _in the space implicitly spanned
by the kernel function_. Intuitively, this means that later vectors, as
ordered by $\mathbf{P}$, are likely to be nearly linearly dependent (in the
kernel space) on the earlier vectors. Therefore, by picking the first $N$
training points as ordered by $\mathbf{P}$ and using the remaining data for
testing, we can get a training set whose errors correlate as much as possible
with all the test data, effectively making the training set more diverse and
potentially preventing out-of-distribution error of the ML model.

We follow this procedure and select the first 225 pivots for the training set,
the next 75 pivots for the validation set, and the remaining 1205 pivots for
the test set. With this dataset partitioning, more than 99% of the GMTKN55
test set data points contain at least one system not present in the train or
validation sets, so the test set is still sufficiently independent to assess
the transferability of the functionals to new systems. All training points not
contained in GMTKN55 were kept the same as in Appendix A.4.

###  A.6 Control Point Selection for Sparse Gaussian processes

In Eqs. 35 and 36, a set of control points $\\{\mathbf{\tilde{x}}_{a}\\}$ is
required to evaluate the Gaussian process covariance kernel in the sparse
approximation. These control points were initially selected by randomly
sampling feature vectors from the PySCF grids for the training set (Appendix
A.2). All training data was sampled except the SOL62 data, since the features
for these systems were evaluated in GPAW.

Due to the size of the numerical grids from which the features were sampled,
there were a large number of potential control points for the sparse Gaussian
process, many of which were nearly linearly dependent in the implicit feature
space of the covariance kernel function. To reduce the size of the control
point set, a Cholesky decomposition with pivoting (Eq. 108) was performed on
an initial sample size of roughly $10^{4}$ feature vector samples with a
tolerance of $10^{-5}$, and only control points corresponding to signficant
pivot indices were kept. This reduced the number of control points from
roughly $10^{4}$ to $10^{2}$.

###  A.7 Selection of Noise Hyperparameters, Training Loss, and Validation
Loss

Chemical and material properties have a wide range of magnitudes, from less
than 1 kcal/mol for some van der Waals binding energies to greater than 100
kcal/mol for many covalent bond dissociation energies. Correspondingly,
different magnitudes of XC errors are expected and acceptable for these
properties; an error of 1 kcal/mol might make a van der Waals binding energy
prediction useless but be acceptable for an atomization energy. To train the
exchange functionals, it was therefore necessary to weight the noise
hyperparameter of each training point to reflect the feasible and desired
precision for the property of interest. To do so, we used the data provided
with GMTKN55 [10] to compute the MAD of four hybrid functionals (PBE0-D3(BJ),
B3LYP-D3(BJ), PW6B95-D3(BJ), and $\omega$B97X-V) compared to the reference
data on each sub-database (denoted $\text{MAD}^{\text{hybrid}}_{\text{DB}}$,
where DB is the sub-database of interest), and then created characteristic
uncertainties for each sub-database as follows (in Hartree atomic units):

$\tilde{\sigma}_{\text{DB}}=\text{MAD}^{\text{hybrid}}_{\text{DB}}\left(\frac{\sigma_{0}}{\text{MAD}^{\text{hybrid}}_{\text{W4-11}}}\right),$
(109)

where W4-11 is a small-molecule atomization energy sub-database in GMTKN55 and
$\sigma_{0}=0.03$ Eh. One drawback of the above approach is that it favors
high accuracy on systems that are already described well by conventional,
dispersion-corrected DFT (like simple van der Waals dimers) and assumes a
large uncertainty on systems like self-interaction-dominated bond energies,
for which even hybrid functionals with small exchange mixing perform poorly.
To remedy this, the noise parameter was augmented with a constant weighting
for every system:

$\sigma_{\text{DB}}=\sqrt{\frac{2}{1/\sigma_{0}^{2}+1/\tilde{\sigma}_{\text{DB}}^{2}}}.$
(110)

With this definition, the matrix $\Sigma_{\text{noise}}$ in Eq. 41 is a
diagonal matrix with $\left(\Sigma_{\text{noise}}\right)_{ii}$ being the
$\sigma_{\text{DB}}^{2}$ for the database to which training point $i$ belongs.
Constructing the Gaussian process for this noise parameter is equivalent to
regression with the loss function

$\mathcal{L}=\frac{1}{2}\sum_{\text{DB}}N_{\text{DB}}\text{MSD}_{\text{DB}}\left(\frac{1}{\sigma_{0}^{2}}+\frac{1}{\tilde{\sigma}_{\text{DB}}^{2}}\right),$
(111)

where $\text{MSD}_{\text{DB}}$ stands for the mean squared deviation of a
given GMTKN55 sub-database labeled by “DB,” and $N_{\text{DB}}$ is the number
of points in the sub-database. This weighting ensures that training points
with small expected errors are appropriately weighted while those with larger
expected errors for conventional functionals do not become under-weighted. For
systems outside the GMTKN55 database that were in the training set,
$\tilde{\sigma}_{\text{DB}}$ was set by hand.

###  A.8 Hyperparameter Selection and Validation Procedure

The models tested in this work have several discrete and continuous
hyperparameters to be optimized, which are listed below.

  * •

The covariance prefactor hyperparameter $\Sigma_{\text{Cov}}$ in eqs 54–57.

  * •

The kernel length-scale ($l_{i}$, where $i$ is the feature index) in eq 58.

  * •

For the nonlocal functionals only, the feature exponent parameters $B_{i}$ and
$C_{i}$ in Eqs. 20–23.

  * •

The baseline exchange functional on top of which the CIDER model is trained
via $\Delta$-learning. We tested PBE [23] and the Chachiyo GGA [70].

It would be quite challenging and expensive to implement differentiation of
the marginal likelihood with respect to the second parameter, and nearly
impossible for the third since backpropagation would need to be performed
through the feature evaluation step. Therefore, we used various combinations
of hyperparameters to train the models, and then for each functional type, we
selected the best-performing model on the validation set as the final
candidate functional.

To obtain a discrete set of exponent hyperparameters to test for the nonlocal
functionals, we developed a reduced number of parameters to specify the
feature exponents. Similarly to the previous CIDER work [35], a control
parameter $A$ tunes $B_{2}$ and $C_{2}$ simultaneously, in two different
schemes. The first scheme (S1) is

$\displaystyle B_{2}$ $\displaystyle=A$ $\displaystyle C_{2}$
$\displaystyle=\frac{A}{32}\frac{6}{5\pi}(6\pi^{2})^{2/3}\approx 0.18A.$ (112)

The second scheme (S2) provides a different mixture of density-based and
gradient/kinetic energy-based terms:

$\displaystyle B_{2}$
$\displaystyle=\frac{A}{16}\frac{6}{5\pi}(6\pi^{2})^{2/3}\approx 0.36A$
$\displaystyle C_{2}$ $\displaystyle=B_{2}.$ (113)

In both schemes, the other constants are determined by $B_{2}$ and $C_{2}$:

$\displaystyle
B_{1}=\frac{1}{2}B_{2},\,\,B_{3}=2B_{2},\,\,B_{0}=\frac{D}{A}B_{2}$ (114)
$\displaystyle
C_{1}=\frac{1}{2}C_{2},\,\,C_{3}=2C_{2},\,\,C_{0}=\frac{D}{A}C_{2},$ (115)

where $D$ is an additional parameter to tune the $B_{0}$ and $C_{0}$
coefficients separately from the other coefficients. With these schemes, all 8
parameters $\\{B_{i}\\},\\{C_{i}\\}$ that define the three-feature CIDER model
were reduced to two parameters per scheme. For each scheme, we performed a
hyperparameter search over the grid $A=\\{1/2,1,2,4,8\\}$ and
$D=\\{1/2,1,2,4,8\\}$, resulting in a total of 50 feature constructions to
test.

Because there were already 50 features from which to select the nonlocal
functionals, we set the covariance and feature length-scales for these models
heuristically:

$\displaystyle\Sigma_{\text{Cov}}$
$\displaystyle=R_{1}\expectationvalue{\Delta F_{\text{x}}^{2}}_{\text{train}}$
(116) $\displaystyle l_{i}$
$\displaystyle=R_{2}\sqrt{\expectationvalue{x_{i}^{2}}_{\text{train}}},$ (117)

where $\expectationvalue{}_{\text{train}}$ is an expectation value over all
the grid points in the exchange energy training set, $\Delta F_{\text{x}}$ is
the deviation of the exact exchange enhancement factor at a given point from
the baseline functional enhancement factor, $x_{i}$ is the $i$-th regularized
feature (see Section II.3.3), and $R_{1}$ and $R_{2}$ are constants. $R_{1}$
and $R_{2}$ were set based on the criteria that the resulting functionals
converge consistently in both GTO and plane-wave calculations and exhibit a
small total energy deviation between the PAW and all-electron GTO atomization
energies for a few small molecules. $R_{2}=1$ worked best for both the
nonlocal GGA and meta-GGA, while the most practical values of $R_{1}$ depended
on the baseline functional and model type, as shown in Table 6. The meta-GGA
models required a smaller covariance prefactor to achieve good stability than
the GGA models. We set $R_{1}$ to be larger by a factor of 20 for functionals
with the Chachiyo baseline compared to the PBE baseline because
$\expectationvalue{\Delta F_{\text{x}}^{2}}_{\text{train}}$ was roughly 20
times larger for PBE than Chachiyo.

Table 6: $R_{1}$ values for each functional type and baseline functional. Type/Baseline | $R_{1}$
---|---
NL-GGA/PBE | $1.0$
NL-GGA/Chachiyo | $20.0$
NL-MGGA/PBE | $0.05$
NL-MGGA/Chachiyo | $1.0$

For the semilocal functionals, there was no need to validate over the possible
feature constructions, leaving more flexibility to test $R_{1}$ and $R_{2}$
systematically. For the functionals with Chachiyo baseline, we tested
$R_{1}=\left\\{1,2,4,8,16,32,64,128\right\\}$, while for the PBE baseline each
$R_{1}$ was divided by 20 to give a roughly equivalent $\Sigma_{\text{Cov}}$.
For the semilocal functionals,
$R_{2}=\left\\{1,\frac{1}{2},\frac{1}{4}\right\\}$ were tested, resulting in
24 functionals to validate. Note that for NL-MGGA-DTR, we determined from
validation set SCF calculations that $R_{1}=0.1$ with the PBE baseline
resulted in sufficiently stable models, so this larger value was used. In
addition, we did not search over feature length-scale for NL-MGGA-DTR, and
instead used the features determined to be most accurate for NL-MGGA and NL-
MGGA-PBE.

After the functionals were trained, we selected the top 10 functionals of each
type and baseline functional, as classified by non-self-consistent performance
on the validation set via Eq. 111. To validate these functionals, the GMTKN55
train and validation set reaction energies were predicted using self-
consistent calculations with a convergence threshold of $10^{-6}$ Eh and a
gradient threshold of $2\times 10^{-3}$ Eh. For a given functional, if any
calculation failed to converge when starting from both the default initial
guess of PySCF and from the PBE orbitals, the functional was deemed
insufficiently numerically stable and disqualified from further consideration.
The non-self-consistent cohesive energies of the SOL62 validation set were
also computed. Of the functionals for which all calculations converged, the
validation loss consisted of

$\displaystyle\mathcal{L}_{\text{val}}$
$\displaystyle=\mathcal{L}_{\text{MOL}}+\mathcal{L}_{\text{SOL}}$ (118)
$\displaystyle\mathcal{L}_{\text{MOL}}$
$\displaystyle=\frac{1}{2}\sum_{\text{DB}}N_{\text{DB/VAL}}\times\text{MSD}_{\text{DB/VAL}}^{\text{SCF}}\times\left(1+\frac{\sigma_{0}^{2}}{\tilde{\sigma}_{\text{DB}}^{2}}\right)$
(119) $\displaystyle\mathcal{L}_{\text{SOL}}$
$\displaystyle=\text{MSD}_{\text{SOL62/VAL22}}^{\text{NSCF}},$ (120)

where $\text{MSD}_{\text{DB/VAL}}^{\text{SCF}}$ is the self-consistent mean
square deviation on the validation partition of GMTKN55 subset DB,
$N_{\text{DB/VAL}}$ is the number of validation data in DB, and
$\text{MSD}_{\text{SOL62/VAL22}}^{\text{NSCF}}$ is the non-self-consistent
mean square deviation of the 22 validation-set cohesive energies per atom. The
final functional of each type was selected based on which candidate had the
lowest loss for Eq. 118. The validation results and selected hyperparameters
for the functionals are provided in the Supplemental Material (Section S4).

###  A.9 Difficult Convergence Cases

The default convergence tolerance for the calculations performed with CIDER
functionals was $10^{-8}$ Eh, except for the NL-MGGA-PBE and NL-MGGA-DTR
calculations, which were converged to $10^{-9}$ Eh when possible. Some CIDER