# Exact constraints and appropriate norms in machine learned exchange-
correlation functionals

Kanun Pokharel<EMAIL_ADDRESS>Department of Physics and Engineering
Physics, Tulane University, Louisiana 70118 New Orleans, USA  James W. Furness
<EMAIL_ADDRESS>Department of Physics and Engineering Physics,
Tulane University, Louisiana 70118 New Orleans, USA  Yi Yao  Thomas Lord
Department of Mechanical Engineering and Material Science, Duke University,
Durham, North Carolina 27708, USA  University of North Carolina at Chapel
Hill, Chapel Hill, NC 27599, USA  Volker Blum  Thomas Lord Department of
Mechanical Engineering and Material Science, Duke University, Durham, North
Carolina 27708, USA  Tom J.P. Irons  School of Chemistry, University of
Nottingham, Nottingham NG7 2RD, UK  Andrew M. Teale  Department of Chemistry,
University of Oslo, Nottingham NG7 2RD, UK  Hylleraas Centre for Quantum
Molecular Sciences, Department of Chemistry, University of Oslo, P.O. Box
1033, Blindern, N-0315 Oslo, Norway  Jianwei Sun<EMAIL_ADDRESS>Department
of Physics and Engineering Physics, Tulane University, Louisiana 70118 New
Orleans, USA

###### Abstract

Machine learning techniques have received growing attention as an alternative
strategy for developing general-purpose density functional approximations,
augmenting the historically successful approach of human designed functionals
derived to obey mathematical constraints known for the exact exchange-
correlation functional. More recently efforts have been made to reconcile the
two techniques, integrating machine learning and exact-constraint
satisfaction. We continue this integrated approach, designing a deep neural
network that exploits the exact constraint and appropriate norm philosophy to
deorbitalize the strongly constrained and appropriately normed SCAN
functional. The deep neural network is trained to replicate the SCAN
functional from only electron density and local derivative information,
avoiding use of the orbital dependent kinetic energy density. The performance
and transferability of the machine learned functional are demonstrated for
molecular and periodic systems.

##  I Introduction

The density functional theory (DFT) of Hohenberg, Kohn, and Sham [1; 2] allows
for efficient computation of material properties by avoiding the complicated
many-electron wave function in favor of the computationally convenient
electron density when solving the electronic structure problem. Due to its
useful accuracy and efficiency, DFT has become the most widely used
computational approach for solving problems in chemistry and condensed matter
physics/electronic structure.

In the Kohn–Sham (KS) formulation of DFT, the majority of the energy is
calculated exactly, leaving only a small portion of the energy, known as the
exchange-correlation (XC) energy, to be approximated. There has been extensive
research on improving approximations to the XC energy and the resulting
functionals are roughly categorized into a hierarchy of increasing complexity
and expected accuracy [3]. The meta-generalized gradient approximations
(mGGAs), the highest category that depends only on semi-local ingredients, are
becoming increasingly popular for allowing high accuracy at favorable
computational cost. The mGGA functionals commonly consist of a per-particle XC
energy density, $\epsilon_{\mathrm{xc}}$, built from three ingredients: the
electron density, $n(\mathbf{r})$, its gradient, |$\nabla{n(\mathbf{r})}$|,
and the kinetic energy density
$\tau(\mathbf{r})=\sum_{i}^{\mathrm{occ.}}|\nabla\varphi_{i}(\mathbf{r})|^{2}$,
where $\varphi_{i}$ are the occupied KS orbitals. Though less common, density
functionals depending on the density Laplacian $\nabla^{2}{n(\mathbf{r})}$
instead of (or in addition to) $\tau$ are also included at the mGGA level. The
total XC energy for the system is calculated by multiplying this per-particle
XC energy density by the local electron density $n(\mathbf{r})$ and
integrating over all space,

$\displaystyle E_{\mathrm{xc}}$ $\displaystyle=\int
d\mathbf{r}\epsilon^{\mathrm{mGGA}}_{\mathrm{xc}}(n(\mathbf{r}),|\nabla{n(\mathbf{r})}|,\tau(\mathbf{r}),\nabla^{2}{n(\mathbf{r})})n(\mathbf{r}),$
(1) $\displaystyle=\int
d\mathbf{r}F_{\mathrm{xc}}^{\mathrm{mGGA}}(n(\mathbf{r}),|\nabla{n(\mathbf{r})}|,\tau(\mathbf{r}),\nabla^{2}{n(\mathbf{r})})n(\mathbf{r})\epsilon_{\mathrm{x}}^{\mathrm{LDA}}(\mathbf{r}),$
(2)

where $F_{\mathrm{xc}}^{\mathrm{mGGA}}$ is the XC enhancement factor and
$\epsilon_{\mathrm{x}}^{\mathrm{LDA}}(\mathbf{r})=-(3/4\pi^{2})(3\pi^{2}n)^{1/3}$
is the exchange energy per-particle of the uniform electron gas.

The kinetic energy density, $\tau$, is commonly used in mGGAs to recognise
different chemical environments through iso-orbital indicator variables [4; 5;
6], and as a component of the spherical exchange hole expansion[7]. While
theoretically convenient, $\tau$ introduces an implicit dependence on the KS
orbitals, which brings some complications. 1) It reduces computational
efficiency by requiring additional basis function derivatives to be computed
on the numerical quadrature grid, which can be more costly for Fourier
transform based periodic codes. 2) It prevents the functional being used in
orbital-free DFT calculations. 3) Evaluation of the XC potential for
$\tau$-dependent functionals requires either optimised effective potential
(OEP) techniques [8; 9], or a generalised KS scheme[10; 11; 12]. While a
generalised KS treatment can be computationally convenient, the effective XC
potential operator of a $\tau$ dependent mGGA is no longer a multiplicative
function, $v_{\mathrm{xc}}(\mathbf{r})$, and is instead a non-local operator,
$\hat{v}_{\mathrm{xc}}$.

Despite the potential advantages offered avoiding the use of orbital dependent
ingredients such as $\tau$, $\nabla^{2}{n(\mathbf{r})}$ remains a less
explored ingredient and its physical significance for the XC energy is
unclear. Recently, Mejia-Rodriguez and Trickey [13; 14] replaced
$\tau(\mathbf{r})$ with functions of $\nabla^{2}{n(\mathbf{r})}$ in many mGGA
XC functionals to recover similar (but not identical) performance to the
parent functionals. This suggests an intriguing but unclear relationship
between $\tau(\mathbf{r})$ and $\nabla^{2}{n(\mathbf{r})}$, though an explicit
relationship remains elusive despite significant effort [15].

Machine learning (ML) has proven to be a powerful tool for building
complicated non-linear mappings for which little theoretical guidance exists.
It has proved successful in building complex models across a wide variety of
fields including robotics [16; 17], pattern recognition[18; 19], drug design
[20; 21; 22], and gaming [23]. Within DFT research, there has been a recent
practice of applying ML to construct density functionals. In 2012, Synder et.
al. used a ML approximation to construct an orbital-free non-interacting
kinetic energy functional $T_{\mathrm{s}}[n]$ for spinless fermion systems
[24; 25]. Brockherde et. al. used ML to learn the Hohenberg-Kohn (HK) map
between electron density and external potential to give a mechanism that
bypasses solving the KS equations [26]. Several other works have focused on
the XC potential problem [27; 28; 29; 30]. A perspective surveying the current
state of ML in computational chemistry and materials science was recently
published by Westermayr, Gastegger, and Schütt in Ref. [31].

The Strongly Constrained and Appropriately Normed (SCAN) functional [5] has
proven to be effective for describing a wide variety of systems [32], such as
liquid water and ice [33], semiconductor materials [34] and metal oxides [35],
as well as for key properties of correlated materials like cuprates [36; 37;
38]. SCAN’s success is credited to its adherence to all of the known exact-
constraints applicable to a meta-GGA functional, along with the philosophy of
using “appropriate norms” to set free parameters with minimal empiricism.
These appropriate norms are provided by the systems for which a semi-local
density functional approximation can be expected to be highly accurate, that
is: the total energies of systems with highly localized exchange-correlation
holes [5].

In light of this dual success of both constraint driven design and ML
techniques, a question arises. Is the philosophy of exact constraints and
appropriate norms compatible with ML for functional design? Growing evidence
that exact constraint adherence can improve ML transferability suggests it is.
Indeed, earlier works from Hollingsworth et. al. show that enforcing
coordinate scaling constraints can improve machine-learned functionals [39].
Nagai and collaborators [40] recently introduced a method to analytically
impose asymptotic constraints on an ML XC functional, finding generally
improved accuracy. More recently, Kirkpatrick and co-workers [41] developed a
functional DM21 (DeepMind21) that for the first time, obeys two classes of
constraints on systems with fractional electrons, which are fractional charge
systems and fractional spin systems[41].

In our work, we explore this idea of exact constraints and appropriate norms
satisfaction by training a deep artificial neural network (ANN) to reproduce
the XC energy density of the SCAN functional using $\nabla^{2}{n(\mathbf{r})}$
instead of $\tau(\mathbf{r})$, a similar goal to the SCAN-L functional [14].
The de-orbitalization of SCAN stands as a convenient application for exploring
the idea of constraint satisfaction in ML functionals, with SCAN-L providing
an analytical benchmark for the task.

The ML models will be trained to perform the transformation,

$\displaystyle
F_{\mathrm{xc}}^{\mathrm{SCAN}}\left(n(\mathbf{r}),|\nabla{n(\mathbf{r})}|,\tau(\mathbf{r})\right)$
$\displaystyle\Big{\downarrow}$ (3) $\displaystyle
F_{\mathrm{xc}}^{\mathrm{SCAN-
ML}}\left(n(\mathbf{r}),|\nabla{n(\mathbf{r})}|,\nabla^{2}{n(\mathbf{r})}\right)$

We approach this mapping of orbital free ingredients onto the SCAN XC energy
density using two different ML models adhering to different numbers of exact
conditions. One model is a single totally-connected ANN trained for Eq. 3
directly, termed the “combined model”. The other model is built as two
complementary exchange and correlation-like ANNs designed to obey exact spin-
scaling constraints, termed the “spin-scaled model”. We also impose the
general Lieb–Oxford bound [42] on these models to create models that also
satisfy this constraint.

## Exact Constraints

While the exact XC functional remains unknown, it is known to obey many
mathematical conditions, commonly called the “exact constraints” of XC
functionals. Currently, 17 exact constraints [5], are known to apply at the
semi-local functional level. These can be broken down as conditions for the
exchange energy: (1) negativity, (2) spin-scaling [43], (3) uniform density
scaling [44], (4) the slowly-varying density gradient expansion (to fourth
order) [45], (5) non-uniform density scaling [46; 47], and (6) a tight bound
for two-electron densities [48; 42]. For correlation: (7) non-positivity, (8)
the slowly-varying density gradient expansion (to second order) [49], (9)
uniform density scaling to the high-density limit [44], (10) uniform density
scaling to the low-density limit [44], (11) zero correlation energy for any
one-electron spin-polarized density, and (12) nonuniform density scaling [46;
47]. Finally, there are constraints known for the exchange and correlation
together: (13) size extensivity, (14) the general Lieb–Oxford bound [50; 51;
42], (15) weak dependence upon relative spin polarization in the low-density
limit [52; 53], (16) static linear response of the uniform electron gas [54],
and (17) the Lieb–Oxford bound for two-electron densities [42].

Here, we consider a subset of the 17 constraints that are easy to enforce in
an ML model. The first of these is the behavior of the exchange energy under
constraint (3): uniform density scaling,

$n_{\gamma}(\mathbf{r})=\gamma^{3}n(\gamma\mathbf{r}),$ (4)

where $\gamma$ is a positive real number. The exact exchange energy is known
to scale as,

$E_{\mathrm{x}}[n_{\gamma}]=\gamma E_{\mathrm{x}}[n],$ (5)

under this transformation. The effect of this condition on kernel ridge
regression models was investigated by Hollingsworth, Baker, and Burke in Ref.
39 for Hooke’s atom model systems, concluding that its inclusion improved ML
functional performance.

Constraint (2), the spin-scaling relation for exchange energy,

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

is simple to enforce for ML exchange models by requiring separate exchange and
correlation models, and that the same exchange model handle each spin channel
independently.

Constraint (14), the general Lieb–Oxford bound on the XC enhancement factor,
states that,

$0\leq F_{\mathrm{xc}}(\mathbf{r})\leq 2.215.$ (7)

These bounds can be enforced on ML models by including a post-processing step
that maps the ML model output, denoted
$\mathrm{ANN}_{\mathrm{xc}}(\mathbf{r})$, to the desired domain, e.g.

$F_{\mathrm{xc}}^{\mathrm{ML}}(\mathbf{r})=\frac{2.215}{1+\mathrm{ANN}_{\mathrm{xc}}(\mathbf{r})^{2}}.$
(8)

A similar approach can be applied to impose constraint (6), the tight bound
for the exchange enhancement factor $F_{\mathrm{x}}(\mathbf{r})\in[0,1.174]$,
if the exchange and correlation models are separated. Conveniently, such post-
processing also enforces constraints (1) and (7), non-positivity, by
constraining $F_{\mathrm{xc}}\geq 0$. It appears more challenging to enforce
the exact constraints outside this subset in ML models. For example, while
enforcing the second (and fourth) order gradient expansions for correlation
(and exchange) is relatively straightforwards in analytical functionals, the
ML design contains thousands of parameters which cannot be fully controlled.
Thus, it is non-trivial to enforce such gradient expansion constraints on the
model _a priori_. Despite this, the nature of supervised training against
methods that obey such constraints (such as SCAN) will result in the trained
model effectively learning the constraints to some degree. However, without
the rigorous enforcement described above it is unclear how well such adherence
will transfer out of the training domain.

## Input domain

Identifying the input domain is a critical part of ML model design as a
model’s performance can be strongly dependent upon the nature of its inputs.
Since our central interest is to replace the kinetic energy density
$\tau(\mathbf{r})$ dependence we will only consider orbital free ingredients.
Four density inputs were initially identified,

$\displaystyle r_{\mathrm{s}}=\left(\frac{3}{4{\pi}n}\right)^{1/3}$
$\displaystyle-\text{Wigner-Seitz radius}$ (9) $\displaystyle s=\frac{|\nabla
n|}{2(3\pi^{2})^{1/3}n^{4/3}}$ $\displaystyle-\text{Reduced density gradient}$
(10)
$\displaystyle\zeta=\frac{n_{\uparrow}-n_{\downarrow}}{n_{\uparrow}+n_{\downarrow}}$
$\displaystyle-\text{Spin polarization}$ (11) $\displaystyle
q=\frac{\nabla^{2}n}{4(3\pi^{2})^{2/3}n^{5/3}}$ $\displaystyle-\text{Reduced
density laplacian}$ (12)

The Weigner–Seitz radius is the radius of a sphere which on average contains
for a uniform density n. The reduced density gradient introduces inhomogeneity
which measures how fast and how much the density varies on the scale of the
local Fermi wavelength $2\pi/k_{\mathrm{F}}$ where
$k_{\mathrm{F}}=(3\pi^{2}n)^{1/3}$ is the Fermi wavevector. The reduced
density Laplacian also measures the density inhomogeneity and can distinguish
bonds, in contrast to reduced gradient which vanishes at the middle of the
bond [15]. The above dimensionless ingredients are preferred for XC
functionals rather than using the density variables directly as the correct
uniform coordinate density-scaling behaviour can be satisfied with them in
conventional XC functionals [44]. We note however that the non-linear
manipulations of these variables made by the ML models likely break the formal
scaling limits in practice. Particularly, we supply density information to the
exchange model through $r_{s}$, thereby violating the exchange coordinate-
scaling condition.

Further exploration revealed that including additional ingredients directly
from the SCAN exchange and correlation functionals (see supplemental material
of Ref. 5) could improve model performance:

$\displaystyle\epsilon_{\mathrm{c}}^{0}(r_{\mathrm{s}},s,\zeta)$
$\displaystyle-\text{Single orbital correlation}$ (13)
$\displaystyle\epsilon_{\mathrm{c}}^{1}(r_{\mathrm{s}},s,\zeta)$
$\displaystyle-\text{Slowly varying correlation}$ (14) $\displaystyle
g_{\mathrm{x}}(s)$ $\displaystyle-\text{Exchange inhomogeneity}$ (15)
$\displaystyle h_{\mathrm{0}}=1.174$ $\displaystyle-\text{Single orbital
exchange}$ (16) $\displaystyle h_{\mathrm{1}}={1.065}$
$\displaystyle-\frac{0.065}{\left(1+\frac{10s^{2}/81}{0.065}\right)}-\text{X
$2^{\mathrm{nd}}$ order gradient expansion}$ (17)

These additional inputs are combinations of the original density ingredients
and do not provide any new information directly. Their inclusion makes
learning more efficient however, as it reduces the manipulations that the
network must learn. This limits dependence on the network to only learning
$\tau$ dependent aspects, rather than requiring it to learn every detail of
the SCAN functional.

The possible range of the input parameters is very different from the desired
range of the model outputs: $0\leq F_{\mathrm{xc}}^{\mathrm{ML}}\leq 2.215$.
For example, the domain of $r_{s}$ and $s$ is $[0,+\infty)$ while $q$ is
$(-\infty,\infty)$. Such a mismatch in the magnitude of input and output is
known to be challenging for ML models. To correct for this we pre-processed
the unbounded inputs using the hyperbolic tangent function, $\tanh(x)$ [55],
to smoothly map the unbounded quantities to $(-1,+1)$. With pre-processing,
the inputs are defined as,

$\displaystyle\tilde{r}_{\mathrm{s}}$
$\displaystyle=\mathrm{tanh}(r_{\mathrm{s}}),$ (18) $\displaystyle\tilde{s}$
$\displaystyle=\mathrm{tanh}(s),$ (19) $\displaystyle\tilde{q}$
$\displaystyle=\mathrm{tanh}(q).$ (20)

Finally, we pre-process $\zeta$ as,

$\tilde{\zeta}=\frac{1}{2}\left[(1+\zeta)^{4/3}+(1-\zeta)^{4/3}\right],$ (21)

to ensure the model is a symmetric function of spin polarisation [56]. The
inputs of Eqs. 13-17 are unprocessed as their ranges are already correctly
bounded. Note that the pre-processed variables (Eqs. 18-21) are only supplied
to the network and are not used to generate the additional inputs of Eqs.
13-17.

Having identified the input domain, a training data set consisting of 20 atoms
was generated using accurate spherical Hartree–Fock Slater-type orbitals [57;
58; 59] containing open and closed shell atoms (He, Li, Be, B, C, N, O, F, Ne,
Na, P, Cl, Ar, K, Cr, Cu, Cu+, As, Kr, Xe) with $s$, $p$, and $d$ valence
shells. Spherical Hartree–Fock Slater-type orbitals were chosen to match the
data used to fit parameters in SCAN, though we expand the set of atoms
significantly beyond the rare gas atoms used to fit SCAN. The energetically
important region of the atomic density is typically between 0 and 4 $a_{0}$
and the density of each atom was therefore sampled in shells of decreasing
sample density for models, with 2500 radial points uniformly sampling in $r<1$
bohr (core), 1500 in $1\leq r<4$ bohr, and 500 in the tail region $4\leq r<10$
bohr.

This atomic training data was augmented with densities from the “appropriate
norm” systems used in SCAN’s construction. The first norm is the one electron
hydrogen atom, which is used to ensure that SCAN is one-electron self
interaction free (constraint 11). The second and third norms are the jellium
surface densities for ${r}_{\mathrm{s}}$ = 2, 3, 4 and 6 [60; 61], and the
converged SCAN orbitals of the compressed argon dimer with nuclear separations
of 1.6, 1.8, and 2.0 Å, which were used to fix SCAN’s interpolation function
parameters [5]. We restrict the training data set to only these appropriate
norms and avoid including chemically bonded systems. While increasing the
domain of the training data is expected to give the resulting ML models higher
accuracy for a wider range of problems, our intent here is to observe how
exact constraint satisfaction can transfer knowledge from minimal training
data onto diverse problems.

Figure 1: ML model architecture and workflow for $(\textbf{a})$ combined model
$(\textbf{b})$ spin-scaled model. For the combined model, total density
$(n_{\uparrow}+n_{\downarrow})$ is given as the input. In the spin-scaled
model, the upper architecture is for the exchange model learning while the
lower is the correlation-like model. Spin-scaling is satisfied in the exchange
model represented by $(2n_{\sigma})$ where $\sigma=\uparrow\downarrow$ spin
channels.

## Neural Network architecture and Training

In this work, all the networks were built on the basis of the ML framework
Tensorflow [62]. For the models, the non-linearity in the mapping is acquired
by using sigmoid[63], tanh and exponential linear unit (elu) [64] activation
functions, chosen as commonly used continuously differentiable activation
functions. Ensuring smooth activation functions was found to be essential for
obtaining reasonable XC potentials, as discussed below. The data set was
randomly divided into a training (80%) set and a validation (20%) set, using
train-test split feature of sci-kit learn [65]. Hyper-parameter searching
identified a 3 layered model with sigmoid to sigmoid to tanh activation
functions as preferable for the combined model while a 2 layered model with
tanh to elu activation functions was preferred for the spin-scaled model. The
network weights and biases were optimized by stochastic gradient descent with
the Adam optimizer [66] using a learning rate of 0.05, applying a gradient
step after each sample in the training data set. The optimized model was
chosen as that which minimizes the error for the validation set, generally
found after one complete pass of the training data.

### Combined Model

Figure 1 a) presents a schematic for the simple neural network(NN)
architecture, termed the “combined model”. As the name suggests, the combined
model receives inputs constructed from total density
($\mathrm{n}_{\uparrow}+\mathrm{n}_{\downarrow}$) as the features, and targets
SCAN’s exchange-correlation enhancement factor. Different numbers of hidden
layers and neuron counts were tested, with a three layered model with 100, 50
and 20 neurons in the respective layers found to perform best.

For the combined model, the loss function to be optimized in the learning
process is defined as,

$\mathcal{L}_{\mathrm{combined}}=\frac{1}{N}\sum_{i}^{N}\left(F_{\mathrm{xc}}^{\mathrm{ML}}-F_{\mathrm{xc}}^{\mathrm{SCAN}}\right)^{2},$
(22)

where $N$ is the number of training data points, thus minimizing the mean
square difference between SCAN XC enhancement and the learned XC
enhancement.The Lieb-Oxford bound for the combined model is introduced as a
post-processing mechanism following the explanation in Eq. 8.

### Spin-Scaled Model

The spin-scaled model follows a more complex architecture that allows it to
obey the spin-scaling exact constraint by treating exchange and correlation
separately, as discussed above. The overall architecture for the spin-scaled
model is presented in Figure 1 b). The spin-scaled model is comprised of two
separate networks, one for exchange and one for correlation. These networks
are trained separately and later combined to form the complete model. Training
is therefore carried out as a two step process.

This separation by spin channel reduces the input domain for the exchange
network to six features suitable for exchange,
{$\tilde{r}_{\mathrm{s}\sigma},\tilde{s}_{\sigma},\tilde{q}_{\sigma},g_{\mathrm{x}\sigma},h_{\mathrm{0}\sigma},h_{\mathrm{1}\sigma}$},
separately generated for each spin $\sigma$. The exchange network is trained
first to minimize the mean square difference with the SCAN exchange
enhancement defined as,

$\mathcal{L}_{\mathrm{exchange}}=\frac{1}{N}\sum_{i}^{N}\left(F_{\mathrm{x}}^{\mathrm{ML}}-F_{\mathrm{x}}^{\mathrm{SCAN}}\right)^{2},$
(23)

with the spin-scaled exchange enhancement,

$F_{\mathrm{x}}^{\mathrm{ML}}=\frac{F_{\mathrm{x}}(2n_{\uparrow})e_{\mathrm{x}}^{\mathrm{LDA}}(2n_{\uparrow})+F_{\mathrm{x}}(2n_{\downarrow})e_{\mathrm{x}}^{\mathrm{LDA}}(2n_{\downarrow})}{2e_{\mathrm{x}}^{\mathrm{LDA}}(n_{\uparrow}+n_{\downarrow})}.$
(24)

As the exchange energy must be invariant to permutation of spin labels, the
same exchange network is used for both spin channels and should be trained on
both spin channels of the training data. The exchange network has two layers
with 80 and 40 neurons at the first and second layers respectively and the
activation functions are tanh and exponential linear unit (elu) [64].

The correlation energy is not subject to the same spin scaling constraint and
is handled by a separate model. This second model takes a reduced set of the
total density ($\mathrm{n}_{\uparrow}+\mathrm{n}_{\downarrow}$) input
variables suitable for correlation:
{$\tilde{r}_{\mathrm{s}},\tilde{s},\tilde{q},\zeta^{2},\epsilon_{\mathrm{c}}^{\mathrm{0}},\epsilon_{\mathrm{c}}^{\mathrm{1}}$}.
This second network has the same hyper-parameters as the exchange network.

Table 1: Deviation of the spin-scaled exchange network total exchange energy (hartrees) from the SCAN exchange functional against which it was trained. Calculated for the set of spherical atom densities and the molecules of the G3 test set. Mean error (ME) and Mean absolute error (MAE) are presented, full data is given in supplemental materials.  |  | No LO Bound | LO Bound
---|---|---|---
Atoms | ME | -0.147 | -0.541
MAE | 0.153 | 0.543
G3 | ME | -0.209 | -0.448
MAE | 0.209 | 0.448

The loss function for the second network is,

$\mathcal{L}_{\mathrm{correlation}}=\frac{1}{N}\sum_{i}^{N}\left[F_{\mathrm{c}}^{\mathrm{ML}}-(F_{\mathrm{xc}}^{\mathrm{SCAN}}-F_{\mathrm{x}}^{\mathrm{ML}})\right]^{2},$
(25)

where $F_{\mathrm{x}}^{\mathrm{ML}}$ is the output of the previously trained
exchange network. This second network is therefore not a true model of SCAN
correlation as the loss function of Eq. 25 drives it to compensate for
deficiencies in the exchange network, though correlation effects will dominate
if the exchange network is accurate. Table 1 shows the deviation of the
bounded and un-bounded exchange networks in total exchange energy from SCAN
for the spherical atoms which were part of the training and test sets, and the
molecules of the G3 set which were not in training set. From this we see that
deviation is typically $<1\%$ for both networks, though the bounded model
deviates more significantly than the unbounded, indicating that good exchange
models have been learned.

Finally, the total enhancement factor is obtained by summing the exchange
network and second network enhancement factors as,

$F_{\mathrm{xc}}^{\mathrm{spin-
scaled}}=F_{\mathrm{x}}^{\mathrm{ML}}+F_{c}^{\mathrm{ML}}.$ (26)

The Lieb-Oxford bounds for the spin-scaled model are introduced separately for
exchange and the correlation parts as they are trained separately. For the
exchange part, we follow similar mechanism as explained in Eq. 8 where the
bound is introduced to the exchange enhancement factor as,

$F_{\mathrm{x}}^{\mathrm{ML-
bound}}(\mathbf{r})=\frac{1.174}{1+\mathrm{ANN}_{\mathrm{x}}(\mathbf{r})^{2}}.$
(27)

Here we choose a tight bound of 1.174 for the exchange.

For correlation, we introduce the bound by modifying the loss function to
include an additional penalty term,

$\mathcal{L}_{\mathrm{correlation}}^{\mathrm{{bound}}}=\frac{1}{N}\sum_{i}^{N}\\{\left[F_{\mathrm{c}}^{\mathrm{ML}}-(F_{\mathrm{xc}}^{\mathrm{SCAN}}-F_{\mathrm{x}}^{\mathrm{ML}})\right]^{2}\\\
+\mu\times\mathrm{relu}[F_{\mathrm{xc}}^{\mathrm{ML}}-2.215]\\},$ (28)

where $\mu$ value is chosen to be 20, and relu is the rectified linear
unit[67]. Here we see that if the total XC enhancement factor is smaller than
the tight bound 2.215 the penalty term is zero whereas any value of total XC
enhancement factor greater than 2.215 will incur penalty. This total loss
function is minimized to satisfy the Lieb–Oxford bound introduced in the ML
model.

## Results and discussion

### Atomic Performance

Figure 2: XC enhancement plots for the test silicon atom which was not
included in the training set for $(\textbf{a})$ combined model and
$(\textbf{b})$ spin-scaled model without Lieb–Oxford bound while
$(\textbf{c})$ and $(\textbf{d})$ represent same as $(\textbf{a})$ and
$(\textbf{b})$ with Lieb-Oxford bound. All models are compared against SCAN
functional. Density was obtained from accurate spherical HF orbitals [57; 58;
59].

Figure 2 shows the SCAN and ML-model XC enhancement for the silicon atom,
which was not part of the training set. Figures a) and b) show combined and
spin-scaled models respectively without the Lieb–Oxford bound, while c), d)
include the Lieb–Oxford bound constraint. All the ML models show good
agreement with the SCAN’s XC enhancement factor for this system. The combined
model without the Lieb–Oxford bound shows slightly less variation in the
energetically important region between r = 0 and r = 4.

During training we did not target the total XC energy directly, in favor of
learning the XC energy density of an existing mGGA instead. This switch
towards learning the XC energy density, a local property, has three benefits.
Most importantly it reduces the complexity of the training by avoiding summing
the derivatives of many training points in a numerical integration batch.
Secondly, far more training data is available for a given functional’s local
energy density than for total XC energies, as every point in the all-space
integration of any system can now be considered a training point. We note that
all training points were weighted equally. As some points are more
energetically important than others this will implicitly bias training towards
energetically important regions, such as core densities, and away from less
important regions, such as asymptotic densities. This could be balanced by
dividing the training weight of each point by the uniform electron gas
exchange energy for the density at the point, but we did not explore this
here. Thirdly, this avoids the introduction of a gauge freedom in which many
different energy density functions can integrate to the same total energy,
which could result in learning a model that gives reasonable total energy, but
poor local accuracy to the SCAN energy density. Despite this training against
local XC density it will be shown that the models are successful in recovering
the total XC energy for the training atom sets, predicting the global property
from local training.

### Molecular Test Sets

For the models trained against data from atomic systems, a real challenge is
to generalize to problems outside the training domain. We examine this by
looking at model accuracy for the open- and closed-shell molecules of the G3
test set [68]. The input ingredients (density, gradient and laplacian) and the
SCAN $F_{\mathrm{xc}}^{\mathrm{SCAN}}$ target for all molecular calculations
were generated from self-consistent SCAN orbitals in the 6-311++G(3df,3pd)
basis set [69; 70]. All molecular calculations were carried using the QUantum
Electronic Structure Techniques (QUEST) program [71].

Table 2: Mean absolute error (MAE) in kcal mol-1 for G3 set of 226 molecular
atomization energies [68]. All SCAN G3 calculations were performed fully self-
consistently with the 6-311++G(3df,3pd) basis set [69; 70] in the QUEST
program[71]. The ML calculations were performed non self-consistently from
SCAN orbitals. The ML models with Lieb–Oxford bound are denoted by “LO” within
table. Errors are given as relative to the total energy calculated using the
parent SCAN functional.

SCAN Comparison
---
ML models | closed shell | open shell | (G3)
| (MAE) | (MAE) | (MAE)
Combined | 7.49 | 7.89 | 7.69
Combined-LO | 11.88 | 9.26 | 10.58
Spin–scaled | 8.03 | 4.85 | 6.44
Spin-scaled-LO | 7.66 | 5.66 | 6.66

Table 3: Mean absolute error (MAE) in kcal mol-1 for G3 set of 226 molecular
atomization energies [68]. All SCAN G3 calculations were performed fully self-
consistently with the 6-311++G(3df,3pd) basis set [69; 70] in the QUEST
program[71]. The ML calculations were performed non self-consistently from
SCAN orbitals. The ML models with Lieb–Oxford bound are denoted by “LO” within
table. Errors are given relative to standard reference values for the G3 set.

Reference Comparison
---
ML models | closed shell | open shell | (G3)
| (MAE) | (MAE) | (MAE)
Combined | 7.80 | 5.17 | 6.48
Combined-LO | 13.24 | 7.26 | 10.25
Spin–scaled | 8.55 | 3.73 | 6.14
Spin-scaled-LO | 11.14 | 6.35 | 8.75

Figure 3: Exchange-correlation enhancement plots for O2 (open shell) and CO
(closed shell) from $(\textbf{a})$ combined model and $(\textbf{b})$ spin-
scaled model without Lieb–Oxford bound. $(\textbf{c})$, $(\textbf{d})$ are
same as $(\textbf{a})$ and $(\textbf{b})$ but for models including Lieb–Oxford
bound. The gray vertical lines mark atomic positions.

The mean absolute error in atomization energy for the G3 set are summarized in
Tables 2 and 3 for open-shell, closed-shell, and total collections. Table 2
shows errors relative to total energies calculated using the parent SCAN
functional to illustrate model fidelity. Table 3 shows errors relative to
standard reference total energies for the set to illustrate the models’
chemical accuracy. Overall all models showed comparable performance giving
accuracy close to the SCAN functional (MAE 6.53 kcal mol-1 ). In particular,
enforcing the correct spin-scaling relation in the ML models improved accuracy
for the open-shell systems, though this was accompanied by a small reduction
in accuracy for the closed shell systems. Introduction of the Lieb–Oxford
bound to the combined model deteriorated performance for both open and closed
shell systems. This unexpected poor performance suggests that enforcing the
Lieb–Oxford bound through Eq. 8 does not seem to be a successful strategy. The
performance of Eq. 27 is more successful, but remains less accurate than the
unconstrained equivalent.

It is possible that the reduced performance found is a result of our locally
imposing the Lieb–Oxford bound on the model. As the Lieb–Oxford bound is a
constraint defined for the total energy imposing it as a constraint on the
integrated energy may be more successful, though this would require a
significantly different training scheme.

While such post-processing is theoretically convenient, the results from the
G3 set and lattice constant below show that it significantly limited model
learning during training. We understand this as an effect of the non-linear
normalization of Eqs. 8 and 27 requiring the raw ANN output,
$\mathrm{ANN(\mathbf{r})}$, to be simultaneously very large when the target
$F_{\mathrm{xc}}(F_{\mathrm{x}})$ is small, and very small in energetically
significant regions where $F_{\mathrm{xc}}(F_{\mathrm{x}})$ approaches $2.215$
(or equivalently $1.174$ for spin-scaled models). The spin-scaled model
performs better than the combined model because the number to return (1.174)
is relatively smaller compared to 2.215 for combined model and hence error is
smaller. An alternative solution for enforcing this constraint could be to
simply truncate the range of the network output,

$F_{\mathrm{xc}}^{\mathrm{ML}}=\max(0,\min(2.215,\mathrm{ANN}(r))),$ (29)

and

$F_{\mathrm{x}}^{\mathrm{ML}}=\max(0,\min(1.174,\mathrm{ANN}(r))),$ (30)

however this may introduce undesirable discontinuities in the partial
derivatives of the model, transferring into a non-physical XC potential.

In order to better understand how faithfully the ML models are reproducing the
SCAN functional, Figure 3 compares the XC enhancement factor,
$F_{\mathrm{xc}}$, for the open shell O2 and closed shell CO molecules. We see
that all the models are accurate in the areas immediately around the nuclei.
This is expected from the good atomic performance as these core regions are
relatively unchanged by the covalent bonding. The models deviate more severely
around the bond center, with all models underestimating
$F_{\mathrm{xc}}^{\mathrm{SCAN}}$. In these regions $|\nabla n|\to 0$ and
hence $s\to 0$, while the density $n$ remains significant. Such regions are
under-represented in the training set, appearing only in small regions at the
center of the compressed Ar2 diatomic. This suggests that the spin scaling
relation enforced in the spin-scaled model as well as the Lieb–Oxford bound
are insufficient to transfer learning from a training set that does not
include chemical bonding, onto systems that are chemically bound. Further
constraint satisfaction, or inclusion of bonding data into the training set,
is likely necessary to improve model accuracy at these important points. We
also note that the combined model with no Lieb–Oxford bound exhibits sharp
spikes in $F_{\mathrm{xc}}$ in the bonding regions that are not seen for other
ML models or the SCAN functional.

Figure 4: Mean Absolute Error in predicted lattice constant

(MAE, Å) for SCAN, whose results are obtained from Ref.[5], SCAN-L, whose
results are obtained from Ref. [14], combined, combined with Lieb–Oxford bound
(denoted “LO” within figure), spin-scaled and spin-scaled with Lieb–Oxford
bound for the LC20 set of 20 solids with lattice constants ranging from 3.415
to 6.042 Å[72]. The ML calculations were performed using fixed self-consistent
PBE densities generated using FHI-aims [73].

### Lattice Constants of solids

Figure 5: Derivatives of SCAN and ML model XC energy density (combined and
spin-scaled) with respect to input ingredients for the silicon atom (not in
training set). ML derivatives are represented by solid lines while SCAN are
dotted.

The transferability of the ML functionals was further tested by calculating
the lattice constants of 20 solids from the LC20 test set[72]. This tests the
model’s ability to generalise into further unseen chemical environments in
periodic systems, as well as requiring description of the energy as a function
of nuclear displacement. The LC20 set is therefore a sensitive test of model
transferability. The equilibrium lattice constants were determined by a nine-
point fit of total energy per unit cell to the Vinet equation of state around
the SCAN equilibrium unit cell volume ($V_{0}$) in a range of $V_{0}\pm 10\%$
[74; 75].

Figure 4 compares LC20 results from ML models with the SCAN and SCAN-L
functionals. The SCAN results for lattice constants are taken from Ref. [5]
and SCAN-L results from Ref. [14]. The ML calculations were obtained non self-
consistently using PBE densities, not SCAN densities as was used for molecular
calculations. This is because SCAN suffered numerical instabilies in
calculations using FHI-aims[73], with some systems failing to converge. FHI-
aims was used for periodic calculations rather than VASP (which was used for
Ref. 5 and 14) as all electron calculations were required to get usable input
data for lattice constant calculations. We should therefore be aware of a
possible density driven error for the ML models resulting from the PBE
densities on which they were evaluated, which could be significant.

The results of Figure 4 show that the ML models do not perform as well for
solids as they do for gas phase atoms and molecules, suggesting difficulty in
generalizing knowledge into this untrained domain. The spin-scaled model
preformed significantly better than the combined model both with and without
enforcing the Lieb–Oxford bound. While the performance of ML models is poor
compared to SCAN and SCAN-L, they improve slightly on the PBE GGA [50] (MAE
0.060Å). Having said that, it suggests a similar conclusion to the molecular
tests; that the spin-scaling and the Lieb–Oxford bound are insufficient to
transfer learning into the periodic systems. Inclusion of such information in
the training set is likely necessary to improve model accuracy.

### Exchange-correlation potential

The XC potential is defined as,

$v_{\mathrm{xc}}(\mathbf{r})=\frac{\delta E_{\mathrm{xc}}}{\delta
n(\mathbf{r})}.$ (31)

This is constructed for the ML models from the partial derivatives of the
model with respect to its ingredients,

$\left\\{\frac{\partial E_{\mathrm{xc}}^{\mathrm{ML}}}{\partial
n(\mathbf{r})},\frac{\partial E_{\mathrm{xc}}^{\mathrm{ML}}}{\partial|\nabla
n(\mathbf{r})|},\frac{\partial
E_{\mathrm{xc}}^{\mathrm{ML}}}{\partial\nabla^{2}n(\mathbf{r})}\right\\},$
(32)

through repeated application of the chain rule. In practice this is achieved
using the back-propagation mechanism of the machine learning framework [76].
Here the choice of the activation functions for the neuron layers is
critically important for obtaining a smooth XC potential appropriate for SCF
applications. When the activation functions are differentiated during back-
propagation the use of activation functions with discontinuous derivatives,
such as the popular rectified linear unit (relu) [77], may introduce
discontinuities into the XC potential that can harm SCF convergence and
computational efficiency. We therefore used only smooth sigmoid, tanh and elu
activation functions within the present models.

Figure 5 shows partial derivatives of the combined and spin-scaled models with
respect to input ingredients, compared against the equivalent for the SCAN
functional where they exist. Note that
$\epsilon_{\mathrm{xc}}^{\mathrm{SCAN}}(n,|\nabla n|,\tau)$ and
$\epsilon_{\mathrm{xc}}^{\mathrm{ML}}(n,|\nabla n|,\nabla^{2}n)$ are
necessarily different functions, even though
$\epsilon_{\mathrm{xc}}^{\mathrm{ML}}$ has been trained to reproduce
$\epsilon_{\mathrm{xc}}^{\mathrm{SCAN}}$. Hence, we should not expect their
partial derivatives in $n$ and $|\nabla n|$ to match. Figure 5 (a-d) shows the
combined model partial derivatives for the test silicon atom, while (e-h)
shows the same for the spin-scaled model.

Figures 5 (a), (e), (i) and (m) show that the density partial derivative is
comparable to SCAN for both models. Both models exhibit oscillations in this
derivative, though these are less severe for the spin-scaled models than the
combined models. Figures 5 (b-c) and (j-k) for the combined model without and
with Lieb–Oxford bound constraint and (f-g) and (n-o) for spin-scaled model,
again with and without Lieb–Oxford bound, show the partial derivatives with
respect to the same spin and cross spin gradient components. While the
combined model does not distinguish between different spins
$\partial\epsilon_{\mathrm{xc}^{\mathrm{ML}}}/\partial\nabla
n_{\sigma}^{2}=\partial\epsilon_{\mathrm{xc}^{\mathrm{ML}}}/\partial\nabla
n_{\uparrow}\cdot\nabla n_{\downarrow}$, the spin-scaled model is able to
distinguish spin channels so the same-spin and cross-spin partial derivatives
are different. Finally, Figures (d), (h), (i) and (k)show the Laplacian
partial derivatives.

It is not clear to what degree the oscillations of the present partial
derivatives may affect the SCF performance of the ML models without performing
such calculations, which is beyond the scope of the present study. We can
reasonably expect the smoother spin-scaled model to outperform the combined
model in this regard however, as a result of the reduced oscillations seen in
Figure 5. The harmful impact of such oscillations may be significantly damped
if the regions they occur in are not energetically important. We note that the
SCAN functional is known to have problematic oscillations in its XC potential
which reduce its computational efficiency but generally do not prevent SCF
convergence [78; 79; 80; 81; 82]. The similarity between the ML and SCAN
partial derivatives of Figure 5 therefore suggest that self-consistency can
likely be reached, though this may be sensitive to the choice of starting
guess orbitals if they are different from the converged orbitals used in
training.

## Discussion

A broad goal of this work has been to explore how exact constraint adherence
can enhance the transferability of ML density functionals beyond a limited set
of training data for which their accuracy is theoretically well justified. For
the semi-local meta-GGA functionals these are systems for which the XC hole
remains well localized (atoms and compressed non-bonding diatoms), or for
which the long-range parts of the X and C holes cancel (jellium) [83; 84]. It
is unclear to what extent our models’ overall fidelity to SCAN is a limitation
of this restricted training data, and to what extent it reflects a fundamental
difficulty in representing SCAN from orbital free ingredients. The success of
the SCAN-L and related functionals [13] suggests the former, though repeating
the present study with a significantly expanded training set that covers more
of the chemical space could provide compelling evidence one way or the other,
though is beyond the present scope.

## Conclusion

We have explored how a philosophy of exact constraints and appropriate norms
can be combined with ML techniques in functional design. We have shown a
simple test of this idea as a de-orbitalization of the SCAN functional, using
the density Laplacian, $\nabla^{2}{n(\mathbf{r})}$, in place of kinetic energy
density $\tau(\mathbf{r})$. Four ML functionals were developed, enforcing a
spin-scaling constraint, the Lieb–Oxford bound, both, and neither. These
models were trained using a dataset with no chemical bonding, preferring the
norms appropriate for semi-local functionals [5]. The model satisfying both
the spin-scaling constraint and the Lieb–Oxford bound generally achieved a
more balanced performance across the properties tested, though performance was
worse than that achieved by the analytical SCAN-L deorbitalization for solids.
Given that model performance was generally improved when both constraints were
imposed, it is reasonable to believe that engineering in further constraints
can enhance robustness of the ML models out of the training domain.

## Data Availability

The data that supports the finding of this study are available from the
corresponding author upon request. Parameters for the ML models presented are
provided as supplemental Tensorflow metagraph parameter files.

## Acknowledgments

J.F. thanks Michael Withnall for his useful discussion and guidance on ML
techniques. J.F., K.P., and J.S. acknowledge discussions with John Perdew on
semi-local XC functionals. J.S. thanks Kieron Burke for his useful discussion
on ML techniques. This work was partially supported by the Research Council of
Norway through its Centres of Excellence scheme, project number 262695. We
acknowledge financial support from the European Research Council under
H2020/ERC Consolidator Grant top DFT (Grant No. 772259).

## Author contributions

J.S and J.F. designed the project. K.P., J.F., Y.Y., V.B., T.I., and A.T.
developed the computational methods, K.P. performed all calculations. K.P.,
J.F., and J.S. analyzed the data and wrote the manuscript.

## competing interests

The authors declare no competing interests.

## References

  * [1] P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Physical review, vol. 136, no. 3B, p. B864, 1964.
  * [2] W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Physical review, vol. 140, no. 4A, p. A1133, 1965\.
  * [3] J. P. Perdew and K. Schmidt, “Jacob’s ladder of density functional approximations for the exchange-correlation energy,” in AIP Conference Proceedings, vol. 577, pp. 1–20, American Institute of Physics, 2001.
  * [4] J. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao, A. Ruzsinszky, G. I. Csonka, G. E. Scuseria, and J. P. Perdew, “Density functionals that recognize covalent, metallic, and weak bonds,” Phys. Rev. Lett., vol. 111, no. 10, p. 106401, 2013.
  * [5] J. Sun, A. Ruzsinszky, and J. P. Perdew, “Strongly constrained and appropriately normed semilocal density functional,” Physical review letters, vol. 115, no. 3, p. 036402, 2015.
  * [6] J. W. Furness and J. Sun, “Enhancing the efficiency of density functionals with an improved iso-orbital indicator,” Physical Review B, vol. 99, no. 4, p. 041119, 2019.
  * [7] A. D. Becke, “A new inhomogeneity parameter in density-functional theory,” The Journal of chemical physics, vol. 109, no. 6, pp. 2092–2098, 1998.
  * [8] A. Görling and M. Levy, “Exact kohn-sham scheme based on perturbation theory,” Physical Review A, vol. 50, no. 1, p. 196, 1994.
  * [9] Q. Wu and W. Yang, “Algebraic equation and iterative optimization for the optimized effective potential in density functional theory,” Journal of Theoretical and Computational Chemistry, vol. 2, no. 04, pp. 627–638, 2003.
  * [10] A. Seidl, A. Görling, P. Vogl, J. Majewski, and M. Levy, “Generalized Kohn-Sham schemes and the band-gap problem,” Phys. Rev. B, vol. 53, no. 7, pp. 3764–3774, 1996.
  * [11] R. Neumann, R. H. Nobes, and N. C. Handy, “Exchange functionals and potentials,” Mol. Phys., vol. 87, no. 1, pp. 1–36, 1996.
  * [12] Q. Wu and W. Yang, “A direct optimization method for calculating density functionals and exchange–correlation potentials from electron densities,” The Journal of chemical physics, vol. 118, no. 6, pp. 2498–2509, 2003.
  * [13] D. Mejia-Rodriguez and S. Trickey, “Deorbitalization strategies for meta-generalized-gradient-approximation exchange-correlation functionals,” Physical Review A, vol. 96, no. 5, p. 052512, 2017.
  * [14] D. Mejia-Rodriguez and S. Trickey, “Deorbitalized meta-gga exchange-correlation functionals in solids,” Physical Review B, vol. 98, no. 11, p. 115161, 2018.
  * [15] F. Della Sala, E. Fabiano, and L. A. Constantin, “Kinetic-energy-density dependent semilocal exchange-correlation functionals,” International Journal of Quantum Chemistry, vol. 116, no. 22, pp. 1641–1694, 2016.
  * [16] M. Dorigo and U. Schnepf, “Genetics-based machine learning and behavior-based robotics: a new synthesis,” IEEE Transactions on Systems, Man, and Cybernetics, vol. 23, no. 1, pp. 141–154, 1993.
  * [17] M. A. Salichs, R. Barber, A. M. Khamis, M. Malfaz, J. F. Gorostiza, R. Pacheco, R. Rivas, A. Corrales, E. Delgado, and D. Garcia, “Maggie: A robotic platform for human-robot social interaction,” in 2006 IEEE conference on robotics, automation and mechatronics, pp. 1–7, IEEE, 2006.
  * [18] N. MacLeod, M. Benfield, and P. Culverhouse, “Time to automate identification,” Nature, vol. 467, no. 7312, pp. 154–155, 2010.
  * [19] L. Chittka and A. Dyer, “Your face looks familiar,” Nature, vol. 481, no. 7380, pp. 154–155, 2012.
  * [20] E. A. Ashley, “The precision medicine initiative: a new national effort,” Jama, vol. 313, no. 21, pp. 2119–2120, 2015.
  * [21] M. C. Schatz and B. Langmead, “The dna data deluge,” Ieee Spectrum, vol. 50, no. 7, pp. 28–33, 2013.
  * [22] H. Ding, I. Takigawa, H. Mamitsuka, and S. Zhu, “Similarity-based machine learning methods for predicting drug–target interactions: a brief review,” Briefings in bioinformatics, vol. 15, no. 5, pp. 734–747, 2014.
  * [23] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, et al., “Mastering the game of go with deep neural networks and tree search,” nature, vol. 529, no. 7587, pp. 484–489, 2016.
  * [24] J. C. Snyder, M. Rupp, K. Hansen, K.-R. Müller, and K. Burke, “Finding density functionals with machine learning,” Physical review letters, vol. 108, no. 25, p. 253002, 2012.
  * [25] J. C. Snyder, M. Rupp, K. Hansen, L. Blooston, K.-R. Müller, and K. Burke, “Orbital-free bond breaking via machine learning,” The Journal of chemical physics, vol. 139, no. 22, p. 224104, 2013.
  * [26] F. Brockherde, L. Vogt, L. Li, M. E. Tuckerman, K. Burke, and K.-R. Müller, “Bypassing the kohn-sham equations with machine learning,” Nature communications, vol. 8, no. 1, pp. 1–10, 2017.
  * [27] Q. Liu, J. Wang, P. Du, L. Hu, X. Zheng, and G. Chen, “Improving the performance of long-range-corrected exchange-correlation functional with an embedded neural network,” The Journal of Physical Chemistry A, vol. 121, no. 38, pp. 7273–7281, 2017.
  * [28] R. Nagai, R. Akashi, S. Sasaki, and S. Tsuneyuki, “Neural-network kohn-sham exchange-correlation potential and its out-of-training transferability,” The Journal of chemical physics, vol. 148, no. 24, p. 241737, 2018.
  * [29] J. Wellendorff, K. T. Lundgaard, A. Møgelhøj, V. Petzold, D. D. Landis, J. K. Nørskov, T. Bligaard, and K. W. Jacobsen, “Density functionals for surface science: Exchange-correlation model development with bayesian error estimation,” Physical Review B, vol. 85, no. 23, p. 235149, 2012.
  * [30] L. Li, T. E. Baker, S. R. White, K. Burke, et al., “Pure density functional for strong correlation and the thermodynamic limit from machine learning,” Physical Review B, vol. 94, no. 24, p. 245129, 2016.
  * [31] J. Westermayr, M. Gastegger, K. T. Schütt, and R. J. Maurer, “Perspective on integrating machine learning into computational chemistry and materials science,” arXiv preprint arXiv:2102.08435, 2021.
  * [32] J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, et al., “Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional,” Nature chemistry, vol. 8, no. 9, p. 831, 2016.
  * [33] M. Chen, H.-Y. Ko, R. C. Remsing, M. F. C. Andrade, B. Santra, Z. Sun, A. Selloni, R. Car, M. L. Klein, J. P. Perdew, et al., “Ab initio theory and modeling of water,” Proceedings of the National Academy of Sciences, vol. 114, no. 41, pp. 10846–10851, 2017.
  * [34] R. C. Remsing, M. L. Klein, and J. Sun, “Dependence of the structure and dynamics of liquid silicon on the choice of density functional approximation,” Physical Review B, vol. 96, no. 2, p. 024203, 2017.
  * [35] G. S. Gautam and E. A. Carter, “Evaluating transition metal oxides within dft-scan and scan+ u frameworks for solar thermochemical applications,” Physical Review Materials, vol. 2, no. 9, p. 095401, 2018.
  * [36] J. W. Furness, Y. Zhang, C. Lane, I. G. Buda, B. Barbiellini, R. S. Markiewicz, A. Bansil, and J. Sun, “An accurate first-principles treatment of doping-dependent electronic structure of high-temperature cuprate superconductors,” Communications Physics, vol. 1, no. 1, pp. 1–6, 2018\.
  * [37] C. Lane, J. W. Furness, I. G. Buda, Y. Zhang, R. S. Markiewicz, B. Barbiellini, J. Sun, and A. Bansil, “Antiferromagnetic ground state of la 2 cuo 4: A parameter-free ab initio description,” Physical Review B, vol. 98, no. 12, p. 125140, 2018.
  * [38] K. Pokharel, C. Lane, J. W. Furness, R. Zhang, J. Ning, B. Barbiellini, R. S. Markiewicz, Y. Zhang, A. Bansil, and J. Sun, “Sensitivity of the electronic and magnetic structures of high-temperature cuprate superconductors to exchange-correlation density functionals,” arXiv preprint arXiv:2004.08047, 2020.
  * [39] J. Hollingsworth, L. Li, T. E. Baker, and K. Burke, “Can exact conditions improve machine-learned density functionals?,” The Journal of chemical physics, vol. 148, no. 24, p. 241743, 2018.
  * [40] R. Nagai, R. Akashi, and O. Sugino, “Machine-learning-based exchange correlation functional with physical asymptotic constraints,” Physical Review Research, vol. 4, no. 1, p. 013106, 2022.
  * [41] J. Kirkpatrick, B. McMorrow, D. H. Turban, A. L. Gaunt, J. S. Spencer, A. G. Matthews, A. Obika, L. Thiry, M. Fortunato, D. Pfau, et al., “Pushing the frontiers of density functionals by solving the fractional electron problem,” Science, vol. 374, no. 6573, pp. 1385–1389, 2021.
  * [42] E. H. Lieb and S. Oxford, “Improved lower bound on the indirect coulomb energy,” International Journal of Quantum Chemistry, vol. 19, no. 3, pp. 427–439, 1981.
  * [43] G. Oliver and J. Perdew, “Spin-density gradient expansion for the kinetic energy,” Physical Review A, vol. 20, no. 2, p. 397, 1979.
  * [44] M. Levy and J. P. Perdew, “Hellmann-feynman, virial, and scaling requisites for the exact universal density functionals. shape of the correlation potential and diamagnetic susceptibility for atoms,” Physical Review A, vol. 32, no. 4, p. 2010, 1985.
  * [45] P.-S. Svendsen and U. von Barth, “Gradient expansion of the exchange energy from second-order density response theory,” Physical Review B, vol. 54, no. 24, p. 17402, 1996.
  * [46] L. Pollack and J. Perdew, “Evaluating density functional performance for the quasi-two-dimensional electron gas,” Journal of Physics: Condensed Matter, vol. 12, no. 7, p. 1239, 2000.
  * [47] J. P. Perdew and Y. Wang, “Accurate and simple analytic representation of the electron-gas correlation energy,” Physical review B, vol. 45, no. 23, p. 13244, 1992.
  * [48] J. P. Perdew, A. Ruzsinszky, J. Sun, and K. Burke, “Gedanken densities and exact constraints in density functional theory,” The Journal of chemical physics, vol. 140, no. 18, p. 18A533, 2014.
  * [49] J. P. Perdew and Y. Wang, “Erratum: Accurate and simple analytic representation of the electron-gas correlation energy [phys. rev. b 45, 13244 (1992)],” Physical Review B, vol. 98, no. 7, p. 079904, 2018.
  * [50] J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized Gradient Approximation Made Simple.,” Phys. Rev. Lett., vol. 77, no. 18, pp. 3865–3868, 1996\.
  * [51] J. P. Perdew, P. Ziesche, and H. Eschrig, “Electronic structure of solids’ 91,” 1991.
  * [52] M. Seidl, J. P. Perdew, and S. Kurth, “Simulation of all-order density-functional perturbation theory, using the second order and the strong-correlation limit,” Physical review letters, vol. 84, no. 22, p. 5070, 2000.
  * [53] J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, “Climbing the density functional ladder: Nonempirical meta–generalized gradient approximation designed for molecules and solids,” Physical Review Letters, vol. 91, no. 14, p. 146401, 2003.
  * [54] J. Tao, J. P. Perdew, L. M. Almeida, C. Fiolhais, and S. Kümmel, “Nonempirical density functionals investigated for jellium: Spin-polarized surfaces, spherical clusters, and bulk linear response,” Physical Review B, vol. 77, no. 24, p. 245107, 2008.
  * [55] B. L. Kalman and S. C. Kwasny, “Why tanh: choosing a sigmoidal function,” in [Proceedings 1992] IJCNN International Joint Conference on Neural Networks, vol. 4, pp. 578–581, IEEE, 1992.
  * [56] G. L. Oliver and J. P. Perdew, “Spin-density gradient expansion for the kinetic energy,” Phys. Rev. A, vol. 20, no. 2, pp. 397–403, 1979.
  * [57] E. Clementi and C. Roetti, “Roothaan-Hartree-Fock Atomic Wavefunctions,” At. Data Nucl. Data Tables, vol. 14, no. 3, pp. 177–478, 1974.
  * [58] T. Koga, K. Kanayama, S. Watanabe, and A. J. Thakkar, “Analytical Hartree-Fock wave functions subject to cusp and asymptotic constraints: He to Xe, Li+ to Cs+, H- to I-,” Int. J. Quantum Chem., vol. 71, no. 6, pp. 491–497, 1999.
  * [59] J. W. Furness and S. Lehtola, “Hartree-Fock Orbitals for Spherical Atoms - a python toolbox,” 2021.
  * [60] L. M. Almeida, J. P. Perdew, and C. Fiolhais, “Surface and curvature energies from jellium spheres: Density functional hierarchy and quantum Monte Carlo,” Phys. Rev. B, vol. 66, no. 7, p. 075115, 2002.
  * [61] B. Wood, N. D. M. Hine, W. M. C. Foulkes, and P. García-González, “Quantum Monte Carlo calculations of the surface energy of an electron gas,” Phys. Rev. B, vol. 76, no. 3, p. 035403, 2007.
  * [62] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al., “Tensorflow: A system for large-scale machine learning,” in 12th $\\{$USENIX$\\}$ symposium on operating systems design and implementation ($\\{$OSDI$\\}$ 16), pp. 265–283, 2016\.
  * [63] J. Han and C. Moraga, “The influence of the sigmoid function parameters on the speed of backpropagation learning,” in International Workshop on Artificial Neural Networks, pp. 195–201, Springer, 1995.
  * [64] D.-A. Clevert, T. Unterthiner, and S. Hochreiter, “Fast and accurate deep network learning by exponential linear units (elus),” arXiv preprint arXiv:1511.07289, 2015.
  * [65] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” Journal of Machine Learning Research, vol. 12, pp. 2825–2830, 2011.
  * [66] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [67] V. Nair and G. E. Hinton, “Rectified linear units improve restricted boltzmann machines,” in Icml, 2010.
  * [68] L. A. Curtiss, K. Raghavachari, P. C. Redfern, and J. A. Pople, “Assessment of Gaussian-3 and density functional theories for a larger experimental test set,” J. Chem. Phys., vol. 112, no. 17, pp. 7374–7383, 2000.
  * [69] T. Clark, J. Chandrasekhar, G. W. Spitznagel, and P. V. R. Schleyer, “Efficient diffuse function-augmented basis sets for anion calculations. iii. the 3-21+ g basis set for first-row elements, li–f,” Journal of Computational Chemistry, vol. 4, no. 3, pp. 294–301, 1983.
  * [70] M. J. Frisch, J. A. Pople, and J. S. Binkley, “Self-consistent molecular orbital methods 25. supplementary functions for gaussian basis sets,” The Journal of chemical physics, vol. 80, no. 7, pp. 3265–3269, 1984.
  * [71] A. M. Teale, J. W. Furness, T. J. P. Irons, M. S. Ryley, J. Zemen, M. Withnal, and J. Elsdon, “QUantum Electronic Structure Techniques (QUEST),” 2017.
  * [72] J. Sun, M. Marsman, G. I. Csonka, A. Ruzsinszky, P. Hao, Y. S. Kim, G. Kresse, and J. P. Perdew, “Self-consistent meta-generalized gradient approximation within the projector-augmented-wave method,” Phys. Rev. B, vol. 84, no. 3, p. 035117, 2011.
  * [73] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, “Ab initio molecular simulations with numeric atom-centered orbitals,” Computer Physics Communications, vol. 180, no. 11, pp. 2175–2196, 2009.
  * [74] P. Vinet, J. Ferrante, J. Smith, and J. Rose, “A universal equation of state for solids,” Journal of Physics C: Solid State Physics, vol. 19, no. 20, p. L467, 1986.
  * [75] P. Vinet, J. H. Rose, J. Ferrante, and J. R. Smith, “Universal features of the equation of state of solids,” Journal of Physics: Condensed Matter, vol. 1, no. 11, p. 1941, 1989.
  * [76] H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, “Differentiable programming tensor networks,” Physical Review X, vol. 9, no. 3, p. 031041, 2019.
  * [77] J. Brownlee, “A gentle introduction to the rectified linear unit (relu),” Machine learning mastery, vol. 6, 2019.
  * [78] Z.-h. Yang, H. Peng, J. Sun, and J. P. Perdew, “More realistic band gaps from meta-generalized gradient approximations: Only in a generalized Kohn-Sham scheme,” Phys. Rev. B, vol. 93, no. 20, p. 205205, 2016.
  * [79] Y. Yao and Y. Kanai, “Plane-wave pseudopotential implementation and performance of SCAN meta-GGA exchange-correlation functional for extended systems,” J. Chem. Phys., vol. 146, no. 22, p. 224105, 2017.
  * [80] Y. Yamamoto, C. M. Diaz, L. Basurto, K. A. Jackson, T. Baruah, and R. R. Zope, “Fermi-Löwdin orbital self-interaction correction using the strongly constrained and appropriately normed meta-GGA functional,” J. Chem. Phys., vol. 151, no. 15, p. 154105, 2019.
  * [81] A. P. Bartók and J. R. Yates, “Regularized SCAN functional,” J. Chem. Phys., vol. 150, no. 16, p. 161101, 2019.
  * [82] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, “Accurate and numerically efficient r2scan meta-generalized gradient approximation,” The journal of physical chemistry letters, vol. 11, no. 19, pp. 8208–8215, 2020\.
  * [83] J. M. Pitarke and J. P. Perdew, “Metal surface energy: Persistent cancellation of short-range correlation effects beyond the random phase approximation,” Phys. Rev. B, vol. 67, p. 045101, Jan 2003.
  * [84] D. C. Langreth and J. P. Perdew, “Exchange-correlation energy of a metallic surface: Wave-vector analysis,” Phys. Rev. B, vol. 15, pp. 2884–2901, Mar 1977.