# Physics-enhanced neural networks for equation-of-state calculations

Timothy J. Callow<EMAIL_ADDRESS>Center for Advanced Systems Understanding
(CASUS), D-02826 Görlitz, Germany  Helmholtz-Zentrum Dresden-Rossendorf,
D-01328 Dresden, Germany  Jan Nikl<EMAIL_ADDRESS>Center for Advanced Systems
Understanding (CASUS), D-02826 Görlitz, Germany  Helmholtz-Zentrum Dresden-
Rossendorf, D-01328 Dresden, Germany  Eli Kraisler
<EMAIL_ADDRESS>Fritz Haber Center for Molecular Dynamics and
Institute of Chemistry, The Hebrew University of Jerusalem, 9091401 Jerusalem,
Israel  Attila Cangi<EMAIL_ADDRESS>Center for Advanced Systems
Understanding (CASUS), D-02826 Görlitz, Germany  Helmholtz-Zentrum Dresden-
Rossendorf, D-01328 Dresden, Germany

###### Abstract

Rapid access to accurate equation-of-state (EOS) data is crucial in the warm-
dense matter regime, as it is employed in various applications, such as
providing input for hydrodynamic codes to model inertial confinement fusion
processes. In this study, we develop neural network models for predicting the
EOS based on first-principles data. The first model utilizes basic physical
properties, while the second model incorporates more sophisticated physical
information, using output from average-atom calculations as features. Average-
atom models are often noted for providing a reasonable balance of accuracy and
speed; however, our comparison of average-atom models and higher-fidelity
calculations shows that more accurate models are required in the warm-dense
matter regime. Both the neural network models we propose, particularly the
physics-enhanced one, demonstrate significant potential as accurate and
efficient methods for computing EOS data in warm-dense matter.

###### Contents

  1. I Introduction
  2. II Average-atom models: theoretical background
  3. III Pressure in average-atom models
    1. III.1 Functional derivative of the free energy
    2. III.2 Stress-tensor
    3. III.3 Virial theorem
    4. III.4 Ideal approximation
  4. IV Neural network methodology
    1. IV.1 Feature selection and scaling
    2. IV.2 Network architecture and hyper-parameter optimization
  5. V Results
    1. V.1 FPEOS database of Militzer _et al._ : Average-atom results
    2. V.2 FPEOS database of Militzer _et al._ : Neural network results
    3. V.3 FPEOS Be data of Ding and Hu: Average-atom and neural network results
  6. VI Conclusions
  7. A Comparison of neural network models with FPEOS direct interpolation

##  I Introduction

The warm-dense matter (WDM) regime has emerged as a topic of great interest in
recent years [1, 2, 3, 4, 5, 6] for several reasons. Firstly, a plethora of
interesting scientific phenomena occur under WDM conditions: of particular
importance, given the present requirement for clean and abundant energy, is
inertial confinement fusion (ICF) [7, 8, 9, 10]. WDM conditions also present
themselves in a variety of astrophysical systems, such as stars in various
stages of their life cycle [11, 12, 13], planetary cores [14, 15], and more
besides [16]. Secondly, it has become possible to probe matter exposed to WDM
conditions in the laboratory [4, 17, 18, 19, 20, 21], due to the increasing
capabilities of large experimental machines such as the National Ignition
Facility [22], the Linac Coherent Light Source [23], SACLA [24], and the
European XFEL [25]. Thirdly, the WDM regime is notoriously difficult to model,
in part due to the unique position it occupies between the condensed-matter
physics and plasma physics domains [6, 5, 26].

Equation-of-state (EOS) data — the pressure and energy of a material as a
function of its temperature and mass density — is of particular importance in
WDM. For example, hydrodynamic codes that are used to guide ICF experiments
rely on accurate EOS data, typically given in tabular form and interpolated to
the conditions of interest, to close the conservation equations [27]. Many
theoretical techniques exists to compute EOS data. These include “first-
principles” methods, namely density-functional theory molecular-dynamics (DFT-
MD) [28, 29] and path integral Monte Carlo (PIMC) [28, 30], as well as reduced
models such as average-atom (AA) models [31, 32, 33], and extensions thereof
[34, 35]. A review of the application of these methods to EOS calculations can
be found in Ref. [27]. Whilst these models differ in both accuracy and
computational expense, they share the common principle of having only some
fundamental physical properties, i.e. the temperature, density, and material
composition, as inputs to compute the EOS. This is in contrast to alternative
approaches such as the widely-used SESAME database [36], which extrapolates
EOS data from cold curves with a variety of different models.

Assessing the accuracy of EOS models is not straightforward. One of the
principle difficulties in bench-marking models to experimental data is that
temperature cannot be directly measured under WDM conditions [17, 37, 38].
Historically, the temperature has been inferred from other quantities which
can be measured, such as the internal energy [17], and is thus dependent on
the model used for the inference; recent developments in this area hold future
promise to resolve this issue [39, 38, 40]. Therefore, whilst experimental EOS
data exists, it is in short supply and not always suited for comparison with
theoretical models [17].

Consequently, results from the first-principles methods, DFT-MD and PIMC, are
typically taken as the highest quality benchmarks for EOS data in the WDM
regime. These methods are denoted “first-principles” because, in theory, they
exactly describe the many-body electronic structure problem. DFT-MD actually
incorporates two different theories, Kohn–Sham DFT-MD (KS-DFT-MD) [41], which
is often synonymous with DFT-MD due to its more widespread use, and orbital-
free DFT-MD (OF-DFT-MD) [42]. Whilst both flavours are in principle exact,
they rely on approximations. In KS-DFT, we have to approximate the exchange-
correlation (xc) energy functional, for which there exists a wide spectrum of
choices [43, 44, 45], with efforts ongoing to assess the impact of this choice
in the WDM regime [46, 47, 48, 49, 50, 51, 26, 52, 53]. Meanwhile, OF-DFT
relies on approximations for both the kinetic energy density-functional in
addition to the xc-energy [54, 55, 56, 57, 58]. Likewise, PIMC calculations of
real materials have historically relied on the fixed-node approximation [59]
to circumvent the fermion sign problem [60] (although PIMC with no such
approximation has recently been applied to study Hydrogen under WDM conditions
[61]). Of course, both methods also rely on other fundamental assumptions
(such as the Born–Oppenheimer approximation), besides various numerical
approximations.

In recent years, (at least) two substantial EOS databases constructed from
first-principles data have been published. The database of Militzer and co-
workers [28] is composed of KS-DFT (at low temperature) and PIMC (at high
temperatures) data. It includes data from a range of elements and mixtures,
and is constructed from a series of historical papers. The database of Ding
and Hu [29] is composed of KS-DFT (at low temperatures) and OF-DFT (at high
temperatures) data. In both databases, it was established that the different
methods used yielded results in close agreement in the overlap region.

Despite the approximations used, DFT-MD and PIMC are trusted methods for
computation of EOS data and thus comparisons with these methods are typically
made when assessing the accuracy of simpler models, such as AA models.
However, we have typically observed that these comparisons generally involve a
relatively small set of benchmark data points (say in the region of tens of
data points), and the analysis is often a visual comparison of the pressures
on a log-log plot. Historically, this was most likely due to the lack of large
datasets with which to make comparisons. However, the EOS databases just
discussed [28, 29] allow for a more comprehensive analysis of the performance
of AA models.

A second consequence of the existence of larger databases is that they offer
greater possibilities for the application of machine learning techniques such
as neural networks. The use of neural networks in scientific fields such as
materials science [62, 63], quantum physics and chemistry [64, 65, 66, 67,
68], high-energy-density physics [69], and even in warm dense matter [70, 71,
72], is becoming increasingly popular. There are many reasons behind this, but
relevant to this paper (besides the aforementioned growth in the size and
number of databases) is the fact that they are excellent function
approximators, with the ability to learn complex non-linear relationships
between benchmark data and input features [73, 74].

Neural network models can be trained on first-principles data and then applied
to arbitrary densities and temperatures within the range of the training data.
Of course, interpolating EOS data is a well-established practise. However,
machine learning models offer an intriguing path to interpolate EOS data with
only minimal physical constraints from the user, and, for example, modern
practises in machine learning can be used for more robust error estimation.
Applying machine learning techniques to EOS data is still a relatively fresh
topic: most of the applications to date focussed on uncertainty quantification
[75, 76, 77], although in a recent work [78] the authors built a surrogate EOS
model using neural networks.

In this paper, our primary objective is to develop a neural network model
capable of accurately interpolating the EOS by leveraging existing first-
principles datasets [28, 29]. Furthermore, we explore the potential of
physics-enhanced neural networks for improving model performance, by
incorporating additional physical information from an AA model, which has a
low computational overhead for data generation. We also show that these neural
network models yield a significant improvement in performance relative to the
base AA model.

The paper is structured as follows. Firstly, in Section II, we give a brief
overview of the AA model which we use. Next, in Section III, we present
various methods for computing the pressure in AA models, which will all be
compared in the results. After that, in Section IV, we explain the method that
we use to train and evaluate our neural networks. Finally, in Section V, we
compare results from the AA and neural network models against the two first-
principles datasets [28, 29]. We note that this section is split into separate
parts: in the first part, we evaluate the performance of the networks on the
FPEOS database of Ref. [28], which was the database used to train the networks
111Of course, the performance of the network is evaluated on data that was not
seen during training; the point is that the training and test data is drawn
from the same database. In the second part, we evaluate the performance of the
networks on the database of Ref. [29], which was not used in the training of
the neural networks. This second part therefore tests the robustness of the
networks with respect to different sources of data, and also tests their
ability to interpolate to unseen species.

##  II Average-atom models: theoretical background

Average-atom (AA) models have a long history in plasma physics [80, 81, 31,
82, 83, 26], especially in terms of calculation EOS data [84, 32, 85, 33, 35].
There exists a broad range of AA models which use different assumptions and
approximations. However, AA models are united by a common concept, namely the
reduction of a many-body fully-interacting system of electrons and nuclei to
an effective single atom model. The effect of inter-atomic interactions is
coarsely accounted for via the boundary conditions used to solve the
Schrödinger equation for the single atom. A derivation from first-principles
of an AA model similar to the one used in this work can be found in Ref. [26].

We use the AA model proposed by Massacrier _et al_ [86]. In this model, we
solve the Kohn–Sham equations in a spherically symmetric potential,

$\left[\frac{\textrm{d}^{2}}{\textrm{d}r^{2}}+\frac{2}{r}\frac{\textrm{d}}{\textrm{d}r}-\frac{l(l+1)}{r^{2}}\right]X_{\epsilon
nl}(r)+2\left[\epsilon_{nl}-v_{\textrm{s}}[n](r)\right]X_{\epsilon
nl}(r)=0\,.$ (1)

In the above, $X_{\epsilon nl}$ are the radial KS orbitals, and
$\epsilon_{nl}$ their associated energies. The $n,l$ indices in the subscripts
denote the principle and azimuthal quantum numbers. Because the average-atom
represents an extended system, rather than an isolated atom, there is an
energy band (for non-core states) for each pair of $n,l$ values. Within this
energy band, there is a continuous set of KS orbitals with energy $\epsilon$,
which is why the KS orbitals $X_{\epsilon nl}$ have an additional subscript
$\epsilon$.

In the KS equation (1), $v_{\textrm{s}}[n](r)$ is the KS potential, given by

$v_{\textrm{s}}[n](r)=-\frac{Z}{r}+4\pi\int_{0}^{R_{\textrm{VS}}}\textrm{d}{x}\frac{n(x)x^{2}}{\max{(r,x)}}+\frac{\delta
F_{\textrm{xc}}[n]}{\delta n(r)}\,,$ (2)

The three terms in the potential are respectively the electron-nuclear
attraction, the classical Hartree repulsion, and the exchange-correlation (xc)
potential, which is equal to the functional derivative of the xc free energy
$F_{\textrm{xc}}[n]$. In this work, we use exclusively the ground-state (no
temperature dependence) local density approximation (LDA) for
$F_{\textrm{xc}}[n]$ [41], using the parameterization of Perdew and Wang [87].
As ever, due to the dependence of the KS potential on the density $n(r)$, the
KS equations must be solved iteratively until self-consistency is reached.

The KS equations (1) are first solved under the following two boundary
conditions,

$\displaystyle 0$ $\displaystyle=X_{\epsilon nl}(R_{\textrm{VS}})\,,$ (3)
$\displaystyle 0$ $\displaystyle=\frac{\textrm{d}X_{\epsilon
nl}(r)}{\textrm{d}r}\Bigg{|}_{r=R_{\textrm{VS}}}\,,$ (4)

where $R_{\textrm{VS}}$ is the atomic radius (or the Voronoi sphere radius),
defined in terms of the ion number density $n_{i}$ as

$R_{\textrm{VS}}=\left(\frac{3}{4\pi n_{i}}\right)^{1/3}\,.$ (5)

Solving under the above boundary conditions yields a set of eigen-functions
$X^{+}_{nl}(r)$ and $X^{-}_{nl}(r)$ with associated eigen-energies
$\epsilon_{nl}^{+}$ and $\epsilon_{nl}^{-}$. These energies define the upper
and lower limits of an energy band, and the KS equations (1) are then solved
for energies inside the band limits to yield $X_{\epsilon nl}(r)$,
$\epsilon^{-}_{nl}\leq\epsilon\leq\epsilon^{+}_{nl}$. After discretizing the
energy bands, the density $n(r)$ is equal to

$\displaystyle
n(r)=2\sum_{k}^{N_{k}}w_{k}\sum_{nl}(2l+1)f_{knl}(\epsilon_{knl},\mu,\tau)|X_{knl}(r)|^{2}\,,$
(6) $\displaystyle w_{k}=\frac{8}{\pi(N_{k}-1)^{2}}\sqrt{k(N_{k}-1-k)}\,,$ (7)

where $N_{k}$ is the number of points used in the discretization of each
energy band. The derivation of Eqs. (6) and (7) can be found in Ref. [88]. We
note that we have changed the notation of the radial KS orbitals from
$X_{\epsilon nl}(r)$ to $X_{knl}(r)$, to signify the fact that they are
represented on a discrete, rather than continuous, energy grid.

The occupation numbers $f_{knl}$ are given in the usual way according to the
Fermi–Dirac distribution,

$f_{knl}(\epsilon_{knl},\mu,T)=\frac{1}{1+e^{(\epsilon_{knl}-\mu)/\tau}}\,.$
(8)

The chemical potential $\mu$ is determined by fixing the electron number
$N_{\textrm{e}}=4\pi\int_{0}^{R_{\textrm{VS}}}\textrm{d}rr^{2}n(r)$ to be
equal to a pre-determined value (in this paper, $N_{\textrm{e}}=Z$ in all
cases). This is an iterative procedure because the electron density, and hence
electron number, depends on the occupations $f_{knl}$, which themselves depend
on the chemical potential.

All AA calculations were performed using the open-source code atoMEC [89, 90].
Ref. [90] gives a broad overview of the code, along with a more detailed
description of how to set up and run SCF calculations. Since that publication,
the stress-tensor and virial methods for the pressure (see next section) have
been added to atoMEC, and can be found in v1.2.0 onwards. We note that the
following libraries are used extensively by atoMEC: NumPy [91], SciPy [92],
LIBXC [93], mendeleev [94], and joblib [95].

##  III Pressure in average-atom models

The total pressure in AA models is equal to the sum of the ionic and
electronic components. Typically, the ionic pressure is just given by the
ideal gas pressure,

$P_{\textrm{ion}}=\frac{n_{M}RT}{V}\,,$ (9)

where $n_{M}$ is the number of moles, $R$ the ideal gas constant, and $V$ the
volume of the gas. Clearly, this assumption is more accurate for lower
material densities and higher temperatures.

There are a number of methods to compute electronic pressure in AA models. In
this paper, we consider four different methods, presented below. We note from
the start that a number of papers have explored and derived links between
these different formulae [96, 97, 33, 98], and that the stress-tensor method
and the virial method should, in principle, be equivalent to the functional
derivative method, which is the established thermodynamic expression for the
pressure 222The ideal approximation, as implied by the name, is a known
approximation, and cannot be derived from the functional derivative of the
free energy. However, AA models rely on a range of different approximations,
and thus an established relationship for one model is unlikely to hold for a
different one. Moreover, the aim of this paper is to take a practical approach
to EOS calculations, i.e. one that is informed by comparisons to first-
principles datasets. Therefore, in this work we present the methods as we use
them, but we do not attempt to draw theoretical links between them.

###  III.1 Functional derivative of the free energy

The first method involves taking the functional derivative of the free energy
with respect to the volume,

$P^{\textrm{fd}}=-\partialderivative{F}{V}\Bigg{|}_{T}\,,$ (10)

which is a well-established way to compute the pressure in thermodynamics.
However, this method is used less frequently in AA calculations, because it is
only thermodynamically consistent if no distinction is made between bound and
unbound electrons, i.e. the same equations are solved for all orbitals,
regardless of their energy. This is true in the AA model we use, but it is not
the case for AA models in general. For example, in Ref. [26], in which the
bound and unbound orbitals are treated differently, unusual results were
observed for the pressure with this method.

In DFT, the free energy is given by

$F[n]=E[n]-TS[n],$ (11)

where $S[n]$ is the entropy and $E[n]$ is the internal energy functional. In
KS-DFT (and thus our AA model), these terms are given by 333We note that, in
Eq. (12), we have written the xc-free energy $F_{\textrm{xc}}[n]$ as part of
the internal energy. However, in principle, it also contains an entropic
contribution, because the entropy $S[n]$ is approximated by the non-
interacting entropy functional in KS-DFT. We direct readers to Ref. [26] for a
more detailed discussion of the xc energy term in ground-state and finite-
temperature KS-DFT.

$\displaystyle E[n]$
$\displaystyle=T_{\textrm{s}}[n]+E_{\textrm{en}}[n]+E_{\textrm{Ha}}[n]+F_{\textrm{xc}}[n]$
(12) $\displaystyle S[n]$
$\displaystyle=-\sum_{k}w_{k}\sum_{l,n}(2l+1)\big{[}f_{knl}\log(f_{knl})+(1-f_{knl})(\log(1-f_{knl})\big{]}\,.$
(13)

In the above, $T_{\textrm{s}}[n]$ denotes the KS kinetic energy,
$E_{\textrm{en}}[n]$ the electron-nuclear attraction energy,
$E_{\textrm{Ha}}[n]$ the Hartree energy and $F_{\textrm{xc}}[n]$ the xc free
energy. These terms are given by

$\displaystyle T_{\textrm{s}}[n]$
$\displaystyle=-2\pi\sum_{k}w_{k}\sum_{l,n}(2l+1)f_{knl}\int_{0}^{R_{\textrm{VS}}}\differential{r}r^{2}X_{knl}(r)\derivative[2]{X_{knl}(r)missing}{r}\,,$
(14) $\displaystyle E_{\textrm{en}}[n]$ $\displaystyle=-4\pi
Z\int_{0}^{R_{\textrm{VS}}}\differential{r}rn(r)\,,$ (15) $\displaystyle
E_{\textrm{Ha}}[n]$
$\displaystyle=\frac{1}{2}(4\pi)^{2}\int_{0}^{R_{\textrm{VS}}}\differential{r}r^{2}n(r)\int_{0}^{R_{\textrm{VS}}}\differential{x}\frac{n(x)x^{2}}{\max(r,x)}\,,$
(16) $\displaystyle F_{\textrm{xc}}[n]$
$\displaystyle=4\pi\int_{0}^{R_{\textrm{VS}}}\differential{r}r^{2}f_{\textrm{xc}}[n](r)n(r)\,,$
(17)

where $f_{\textrm{xc}}[n](r)$ is the xc energy density, which for LDA is the
xc energy per electron of a uniform electron gas of density $n$.

Besides the issue with thermodynamic consistency, a more minor issue with this
method is that the functional derivative in Eq. (10) is evaluated using finite
differences, which requires two separate SCF calculations at different
volumes.

###  III.2 Stress-tensor

We next consider the stress-tensor method for calculating the electronic
pressure. This method was applied earlier using only the radial component of
the stress-tensor [97, 101, 102], but more recently the full expression for
the stress-tensor was used in an AA model [33], where the AA results were seen
to be in good agreement with DFT-MD calculations and experimental data. In a
follow-up paper [103], excellent agreement was found between an AA model and
DFT-MD simulations using this method, for compressed Carbon at temperature 100
eV.

In spherical co-ordinates, the diagonal components of the stress-tensor matrix
are given by

$\displaystyle T_{rr}=P^{\textrm{st}}_{rr}$
$\displaystyle=\frac{1}{2}\sum_{k}w_{k}\sum_{n,l}2(2l+1)f_{knl}\left[\left(\derivative{X_{knl}(r)missing}{r}\right)^{2}-X_{knl}(r)\derivative[2]{X_{knl}(r)missing}{r}\right]_{r=R_{\textrm{VS}}}$
(18)
$\displaystyle=\frac{1}{2}\sum_{k}w_{k}\sum_{n,l}2(2l+1)f_{knl}\Bigg{[}\left(\derivative{X_{knl}(r)missing}{r}\right)^{2}+\frac{2}{r}X_{knl}(r)\derivative{X_{knl}(r)missing}{r}$
$\displaystyle\hskip
120.00018pt+\left(2[\epsilon_{knl}-v_{\textrm{s}}(r)]-\frac{l(l+1)}{r^{2}}\right)X_{knl}^{2}(r)\Bigg{]}_{r=R_{\textrm{VS}}}$
(19) $\displaystyle T_{\theta\theta}=T_{\phi\phi}$
$\displaystyle=\frac{1}{2}\sum_{k}w_{k}\sum_{n,l}2(2l+1)f_{knl}\left[\frac{l^{2}+l+1}{r^{2}}X_{knl}^{2}(r)-\frac{1}{2r^{3}}\derivative{[rX_{knl}(r)]^{2}}{r}\right]_{r=R_{\textrm{VS}}}$
(20)
$\displaystyle=\frac{1}{2}\sum_{k}w_{k}\sum_{n,l}2(2l+1)f_{knl}\left[\frac{l(l+1)}{r^{2}}X_{knl}^{2}(r)-\frac{1}{r}X_{knl}(r)\derivative{X_{knl}(r)missing}{r}\right]_{r=R_{\textrm{VS}}}\,.$
(21)

In Ref. [97] and Ref. [98], the radial component $T_{rr}$ is used to define
the total electronic pressure, i.e. $P^{\textrm{st}}_{rr}=T_{rr}$. On the
other hand, in Ref. [33], the (average of) the trace is taken, which leads to
the expression

$\displaystyle P^{\textrm{st}}_{\textrm{tr}}$
$\displaystyle=\frac{1}{3}(T_{rr}+T_{\theta\theta}+T_{\phi\phi})$ (22)
$\displaystyle=\frac{1}{6}\sum_{k}w_{k}\sum_{n,l}2(2l+1)f_{knl}\left\\{\left(\derivative{X_{knl}(r)missing}{r}\right)^{2}+\left[2\left(\epsilon_{knl}-v_{\textrm{s}}(r)\right)+\frac{l(l+1)}{r^{2}}\right]X_{knl}^{2}(r)\right\\}\Bigg{|}_{r=R_{\textrm{VS}}}\,.$
(23)

Eq. (23) is essentially the same expression as that presented in Eqs. (17) and
(18) of Ref. [33], but there are a few small differences. Firstly, we use the
band-structure model, which introduces the weightings $w_{k}$, and we make no
distinction between bound and free electrons, hence there is just a single
term for the pressure. Additionally, we don’t enforce a boundary condition on
the potential, i.e. $v_{\textrm{s}}(R_{\textrm{VS}})\neq 0$ in our model; this
introduces the extra potential term in Eqs. (19) and (23). Finally, the
spatial wave-functions are also defined slightly differently: the conversion
between $P_{nl}(r)$ in Ref. [33] and $X_{knl}(r)$ in our work is given by
(noting that the orbitals $X_{knl}(r)$ contain an extra index $k$ because of
the band-structure model used),

$P_{nl}(r)\rightarrow\sqrt{4\pi}rX_{knl}(r)\,.$ (24)

However, although the formula (23) is consistent with Ref. [33], the
individual components $P_{rr}$ and $P_{\theta\theta}$ are, in fact, different.
These differences cancel out when the trace is taken. Of course, this does not
affect the results in Ref. [33] because the individual components are not used
in that paper, but we mention it here to avoid confusion for the reader. When
we present results in Sec. V, we consider both the radial component,
$P^{\textrm{st}}_{rr}$ (19), and the trace, $P^{\textrm{st}}_{\textrm{tr}}$.

###  III.3 Virial theorem

Usually the virial theorem relates the pressure to the kinetic energy $T$ and
potential energy $U$ as follows [104, 105, 106]:

$P^{\textrm{vir}}_{T}=\frac{2T+U}{3V}\,,$ (25)

where the potential energy $U$ denotes the sum of all the interaction energies
in the system. In DFT, the above expression becomes [107]

$P^{\textrm{vir}}_{T}=\frac{2T_{\textrm{s}}+E_{\textrm{en}}+E_{\textrm{ha}}+W_{\textrm{xc}}}{3V}\,.$
(26)

The form of $W_{\textrm{xc}}$ depends on the type of xc-functional used. For
the LDA xc-functional used in this paper, it is given by [107]

$W_{\textrm{xc}}^{\textrm{LDA}}=-3\left[F_{\textrm{xc}}^{\textrm{LDA}}+4\pi\int_{0}^{R_{\textrm{VS}}}r^{2}n(r)v_{\textrm{xc}}^{\textrm{LDA}}(r)\right]\,.$
(27)

However, it is noted, e.g. in Refs. [97] and [98], that the expressions (25)
and (26) are only valid in the case of an infinite system
($R_{\textrm{VS}}\to\infty$), or if _all_ the wave-functions obey one of the
boundary conditions (3) or (4). Neither of these assumptions is true for the
AA model we use because the boundary conditions (3) and (4) only determine the
band-structure limits, and we solve for all wave-functions between these
limits. Instead, Refs. 97 and 98 propose the following form for the virial
pressure,

$P^{\textrm{vir}}_{K_{12}}=\frac{K_{1}+K_{2}+U}{3V}\,.$ (28)

In Eq. (28), the first term $K_{1}$ is equal to the kinetic energy defined in
Eq. (14), $K_{1}=T_{\textrm{s}}$. The second term $K_{2}$ is given by

$K_{2}=2\pi\sum_{k}w_{k}\sum_{l,n}(2l+1)f_{knl}\int_{0}^{R_{\textrm{VS}}}\differential{r}r^{2}\left\\{\left|\derivative{X_{knl}(r)missing}{r}\right|^{2}+\frac{l(l+1)}{2r^{2}}X_{knl}^{2}(r)\right\\}\,.$
(29)

It is straightforward to see that $K_{1}=K_{2}=T_{\textrm{s}}$ for the case
when either of the conditions (3) or (4) holds. In fact, $K_{2}$ is a well-
known alternative expression for the kinetic energy in KS-DFT [108, 109], and
is used, for example, to calculate the electron localization function [110].

In Refs. 98 and 97, it was shown that the expression (28) for the virial
pressure, and the radial component of the stress-tensor (19), are in principle
equivalent. However, as discussed, AA models are based on various different
assumptions which may mean that an established relationship in one model does
not hold for another. For example, in the proof presented in Ref. 97, there is
no explicit exchange-correlation term included in the virial formula.

In this work, we calculate the pressure with both forms of the virial
expression $P^{\textrm{vir}}_{T}$ (25) and $P^{\textrm{vir}}_{K_{12}}$ (28)
and benchmark the results against first-principles simulations.

###  III.4 Ideal approximation

The final method we use to compute the pressure in the AA model is based on
the assumption that the electron density at the sphere boundary is completely
free [111]. We therefore call this approach the ideal approximation, and the
electron pressure is calculated as [26]

$P_{\textrm{e}}^{\textrm{id}}=\frac{2^{3/2}}{3\pi^{2}}\int_{v_{\textrm{s}}(R_{\textrm{VS}})}^{\infty}\differential{\epsilon}\frac{\epsilon^{3/2}}{1+e^{(\epsilon-\mu)/\tau}}\,.$
(30)

The lower bound of the integral is the value of the KS potential at the atomic
radius, $v_{\textrm{s}}(R_{\textrm{VS}})$, instead of zero. This is because
the KS potential is not zero at the boundary, and electrons are considered
“free” if their energy exceeds the KS potential at the boundary.

In Ref. [33], it is shown that the above expression for the electron pressure
is consistent with both the trace of the stress-tensor (23) and the radial
component (19), under the assumption of free electron density. Of course, it
is unclear how well the assumption of free electron density at the sphere
boundary holds, particularly considering the known difficulties of defining
whether electrons are free or bound in the WDM regime [88]. However, since the
AA model is already built on numerous approximations, the accuracy of the
above expression will be assessed in the results section.

##  IV Neural network methodology

Name (formal) | Name (informal) | Definition |
---|---|---|---
$\displaystyle\epsilon_{1}(y_{1},y_{2})$ | (M)APE1 | $\displaystyle 100\times\Bigg{|}\frac{y_{1}-y_{2}}{y_{1}}\Bigg{|}$ |  (31)
$\displaystyle\epsilon_{2}(y_{1},y_{2})$ | (M)ALE2 | $\displaystyle|\log_{10}y_{1}-\log_{10}y_{2}|$ |  (32)
$\displaystyle\epsilon_{3}(y_{1},y_{2})$ | ad(M)ALE3 | $\displaystyle|\log_{10}y_{1}(\log_{10}y_{1}-\log_{10}y_{2})|$ |  (33)
$\displaystyle f_{20}$ | $f_{20}$ | $\displaystyle 100\times\frac{\mathcal{N}(\epsilon_{1}\leq 20)}{\mathcal{N}(\epsilon_{1})}$ |  (34)
$\displaystyle f_{5}$ | $f_{5}$ | $\displaystyle 100\times\frac{\mathcal{N}(\epsilon_{1}\leq 5)}{\mathcal{N}(\epsilon_{1})}$ |  (35)
Table 1: The error metrics that are used in this work. We note that the first
three metrics, $\epsilon_{i}(y_{1},y_{2})$, $i=1..3$, denote an error
measurement between two data points $y_{1}$ and $y_{2}$. If the mean (M) over
a set of data points is computed, this gives an average error. The last two
error metrics, $f_{20}$ and $f_{5}$, which represent the fraction of points
with under $20\%$ and $5\%$ errors respectively, are always aggregate
measures.
1 (Mean) absolute percentage error
2 (Mean) absolute log error
3 Adjusted (mean) absolute log error

There are various flavours of neural network, such as convolutional neural
networks [112, 113, 114], recurrent neural networks [115], and generative
adversarial networks [116]. In our work, we use a standard multi-layer
perceptron (MLP) feed-forward neural network [117]. As mentioned in the
introduction, we have trained two completely separate neural networks. The
first network is trained using only fundamental physical quantities — i.e.,
the temperature, material density, and atomic number $Z$ — together with the
ideal gas pressure, Eq. (9), as input features. The second network is trained
using various outputs of an AA calculation, chosen based on physical
intuition, in addition to the aforementioned physical quantities, as input
features. Aside from this difference, the networks were trained in an
identical manner.

We follow a nested cross-validation workflow to train and evaluate our
networks [118]. In this procedure, there is an inner cross-validation loop
during which feature selection and hyper-parameter optimization are performed;
and an outer loop, in which the generalization error (the performance of the
network on unseen data) is evaluated. Nested cross-validation is a
computationally expensive procedure, however it is recommended for relatively
small datasets (such as the one used in this paper), because it avoids
overfitting to the training sets (bias) and yields a robust estimate for the
generalization error.

In this paper, we work with several error evaluation metrics. These are used
not only in the neural network training and evaluation, but also in evaluation
of the base AA results. As shall become apparent in Sec. V, it is difficult to
compare the performance between different models using just a single error
metric, which is why we consider several in this work. In Table 1, we present
the error metrics used. The (mean) absolute percentage error (MAPE) (1) and
(mean) absolute log error (MALE) (1) are widely-used expressions, well-suited
to this work because the range of output values (pressures) spans many orders
of magnitude (making, for example, the mean-squared error a useless metric).
The adjusted (mean) absolute log error (adMALE) (1) is similar to the MALE but
gives greater weighting to higher pressures, which (as we shall later see) is
a helpful property for network training. Finally, the $f_{20}$ (1) and $f_{5}$
(1) scores are useful metrics because they are not sensitive to the presence
of large outliers.

Below, we outline the workflow we use to train and test our neural networks.
As discussed, the same procedure is used for the networks trained with and
without AA features, with the only difference being the set of features that
can be used in the training. In Secs. IV.1 and IV.2, we discuss in more detail
the feature selection and hyper-parameter optimization procedures. Note that,
below and throughout this paper, we use the notation $P_{\textrm{ref}}$ and
$Y^{0}$ to denote the reference (ground-truth) pressure values 444The reason
we use two different notations is because $P_{\textrm{ref}}$ is used to denote
the reference pressure in general, and is used later when evaluating the raw
AA results, as well as the neural network results. $Y^{0}$ is only used in the
context of the neural network training procedure, and is used to denote a
specific reference pressure for the given stage / subset of the training
workflow. Furthermore, as will be discussed in Section IV.1, $Y^{0}$ is a
transformation of the original pressure via a scaling relation., i.e. the
first-principles datasets of Ref. 28 (used for network training and
evaluation) and Ref. 29 (used for network evaluation only).

  1. 1.

Assign class-labels $C_{i},\ i=1..10$ based on the magnitude of the reference
pressure $P_{\textrm{ref}}$, i.e. the lowest 10% values of $P_{\textrm{ref}}$
are labelled $C_{1}$, and so on.

  2. 2.

Using a five-fold stratified cross-validation (SCV) approach, split the
dataset randomly (preserving equal ratios of the class labels in each
splitting) into training and test sets, $S^{\textrm{tr}}_{i}$ and
$S^{\textrm{te}}_{i},\ i=1..5$. Put the test sets aside.

  3. 3.

For each training set $S^{\textrm{tr}}_{i}$, randomly select a set of features
$\\{f_{m}\\}$

  4. 4.

For each set of features $\\{f_{m}\\}$:

    1. i

Pick a set of hyper-parameters $\\{h_{k}\\}$.

    2. ii

Split the training set $S^{\textrm{tr}}_{i}$ using five-fold SCV into new
training and validation sets, $\tilde{S}^{\textrm{tr}}_{ij}$ and
$\tilde{S}^{\textrm{te}}_{ij},\ j=1..5$.

    3. iii

Train a network for each training set $\tilde{S}^{\textrm{tr}}_{ij}$, with the
chosen hyper-parameters and features. Let us denote these networks
$g_{ijkm}[S]$, where $S$ is the data (or dataset) on which they are evaluated.

    4. iv

Compute the adMALE $\langle\epsilon_{3,ikm}\rangle$ over $j=1..5$ for the
given hyper-parameters and features:

$\langle\epsilon_{3,ikm}\rangle=\frac{1}{N^{\textrm{tr}}_{i}}\sum_{j=1}^{5}\sum_{n=1}^{N^{\textrm{te}}_{ij}}\epsilon_{3}(Y^{0}_{ijn},g_{ijkm}[\tilde{S}^{\textrm{te}}_{ijn}])$
(36)

where $N^{\textrm{tr}}_{i}$ is the total number of samples in the training set
$\tilde{S}^{\textrm{tr}}_{i}$,
$N^{\textrm{tr}}_{i}=\sum_{j=1}^{5}N^{\textrm{te}}_{ij}$, and $n$ denotes a
sample within that set. We note that all the reference pressures $Y^{0}$ are
positive, and because we use the rectified linear activation function in the
neural networks, all the predictions are also strictly positive. Therefore
there are no problems taking the logarithms.

    5. v

Pick a new set of $\\{h_{k}\\}$. Repeat steps ii-vi 10 times, and finally save
the minimum error metric and the associated hyper-parameters:

$\langle\bar{\epsilon}^{0}_{3,im}\rangle=\min_{\\{h_{k}\\}}\langle\epsilon_{3,ikm}\rangle$
(37)

  5. 5.

Repeat steps 3-4 20 times, saving the error metric each time and features used
for each repetition. At the end, select the _three_ lowest error metrics and
save their associated features:

$\displaystyle\langle\bar{\epsilon}_{3,i}^{0}\rangle$
$\displaystyle=\min_{\\{f_{m}\\}}\langle\bar{\epsilon}^{0}_{3,im}\rangle$ (38)
$\displaystyle\langle\bar{\epsilon}_{3,i}^{1}\rangle$
$\displaystyle=(\min_{\\{f_{m}\\}}+1)\langle\bar{\epsilon}^{0}_{3,im}\rangle$
(39) $\displaystyle\langle\bar{\epsilon}_{3,i}^{2}\rangle$
$\displaystyle=(\min_{\\{f_{m}\\}}+2)\langle\bar{\epsilon}^{0}_{3,im}\rangle\,,$
(40)

where $(\min+1)$ and $(\min+2)$ denote the second and third lowest errors
respectively. We shall comment on this step in more detail below, because it
is not a standard part of a nested cross-validation (CV) procedure.

  6. 6.

Following steps 3-5, we end up with three optimal sets of features
$\\{f_{m}^{l}\\}$ and hyper-parameters $\\{h_{k}^{l}\\}$, $l=1,2,3$, for the
$i$-th iteration of the outer CV loop. Using these optimal sets, train three
networks with the _full_ training set $S_{i}^{\textrm{tr}}$. We denote these
networks $g_{il}[S]$.

  7. 7.

Using the networks trained in the above step, make predictions on the
corresponding _test_ dataset $S_{i}^{\textrm{te}}$. The final prediction,
which we denote $\bar{g}_{i}[S]$, is given by a simple average of the three
individual predictions,

$\bar{g}_{i}[S]=\frac{1}{3}\sum_{l=0}^{2}g_{il}[S].$ (41)

With these predictions, evaluate various error metrics
$\langle\bar{\epsilon}_{j,i}\rangle$ on the corresponding _test_ dataset
$S_{i}^{\textrm{te}}$. Finally, take the average over all the test datasets to
compute the final error metrics. For example, the MAPE (1), after averaging
over all test sets, is defined as

$\langle\bar{\epsilon}_{1}\rangle=\frac{1}{N_{\textrm{tot}}}\sum_{i=1}^{5}\sum_{n=1}^{N_{i}^{\textrm{te}}}\epsilon_{1}\left(Y^{0,\textrm{te}}_{in},\bar{g}_{i}[S^{\textrm{te}}_{in}]]\right),$
(42)

where $N_{\textrm{tot}}=\sum_{i=1}^{5}N_{i}^{\textrm{te}}$ is the total number
of samples in the full dataset.

  8. 8.

The above steps tell us the _generalization_ error as opposed to yielding a
specific model. To train the final model, repeat steps 4-6, but instead of
looping over separate training sets $S_{i}^{\textrm{tr}}$, use the _full_
dataset.

We emphasize that steps 1 and 2 are carried out once and only once, before any
network training or feature selection begins. This ensures there is no
contamination of training and test data.

As promised above, we now elaborate a bit on step 5, in which the best three
combinations of features and hyper-parameters were chosen. Typically, only a
single set of features and hyper-parameters (that with the lowest error) would
be carried forward in this stage. However, in order to reduce the bias (the
likelihood of over-fitting to a specific model), we take forward the three
best networks and then later use them to make an average prediction. This idea
is borrowed from the idea of ensemble models in machine-learning [120], where
multiple models are trained and the outputs from those models are used to
build a new model; however, our case is a very simple example because we just
do a basic averaging in the final step, without any optimization on the inner
validations sets.

###  IV.1 Feature selection and scaling

Initially, a set of possible features was chosen for the two networks. For the
AA-free network, we considered only some fundamental physical quantities:
namely, the material density and temperature, the volume of the Voronoi
sphere, the atomic number, and the ideal gas pressure for the ions (9). For
the AA network, various outputs from the AA model were considered as features:
for example, the electronic pressure computed using the four different methods
described in Section III, the value of the density at the Voronoi sphere
boundary, the mean ionization state, and so on. These features were selected
from a physically intuitive viewpoint, because there is reason to expect a
correlation between the feature and the reference pressure.

An initial thinning of the features was performed by hand. Using the first
training set following the outer SCV splitting, $S_{1}^{\textrm{tr}}$, we
first computed the correlation of each feature and the reference pressure,
using the Kendall–Tau correlation measure [121]. The results (in descending
order of correlation) are shown in Table 2. Following this, the three features
with the lowest correlations, namely the atomic number $Z$, the chemical
potential $\mu$, and the free energy $F[n]$, were dropped (there were
additional reasons for dropping $F[n]$ and $\mu$, which are discussed in the
paragraph following Fig. 1.

Next, the correlations between the different features were checked. The
potential $v_{\textrm{s}}(R_{\textrm{VS}})$ and electron density
$n(R_{\textrm{VS}})$ at the cell boundary were found to be fully anti-
correlated ($C=-1$); in fact, this is expected due to the Hohenberg–Kohn
theorem in DFT, which establishes a one-to-one mapping between density and
potential [122]. As a result, $v_{\textrm{s}}(R_{\textrm{VS}})$ was excluded
from the list of potential features. Although the correlations between the
different AA electron pressures were high ($C\geq 0.9$), we kept all of them
as potential features as their correlation with the reference pressure was
also high.

Feature | AA network | AA-free network | Correlation to target
---|---|---|---
$P^{\textrm{st}}_{rr}$ | ✓ | ✗ | 0.928
$P^{\textrm{id}}$ | ✓ | ✗ | 0.927
$P^{\textrm{st}}_{\textrm{tr}}$ | ✓ | ✗ | 0.925
$P^{\textrm{fd}}$ | ✓ | ✗ | 0.919
$P^{\textrm{vir}}_{T}$ | ✓ | ✗ | 0.916
$P^{\textrm{vir}}_{K_{12}}$ | ✓ | ✗ | 0.900
$\derivative*{n(r)}{r}|_{R_{\textrm{VS}}}$ | ✓ | ✗ | 0.844
$P_{\textrm{ion}}$ | ✓∗ | ✓ | 0.834
$\derivative*{v_{\textrm{s}}(r)}{r}|_{R_{\textrm{VS}}}$ | ✓ | ✗ | 0.749
${n(R_{\textrm{VS}})}$ | ✓ | ✗ | 0.687
${v_{\textrm{s}}(R_{\textrm{VS}})}$ | ✗ | ✗ | 0.687
$T$ | ✓ | ✓ | 0.661
$\rho_{\textrm{m}}/Z$ | ✓ | ✓ | 0.540
$V$ | ✓ | ✓ | 0.531
$Z^{*}$ | ✓ | ✗ | 0.502
$\rho_{\textrm{m}}$ | ✓ | ✓ | 0.490
$F[n]$ | ✗ | ✗ | 0.434
$\mu$ | ✗ | ✗ | 0.238
$Z$ | ✗ | ✗ | 0.116
Table 2: Initial features considered for the AA and AA-free neural networks,
and the Kendall-Tau correlation figures between the features and the reference
pressure. We also show which features are used for both networks during the
training.
∗ Summed together with the electron pressures in the network with AA features

In Table 2, we show which features remained and which were discarded following
this manual thinning of the feature space, for the AA and AA-free networks.
For the remaining features, we plotted the relationship between the reference
pressure and feature, to get a better insight into their relationship. We
observed that, for every feature considered, a much better correlation could
be observed between the feature and reference pressure when logarithms of both
were taken. A couple of examples of the features plotted against the pressure,
with and without taking logarithms, can be found in Fig. 1.  We note that, in
this table, $Z^{*}$ denotes the mean ionization state (MIS), which is defined
here as

$Z^{*}=2\sum_{k}w_{k}\sum_{n,l}(2l+1)f_{knl}(\mu,\tau)\Theta\big{(}\epsilon_{knl}-v_{\textrm{s}}(R_{\textrm{VS}})\big{)}\,,$
(43)

where $\Theta\big{(}\epsilon_{knl}-v_{\textrm{s}}(R_{\textrm{VS}})\big{)}$ is
the Heaviside step function. In other words, the MIS is defined as the number
of electrons with energy above the potential energy at the sphere boundary.

We note that the AA pressures sometimes have a negative sign, but in these
limited instances the absolute value of the pressure was taken instead. All
the other features considered had a consistent sign; this further motivated
the dropping of the chemical potential and free energy as features, since they
spanned a wide spectrum of positive and negative signs for the dataset used,
which complicated the use of logarithms.

Figure 1: Features (output from AA calculation) vs reference pressure. Top
row: linear scale. Bottom row: logarithmic scale. Left column: derivative of
KS potential at sphere boundary. Right column: MIS (43). This plot indicates
the potential benefits of using a logarithmic scale for both the features and
the reference pressure.

Furthermore, some preliminary networks were trained using the training set
$S^{\textrm{tr}}_{1}$. These networks were not trained using the full workflow
described in the previous sub-section. Rather, a few simple networks —
consisting of one or two hidden layers and between 10-30 neurons per layer —
were trained (also varying a few other hyper-parameters), just to get an
initial idea of whether the logarithmic scaling improved the network
performance. The outcome of this investigation clearly indicated that taking
the logarithm of the reference pressure and all features yields significantly
lower errors for all the different network architectures. Consequently, this
logarithmic scaling was adopted for all subsequent training of networks.

As described in the previous sub-section, in step 3 of the workflow a random
subset of features is chosen from the full feature space. In total, there were
14 available features for the AA network (after the manual reduction performed
above) and just 5 for the AA-free network. For the AA network, we randomly
selected between 5 and 13 features in step 3. As discussed, following hyper-
parameter optimization for this set of features, the process was repeated; in
total, it was repeated 20 times for the AA network. For the AA-free network,
since the number of features is limited, we randomly chose between 3-5
features, and again repeated the process 20 times. Of course, there is only
one subset of 5 features in this case, but repeating the hyper-parameter
optimization with this same set of features is still beneficial, because the
hyper-parameter search will change every time.

In Fig. 2 we plot the error metrics for the feature-sampling procedure
described above, for the AA network. Note that each of the dark blue points in
this graph denotes the error $\langle\bar{\epsilon}_{3,im}^{0}\rangle$ (37),
i.e. the error for that particular feature subset, following the inner CV
search for the best hyper-parameters. This plot is an amalgamation of all the
outer loop training sets $S_{i}^{\textrm{tr}}$. The dark orange line tracks
the average (mean) error over the number of features selected for the sample.
This information is also summarized in Table 3. We observe that, in general,
more features leads to better performance; however, the best performing models
actual have similar errors, regardless of the number of features used.

Following the logarithmic scaling described above, we scaled both the features
and the reference pressure for each of the training sets
$S_{\textrm{i}}^{\textrm{tr}}$. This was done via the standard feature scaling
formula,

$\tilde{f}_{nm}=\frac{f_{nm}-u_{m}}{s_{m}},$ (44)

where $f_{nm}$ and $\tilde{f}_{nm}$ are the original and scaled features for
the $n$-th data point, and $u_{m}$ and $s_{m}$ are respectively the mean and
standard deviation of the $m$-th feature. The same approach was used for the
reference pressure. Feature scaling is a standard procedure in machine
learning; on the other hand, the target is not usually scaled. However, the
preliminary networks we trained performed much better when the target was also
scaled, and therefore we adopted this approach.

Feature number | MAPE (%) | MALE
---|---|---
5 | 3.8 | 0.015
6 | 2.5 | 0.0099
7 | 2.2 | 0.0088
8 | 2.0 | 0.0085
9 | 2.0 | 0.0086
10 | 1.6 | 0.0069
11 | 1.8 | 0.0079
12 | 1.7 | 0.0071
13 | 1.6 | 0.0066
Table 3: Average mean absolute percentage error (MAPE) and mean absolute log
error (MALE) for different numbers of features, over the _training_ sets (for
the AA neural network). Figure 2: Mean absolute percentage error (MAPE) (top)
and mean absolute log error (bottom) for the AA neural networks trained with
different numbers of features. The average for each number of features is
shown with the orange line.

###  IV.2 Network architecture and hyper-parameter optimization

As mentioned earlier, we use a standard feed-forward neural network
architecture. As discussed in the previous sub-section, we trained some
initial networks using just the training set $S_{1}^{\textrm{tr}}$ to get a
general idea of what features and hyper-parameters to use in our networks.
This initial exploration indicated which hyper-parameters should be fixed, and
of the hyper-parameters that were not fixed, what range of values to consider
during the automated hyper-parameter optimization phase.

Following this initial training, we decided to use the ‘Adam’ stochastic
optimization method [123] and the rectified linear (ReLU) activation function,
$\textrm{ReLU}(x)=\max(x,0)$. Below, we summarize the hyper-parameters that we
allowed to vary during the automated search for the best hyper-parameters:

  * •

The number of hidden layers in the network, $N_{\textrm{hidden}}$

  * •

The number of neurons in the _first_ hidden layer, $N_{\textrm{neuron}}^{(1)}$

  * •

The fraction by which the number of neurons in successive layers are reduced
in size, $\alpha_{r}$; for example, the second hidden layer has
$N_{\textrm{neuron}}^{(2)}=\alpha_{r}N_{\textrm{neuron}}^{(1)}$ neurons

  * •

The number of epochs over which training is performed, $N_{\textrm{epoch}}$

  * •

The size of the mini-batch used for training, $N_{\textrm{batch}}$

  * •

The learning rate $l_{r}$

All remaining hyper-parameters were fixed to their default arguments in the
TensorFlow library [124].

With the hyper-parameters listed above, the Optuna library [125] was used to
search for the best hyper-parameters. For each hyper-parameter, a search
window was specified (based on observations from the preliminary networks that
were trained); the search windows that were used are shown in Table 4. As
mentioned, the optimal hyper-parameters were chosen as those which minimize
the adMALE (36) between the prediction and reference pressures for a given set
of features. This error metric is related to the mean-absolute log error
(MALE), however, we scale the errors with the logarithm of the reference
pressure. This scaling is done because the raw AA models perform best at high
pressures (resulting from high-temperatures and densities), and it would be
undesirable to train a neural network that performs badly under these
conditions. This scaling gives a slightly bigger weighting to higher-pressure
errors.

On the other hand, the loss function used in the network training was the
standard mean-absolute log error (MALE) between the reference and predicted
pressures,

$L_{ijkm}=\frac{1}{N_{ij}^{\textrm{tr}}}\sum_{n=1}^{N_{ij}^{\textrm{tr}}}\epsilon_{2}(Y^{0}_{ijn},g_{ijkm}[\tilde{S}_{ijn}^{\textrm{tr}}])\,.$
(45)

This was chosen (as opposed to the weighted MALE described above) because the
MAE is a regularly-used loss function in neural networks, and therefore is
implemented in an optimal way in neural network libraries. However, there is
an added bonus: by using the MALE as the loss function and the weighted MALE
as the error metric for hyper-parameter optimization and feature selection, we
actually minimize (to some extent) two important errors.

We finish this section by briefly summarizing the training and evaluation
procedures for our neural network models, and their architectures. Our neural
network models are based on a standard feed-forward (MLP) network architecture
[117]. They are shallow networks, with a maximum of two hidden layers used
during the hyper-parameter optimization. We train two kinds of networks, one
whose input features are just basic physical properties, and another which
uses additional features from the output of an AA calculation. We use a nested
CV approach [118]: this means hyper-parameter optimization and feature
selection are performed on an inner CV loop, meanwhile errors are evaluated on
an outer CV loop, which ensures no contamination of training and test data.
Several error metrics, shown in Table 1, are calculated. Borrowing from the
idea of ensembles in machine learning [120], final predictions are constructed
by taking the average predictions of the three best models.

Hyper-parameter | Min value | Max value
---|---|---
$N_{\textrm{hidden}}$ | 1 | 2
$N_{\textrm{neuron}}^{(1)}$ | 20 | 80
$\alpha_{r}$ | 0.5 | 1.0
$N_{\textrm{epoch}}$ | 1000 | 2000
$N_{\textrm{batch}}$ | 20 | 60
$l_{r}$ | $10^{-4}$ | $10^{-3}$
Table 4: Search windows for hyper-parameters

##  V Results

###  V.1 FPEOS database of Militzer _et al._ : Average-atom results

We start by presenting the results of the AA model, using the four different
approximations for the pressure described in Sec. III, compared to the FPEOS
dataset of Militzer _et al._ [28]. From this dataset, we use only single-
element data, with atomic number $1\leq Z\leq 14$. We do not use the mixtures
from this dataset, because at present our AA method is not suited for mixture
calculations. Furthermore, for a small subset of the conditions (high
temperatures and low densities), we were unable to perform AA calculations.
This is a limitation of the atoMEC code, which uses sparse matrix
diagonalization to solve the KS equations [90]; under these conditions, the
number of eigenvalues (KS orbitals) required is so large that the sparse
matrix diagonalization breaks down. In general, however, these conditions are
accessible without issues for AA codes. Nevertheless, this still leaves us
with 2181 total data points, out of a possible 2371 had we been able to
perform all calculations. In Fig. 3, we show the data points which we were
able to compute with the AA model (in green), and those which were
inaccessible (in red). As a brief aside, we note that the same set of 2181
points was used in the training of both neural networks (Section V.2). In
principle, all the data points from the FPEOS database [28] could be used to
train the network which does not use AA features; however, we choose to use
the same subset to facilitate comparison between the two networks and the raw
AA results.

Figure 3: The data-points from the FPEOS database [28] that could be
calculated with the atoMEC code (shown in green), and those that were
inaccessible (in red), in terms of the density and temperature.

In Fig. 4, we plot the absolute percentage error (APE) (1) between the AA and
FP pressures as a function of temperature,
$\epsilon_{1}(P_{\textrm{ref}},P_{\textrm{AA}})$, with the density also
highlighted using a colour scale. Note that the $y$-axis is linear up to 20%,
and after that the scale is logarithmic. The two dashed lines (and associated
shaded areas) denote errors of 5% and 20% respectively.

From this figure, it is clear that the AA model has large errors for low
temperatures ($\lessapprox 10$ eV), independent of the approximation used for
the error. The errors appear to be particularly severe when the density is
also low, at least for the finite-difference and virial methods. This could
partially be attributed to the asymmetry of the MAPE, where over-estimates can
lead to especially large errors. However, under these conditions (low
temperature and density), one would expect the electronic pressure to approach
zero, with the only contribution to the pressure coming from the ions; as mass
density increases for fixed temperature, we expect the electron (and ion)
pressures to also increase.

Figure 4: Absolute percentage errors (APE) in pressure between different AA
methods and first-principles data [28], as a function of temperature. A colour
scale is also included to indicate the effect of mass density. Note that the
$y$-axis scale is linear up to $20\%$ (with the dashed lines representing
errors of $5\%$ and $20\%$ respectively), and logarithmic above $20\%$.

However, this is not (in general) the case: consider, for example, Fig. 5. In
this figure, we plot the total pressure as a function of mass density for
Helium, at fixed temperature $T=0.043$ eV (the lowest temperature in the FPEOS
database). We see a strange behaviour in the finite-difference pressure and
both virial pressures: it is negative and does not monotonically increase, as
should be expected for a fixed temperature and increasing density. On the
other-hand, both stress-tensor methods and the ideal pressure show the correct
qualitative behaviour, with the stress-tensor results also showing good
quantitative accuracy in this example. Of course, the ideal gas expression for
the ions is also less accurate at low temperatures, particularly as the
density increases; it is difficult to say how much the error in the total
pressure is affected by errors in the ion and electron pressures, but it is
likely both are incorrect to some extent.

Figure 5: Pressure as a function of mass density, for Helium at temperature
$T=0.043$ eV, with the various AA methods compared to the first-principles
reference [28]. We observe that, for both the virial methods and the finite-
difference method, the pressure actually _decreases_ as density _increases_
for the first few data points. This is both counter-intuitive, and also does
not agree with the benchmark.

Nevertheless, for all the pressure approximations with the exception of
$P^{\textrm{st}}_{\textrm{tr}}$, the AA model improves steeply with increasing
temperature. This is a positive finding, because we would expect the AA model
to be most accurate at higher temperatures, due to the decreasing importance
of quantum effects. The trend appears to be most strong for the finite-
difference and the $P^{\textrm{vir}}_{T}$ virial methods. A somewhat puzzling
observation is the systematic errors of $\sim 10-20\%$ in the high-temperature
region for the $P^{\textrm{st}}_{\textrm{tr}}$ method. This is in contrast to
the results of Refs. [33] and [103], in which this method agreed very well
with DFT-MD simulations under similar conditions. This could be related to the
differences in the AA model we used compared to the model of Refs. [33] and
[103].

Next, in Fig. 6, we plot the AA pressures against the reference FPEOS
pressure, on a logarithmic scale. In the lower panels of this plot, we plot
the absolute log errors (ALE) (1) between the AA and reference pressures,
$\epsilon_{2}(P_{\textrm{ref}},P_{\textrm{AA}})$. In this figure, it is clear
that, for high pressures (high material density and/or temperatures), the AA
results are in close agreement with the FPEOS reference for almost all the
methods, but in particular for the virial and finite-difference methods. As
was observed in Fig. 4 for the APE, upon close observation, the stress-tensor
results $P^{\textrm{st}}_{\textrm{tr}}$ appear to be systematically lower than
the reference results at high pressures ($P\gtrapprox 10^{5}$). However, this
behaviour is more noticeable when the APE is plotted (compared to the ALE);
this indicates the importance of considering different error metrics when
comparing pressures across a large range of temperatures and densities.

Figure 6: AA pressures (various methods) vs first-principles data [28],
plotted on a $\log-\log$ scale. The absolute log errors (ALEs) (1) are also
shown in the smaller sub-plots. In general, we observe the tendency of the AA
model to perform better at higher-pressures.

In Table 5, we collect some error metrics for all the different AA methods.
From this table, we notice that the $P^{\textrm{fd}}$, $P^{\textrm{st}}_{rr}$,
$P^{\textrm{vir}}_{T}$ and $P^{\textrm{id}}$ methods yield quite similar
results for the majority of the error metrics. The
$P^{\textrm{st}}_{\textrm{tr}}$ and $P^{\textrm{vir}}_{K_{12}}$ show greater
differences; in particular, they both yield a far smaller fraction of results
within a $5\%$ error (by an order of magnitude in the case of
$P^{\textrm{st}}_{\textrm{tr}}$) and also have larger adMALE values,
demonstrating bigger errors at higher pressures. This table further
demonstrates the importance of considering a range of error metrics; for
example, $P^{\textrm{st}}_{\textrm{tr}}$ has a lower MAPE than
$P^{\textrm{fd}}$ and $P^{\textrm{vir}}_{T}$, but performs worse (and
typically far worse) in all the other error metrics.

| $P^{\textrm{fd}}$ | $P^{\textrm{st}}_{rr}$ | $P^{\textrm{st}}_{\textrm{tr}}$ | $P^{\textrm{vir}}_{T}$ | $P^{\textrm{vir}}_{K_{12}}$ | $P^{\textrm{id}}$
---|---|---|---|---|---|---
MAPE | 33 | 27 | 29 | 34 | 48 | 26
MALE | 8.7 | 8.7 | 11 | 8.8 | 13 | 10
adMALE | 21 | 23 | 39 | 21 | 35 | 26
$f_{20}$ | 75 | 73 | 55 | 76 | 72 | 72
$f_{5}$ | 64 | 55 | 5.4 | 60 | 24 | 55
Table 5: Aggregate error metrics for the AA methods, compared to the FPEOS
dataset [28]. See Table 1 for the definitions of the errors. Note that here,
and in all subsequent tables, the MALE and adMALE results are multiplied by a
factor of 100.

To finish our analysis of the pressure results with the AA model, consider
Table 6. This table is analagous to Table 5, however we now only consider data
with temperature above $10$ eV. This quantifies the big improvements at higher
temperatures seen in Fig. 4, and implicitly in Fig. 6: in particular, the
$P^{\textrm{fd}}$, $P^{\textrm{st}}_{rr}$, $P^{\textrm{vir}}_{T}$, and
$P^{\textrm{id}}$ methods all have MAPEs of $\sim 4\%$, which may be
considered within the bounds of experimental uncertainties. Above $10$ eV, all
these methods display very similar results, and significant improvements to
when all temperatures were considered. The $P^{\textrm{vir}}_{K_{12}}$ method
also shows improvements across all metrics, but peforms somewhat worse overall
than the aforementioned methods. However, as expected from Fig. 4, the
$P^{\textrm{st}}_{\textrm{tr}}$ method does not show any improvement (in fact
the results are arguably worse) when temperatures below 10 eV are eliminated.

| $P^{\textrm{fd}}$ | $P^{\textrm{st}}_{rr}$ | $P^{\textrm{st}}_{\textrm{tr}}$ | $P^{\textrm{vir}}_{T}$ | $P^{\textrm{vir}}_{K_{12}}$ | $P^{\textrm{id}}$
---|---|---|---|---|---|---
MAPE | 2.6 | 4.3 | 18 | 3.1 | 8.0 | 3.9
MALE | 1.1 | 1.8 | 8.7 | 1.4 | 3.8 | 1.6
adMALE | 5.2 | 7.4 | 43 | 6.3 | 17 | 6.9
$f_{20}$ | 98 | 97 | 60 | 99 | 94 | 79
$f_{5}$ | 87 | 78 | 1.6 | 83 | 46 | 96
Table 6: Aggregate error metrics for the AA methods, compared to the FPEOS
dataset[28]. The difference between this table and Table 5 is that here we
only consider data with temperature above 10 eV, whereas in Table 5 the whole
dataset was considered.

###  V.2 FPEOS database of Militzer _et al._ : Neural network results

Now, we consider results for the neural networks trained by the procedure
described in Section. IV. As mentioned, we trained two networks, one with
features from the output of AA calculations, and one without. The results for
the neural network are the aggregate of the results over the outer five test
sets of the nested CV procedure.

First, consider Fig. 7, in which we plot the APE of the AA neural network
(left panel) and the AA-free neural network (right panel). What is immediately
clear from this figure, when compared to Fig. 4, is that both neural networks
show a dramatic improvement relative to the raw AA results, at least for low-
to-moderate temperatures ($T<100\ \textrm{eV}$). In particular, there are now
only a handful of results with errors of over 20%, and indeed it seems that
most lie within $5\%$, even at low temperatures.

Figure 7: Absolute percentage errors (APE) for the neural network models
evaluated on the FPEOS dataset [28] as a function of temperature, with a
colour bar to indicate the dependence on mass density. Left: model trained
with AA output as features. Right: model trained without any AA outputs.

Fig. 8 paints a similar picture; on a logarithmic scale, it is challenging to
discern any real differences between the neural-network pressures and the
reference pressure. In Fig. 9, we directly compare the the APE (left) as a
function of temperature, and the ALE (right) as a function of the reference
pressure, between the two neural network models and the $P^{\textrm{fd}}$ raw
AA result. We observe clear improvements in the neural network models relative
to the raw AA results in this figure, at least for $T<100\ \textrm{eV}$ and
$P_{\textrm{ref}}<10^{6}\ \textrm{GPa}$.

Figure 8: Pressures computed from neural network models against the benchmark
pressures from Ref. [28]. The lower panels show the logarithmic errors (1).
Left: network trained with AA output as features. Right: network trained
without any AA outputs. Figure 9: Comparison of errors between a raw AA model
($P^{\textrm{fd}}$), and neural networks trained with and without AA output,
evaluated on the FPEOS dataset [28]. Left: absolute percentage errors (1).
Right: absolute log errors (1).

In Table. 7, we show the aggregate error metrics for the neural network models
and the $P^{\textrm{fd}}$ AA results. These numbers further demonstrate the
stark improvements with the neural network models, as the error metrics are
uniformly (and significantly) lower than for the raw AA results. For example,
respectively only $0.4\%$ and $1.1\%$ of results now show an error of more
than $20\%$ compared to the reference. Interestingly, there is not much
difference between the networks trained with and without the AA outputs,
although the AA network does have slightly lower average errors for all the
metrics considered.

| $P_{\textrm{nn}}^{\textrm{AA}}$ | $P_{\textrm{nn}}^{\textrm{no-AA}}$ | $P^{\textrm{fd}}$
---|---|---|---
MAPE | 1.9 (0.2) | 2.6 (0.5) | 33 (2)
MALE | 0.77 (0.1) | 1.0 (0.2) | 8.7 (0.5)
adMALE | 2.5 (0.4) | 3.1 (0.7) | 21 (0.7)
$f_{20}$ | 99.6 (0.3) | 98.9 (0.5) | 75 (1)
$f_{5}$ | 94 (2) | 90 (3) | 63 (2)
Table 7: Aggregate error metrics for the AA and AA-free neural networks, with
AA results ($P^{\textrm{fd}}$) also shown for comparison. The numbers in
brackets for the MAPE, MALE and adMALE are the standard deviations of these
quantities across the five outer cross-validation folds.

Of course, it would come as a great surprise if the neural networks did not
show a (significant) improvement to the raw AA results, since a neural network
(at least in this case) is essentially a polynomial fitting to the reference
data. Interpolation methods for EOS data are already a well-established tool
in the WDM and plasma-physics communities. In Appendix A, we evaluate the
performance of direct interpolation methods provided in Ref. 28 using a
similar cross-validation approach to the one used to test our neural network
models. One of the conclusions we draw from Appendix A is that our method has
greater capability to extrapolate to densities and temperatures outside of the
training range. Another benefit of our method is that it is more global, in
that it does not fit to a specific element but instead to the whole space of
elements considered. Results supporting this claim will be shown in Sec. V.3.
Our work has further novelty because, to the best of our knowledge, output
from AA calculations has not previously been used to supplement interpolations
of first-principles (or alternative sources of high-fidelity) data. Our work
is similar in spirit to Ref. [78], but it is distinguished by the inclusion of
the AA outputs, besides various other differences such as the reference data,
the neural network architecture, and the training and evaluation framework.

###  V.3 FPEOS Be data of Ding and Hu: Average-atom and neural network results

Next, we consider the FPEOS database from Ref. [29], which we henceforth
denote as the FP-Be dataset. In this case, all the data comes from DFT-MD
calculations, with KS-DFT-MD used for temperatures up to 250,000 K ($21.5\
\textrm{eV}$), and OF-DFT-MD used for the calculations with temperatures above
that. They observed that, at the temperature transition point, the results
between the two methods were consistent to within $1\%$. Beryllium is the only
element in this dataset. We note that, in the results we shall present, we use
only the raw data and not any extrapolated or interpolated values; and like
for the FPEOS data, there were a small number of high-temperature and low-
density calculations which we are currently not able to perform with atoMEC.

Figure 10: Distribution of temperatures (left), mass densities (middle) and
reference pressures (right) in the FPEOS [28] and FP-Be [29] datasets.

We emphasize that this dataset is completely separate from the FPEOS data used
to train the neural network models, and furthermore, that Beryllium does not
feature in the FPEOS dataset. Therefore this tests the ability of the neural
network models to extrapolate to elements on which they were not trained,
instead of the more straightforward interpolation considered in the previous
section. Of course, Beryllium has an atomic number of 4, which lies inside the
range $1\leq Z\leq 14$ spanned by the FPEOS dataset. We also note that the FP-
Be data has a similar temperature and density distribution as the FPEOS
dataset, as shown in Fig. 10. In this figure, we see the temperature profiles
are almost identical, and while the FP-Be data is slightly more skewed towards
higher densities (and hence higher pressures), the difference is quite small.

In Fig. 11, we compare the APEs across all the different AA models. As with
Fig. 4, we see that all methods with the exception of
$P^{\textrm{st}}_{\textrm{tr}}$ tend to approach the reference pressure as
temperature increases; and we again observe very large errors for
$T\lessapprox 10$ eV. This observation is especially stark for the
$P^{\textrm{st}}_{rr}$ method, in which the errors drop sharply from $\sim
100\%$ to $<5\%$ almost immediately as temperature increases beyond 10 eV.

Figure 11: Absolute percentage errors (APE) in pressure between different AA
methods and Be first-principles data [29], as a function of temperature, with
the colour scale indicating the dependence on mass density.

Next, consider Fig. 12, in which we plot the APE for the two neural network
models. Both models show a clear improvement relative to the raw AA results in
the low-to-mid temperature range; this is particularly true for the model
trained with AA features, which has very few results with an error $>20\%$. On
the other hand, under high temperatures at which the AA models perform
strongly, there is no evidence of improvement for the neural networks.
However, it is worth noting that the AA results are themselves highly accurate
at these temperatures, so it would be difficult for the neural networks to
demonstrate much improvement.

Figure 12: Absolute percentage errors (APE) for the neural network models
evaluated on the FP-Be dataset [29] as a function of temperature, with a
colour bar to indicate the dependence on mass density. Left: model trained
with AA output as features. Right: model trained without any AA outputs.

The above observations are further borne out when considering Figs. 13 and 14.
We can see, by direct comparison of the errors between the neural network and
AA model (which was chosen as the virial method $P^{\textrm{vir}}_{T}$ in this
case), that the AA result has larger errors and lower temperatures and
pressures. However, it rapidly converges to the reference. Both neural
networks perform relatively strongly across the full range of conditions
spanned, albeit with evidence of slightly larger errors at the highest
temperatures and pressures. In these figures, the network trained with AA data
seems markedly more accurate than the raw AA results or the network trained
without AA data.

Figure 13: Pressures computed from neural network models and the
$P^{\textrm{vir}}_{T}$ AA method, against the benchmark pressures from Ref.
[29]. The lower panels show the logarithmic errors (1). Left: network trained
with AA output as features. Right: network trained without any AA outputs.
Figure 14: Comparison of errors between an AA model ($P^{\textrm{vir}}_{T}$),
and neural networks trained with and without AA output, evaluated on the FP-Be
dataset [29]. Left: absolute percentage errors (1). Right: absolute log errors
(1).

Now consider Table 8, in which we compare error metrics for the whole set of
AA pressures and also the neural network predictions. Of the AA methods,
$P^{\textrm{fd}}$ and $P^{\textrm{vir}}_{T}$ seem to overall perform best, and
also perform better than they had done on the FPEOS data. One cause could be
the slightly different mass density distributions between the two datasets;
the MAPE (and to a lesser extent the MALE and adMALE) can be easily influenced
by large errors. The $f_{5}$ and $f_{20}$ scores, which are not unduly
influenced by outliers, demonstrate the most consistency between the two
datasets. The predictions from the neural network trained with AA data show
significant improvement in all error metrics relative to the raw AA results;
meanwhile, the predictions from the neural network trained without AA data
also show (more moderate) improvement compared to the raw AA results in all
metrics aside from the $f_{5}$ score. However, the improvement of the neural
network models relative to the raw AA results for the FP-Be dataset is not as
strong as the equivalent improvement for the FPEOS dataset. We can see this,
for example, by comparing Table 5 and Table 8. There is a bigger difference in
the various error metrics between the $P^{\textrm{fd}}$ AA pressure in Table 5
and the neural network pressures, than there is in Table 8. This is an
indication that the neural networks can interpolate to elements to which they
were not trained on, but that the performance will drop slightly as a result.

| $P^{\textrm{fd}}$ | $P^{\textrm{st}}_{rr}$ | $P^{\textrm{st}}_{\textrm{tr}}$ | $P^{\textrm{vir}}_{T}$ | $P^{\textrm{vir}}_{K_{12}}$ | $P^{\textrm{id}}$ | $P_{\textrm{nn}}^{\textrm{AA}}$ | $P_{\textrm{nn}}^{\textrm{no-AA}}$
---|---|---|---|---|---|---|---|---
MAPE | 18 | 68 | 44 | 15 | 11 | 36 | 4.5 | 7.0
MALE | 5.4 | 13 | 13 | 4.8 | 5.2 | 9.9 | 1.9 | 3.1
adMALE | 15 | 37 | 46 | 13 | 19 | 31 | 6.6 | 8.9
$f_{20}$ | 79 | 67 | 46 | 83 | 88 | 65 | 98.8 | 95
$f_{5}$ | 59 | 59 | 5.5 | 57 | 38 | 55 | 70 | 47
Table 8: Aggregate error metrics for the AA methods, compared to the FP-Be
dataset [29]. See Table 1 for definitions of the error metrics.

As a final analysis, consider Table 9, in which we compare error metrics for
the AA methods against the neural networks, but this time including only
temperatures above 10 eV in our analysis. We see that the AA results are very
comparable to Table 6; here, the $P^{\textrm{fd}}$ and $P^{\textrm{vir}}_{T}$
methods show the best performance when all metrics are considered, with
$P^{\textrm{id}}$ and $P^{\textrm{st}}_{rr}$ slightly behind. Again, the
$P^{\textrm{st}}_{\textrm{tr}}$ method is a clear outlier. Interestingly, the
neural network results do not seem to show improvement relative to the best AA
methods in Table 9. This suggests that, when considering elements on which
training has not been performed, and temperatures above 10 eV, the base AA
model (with the $P^{\textrm{vir}}_{T}$ or $P^{\textrm{fd}}$ method) may be
more suitable than either of the neural network models. Of course, this
conclusion does not apply when considering the full temperature range, or for
elements present in the training data.

| $P^{\textrm{fd}}$ | $P^{\textrm{st}}_{rr}$ | $P^{\textrm{st}}_{\textrm{tr}}$ | $P^{\textrm{vir}}_{T}$ | $P^{\textrm{vir}}_{K_{12}}$ | $P^{\textrm{id}}$ | $P_{\textrm{nn}}^{\textrm{AA}}$ | $P_{\textrm{nn}}^{\textrm{no-AA}}$
---|---|---|---|---|---|---|---|---
MAPE | 2.8 | 5.2 | 18 | 2.9 | 8.7 | 6.2 | 3.4 | 5.1
MALE | 1.2 | 2.1 | 8.7 | 1.3 | 4.1 | 2.4 | 1.5 | 2.4
adMALE | 4.3 | 7.7 | 4.1 | 4.7 | 18 | 11 | 6.1 | 7.1
$f_{20}$ | 100 | 93 | 57 | 100 | 95 | 89 | 100 | 96
$f_{5}$ | 80 | 82 | 6.0 | 76 | 47 | 76 | 76 | 63
Table 9: Aggregate error metrics for the AA methods, compared to the FP-Be
dataset [29]. The difference between this table and Table 8 is that here we
only consider data with temperature above 10 eV, whereas in Table 8, the whole
FP-Be dataset was considered. See Table. 1 for definitions of error metrics.

##  VI Conclusions

In this paper, we sought to develop methods to accurately and efficiently
compute equation-of-state (EOS) data, specifically the pressure, in the warm
dense matter (WDM) regime. First-principles methods to compute accurate EOS
data under WDM conditions — namely, density-functional theory molecular-
dynamics (DFT-MD) and path integral Monte Carlo (PIMC) — are well known.
However, these methods are far too expensive to generate on-the-fly data for
hydrodynamic simulations. AA calculations, on the other hand, often take only
a few seconds [86], and are fast enough to run in-line with hydrodynamic
simulations [126, 127].

In the first part of the paper, we first compared the performance of an AA
model to the FPEOS dataset of Militzer _et al._ [28]. Within this model, we
considered various methods for calculating pressure, which were detailed in
Section III. We observed that, whilst all the methods have significant errors
at relatively low temperatures ($\lessapprox 10$ eV), they almost unanimously
converge to the reference results as temperature increases. For example, as
seen in Table 6, if only temperatures above $10$ eV are included, the mean
absolute percentage error (MAPE) is just 2.6%, and 98% of the results lie
within 20% of the reference. These errors are of the same order of magnitude
as could be expected from experimental measurements. A similar performance was
observed when the AA pressures were compared with the FP-Be database of Ding
and Hu [29].

However, the primary goal of this paper was to investigate the application of
neural networks to EOS calculations: to this end, we developed two neural
networks, using the procedure described in Section IV. The difference between
the two networks was that one used output data from AA calculations — see
Table 2 — as input features to the network. The other just used fundamental
physical data such as the temperature and mass density. Both these models
demonstrated the ability to accurately interpolate the pressure, with MAPEs of
1.9% (AA network) and 2.6% (AA free network) across the whole temperature
range respectively.

There are, of course, classical interpolation schemes which can be used to
rapidly compute EOS data based on the first-principles datasets used in this
paper. However, our method has several advantages relative to these schemes.
Most importantly, our neural network models are global (rather than being
confined to a specific element). This means that, firstly, our methods can be
applied to elements outside of the original training set. We tested our pre-
trained networks on the FP-Be dataset [29] to evaluate their performance on an
element not in the training set. The network trained with AA features
demonstrated significant improvement relative to the raw AA results in the
low-to-mid temperature range, and maintained a roughly equal performance with
the AA model at temperatures above 10 eV, when the AA model is itself highly
accurate. The second advantage of our methods being global is that they can be
extrapolated beyond the boundaries of the interpolation tables. This is
discussed in detail in Appendix A, in which we compare direct interpolation of
the FPEOS tables with our neural network models. Furthermore, though beyond
the scope of this paper, we note that neural network models hold the promise
of extensibility to non-local thermal equilibrium (NLTE) modelling, as
outlined in Appendix A.

Overall, the results in this paper indicate that — with a good choice of
method for calculating the pressure — AA models generally show accurate
agreement with DFT-MD and PIMC benchmarks at moderate-to-high temperatures
($\gtrapprox 10$ eV). In fact, it is at these temperatures that AA models are
most useful, because DFT-MD calculations are more computationally feasible at
low temperatures.

Moreover, we have developed neural network models which, especially if trained
with output data from AA calculations, are highly effective interpolation
tools for first-principles EOS data. There is even evidence that these models,
in particular the model trained with AA outputs, can be applied successfully
to elements on which they were not trained: for example, the networks
demonstrate significant improvements relative to raw AA results at low
temperatures ($<10$ eV). This offers a promising solution for a global EOS
method that is both accurate and efficient throughout the WDM regime.

## Data Availability

The AA results and final trained networks can be downloaded from this link:
rodare.hzdr.de/record/2289. The code required to run the AA calculations, and
train and test the neural networks, can be accessed here: github.com/atomec-
project/neuralEOS. The FPEOS database of Militzer _et al._ can be downloaded
from the Supplementary Material of Ref. 28. The FPEOS Be data from Ding and Hu
can be requested from the authors of Ref. [29], as stated in that paper.

## Acknowledgements

We thank S. Hu for providing us with the FP-Be dataset [29]. This work was
partially supported by the Center for Advanced Systems Understanding (CASUS),
which is financed by Germany’s Federal Ministry of Education and Research
(BMBF) and by the Saxon state government out of the State budget approved by
the Saxon State Parliament. We declare that there are no conflicts of
interest.

##  Appendix A Comparison of neural network models with FPEOS direct
interpolation

In hydrodynamic simulations, EOS data is required on-the-fly, and therefore
interpolation of pre-computed EOS tables is a common practice. In the FPEOS
database used in this paper [28], classical interpolation routines are
provided with the data. In this appendix, we compare the accuracy and
computational load of these interpolation routines with our neural network
models.

First, let us consider the computational load of the interpolation routines
relative to our neural network approaches. After the neural networks have been
trained — which is a one-time procedure as described in Section IV —
performing inference with the networks, i.e. applying the trained networks to
new data, is effectively instantaneous, just like standard interpolation
routines.

However, the network which uses AA data as input features has an extra
computational load, namely the computation of the AA data. As mentioned in the
conclusions of this paper, there exist AA codes which are fast enough to be
run inline with hydrodynamic simulations [86, 126, 127]. On the other hand,
the atoMEC code we have used is not (in general) fast enough to generate on-
the-fly data for hydrodynamic simulations. Fig. 15 shows the AA SCF
calculation times for all the calculations in the FPEOS dataset. The median
time is 16 seconds, with the 25% and 75% quartiles being 10 and 31 seconds
respectively. Some calculations, however, are much more expensive (taking over
600 seconds, i.e. 10 minutes). The cost increases when more orbitals are
required, i.e. when the density is lower or the temperature higher (the cost
is also dependent on convergence parameters such as the grid size).

Firstly, we note that development of the atoMEC code only began two years ago,
and so far was not focussed on performance. We imagine that the code could be
sufficiently optimized so that it performs as well as the codes mentioned in
Refs. [86, 127]. Secondly — and more importantly — it is in fact not necessary
to run an AA code inline to generate the input data for the neural networks.
In fact, we are developing an approach to interpolate AA data across the full
thermodynamic spectrum that can be covered by the atoMEC code, which will
grant instantaneous access to AA data.

As an aside, we note that this approach has the additional advantage that it
can be extended to non-local thermal equilibrium (NLTE) conditions. Solving
the collisional-radiative equations requires spectral information as input
data, which is most accurately obtained from a fully-quantum method such as an
AA code. Under NLTE conditions, multiple configurations must be accounted for,
and it is therefore not possible to run even the fastest AA codes inline; the
number of configurations must be truncated, severely limiting the accuracy.
Furthermore, by using a neural network to interpolate the AA data, rather than
a classical interpolation scheme, the complexity of the system is reduced and
we expect to be able to scale the method to NLTE.

Figure 15: Average-atom SCF calculation times. The left panel is a box plot
(showing the median time, interquartile ranges, limits (without outliers), and
outliers. The right panel shows the time taken per calculation, as a function
of temperature and density.

Next, we compare the accuracy of direct interpolation of the FPEOS data with
the AA neural network model. Firstly, we recall that the FPEOS data can only
be interpolated on a species-by-species basis, and therefore the interpolation
routines can only be used for elements present in the dataset. A major
advantage of our neural network approach is that it was succesfully applied,
for example, to Beryllium, which cannot be done with the FPEOS routines.

Figure 16: Comparison of predicted vs reference pressures for direct
interpolation of FPEOS data. Left panels: Predicted vs reference pressures on
a log-log scale. Right panels: Absolute percentage errors between predicted
and reference pressures, as a function of density and temperature. Top panels:
Evaluation is performed on the full test sets. Bottom panels: Evaluation is
performed on the restricted test sets.

To assess the accuracy of direct FPEOS interpolation, we used a cross-
validaiton testing strategy. That is to say, we randomly split the data for
each element into 5 non-overlapping “training” and test sets. The data was
fitted (i.e. the interpolation functions constructed) using the “training”
sets, and then the accuracy was assessed on the test sets. We note that, with
this random splitting, it is possible to end up with test data that falls
outside of the training range, i.e. with densities or temperatures in the test
set that are fall outside of the bounds in the training set. This is no
different from the neural network training, but in that case the bounds are
much wider for each element, because the training is performed over all the
elements present in the dataset. To fairly assess the accuracy of the FPEOS
direct interpolation, we introduce additional test sets. These are restricted
subsets of the original test data sets, in which data that falls outside of
the training temperature and density range is removed from the test set. We
test on both the full and restricted sets, and use the notation
$P_{\textrm{pred}}^{\textrm{f}}$ or $P_{\textrm{FPEOS}}^{\textrm{f}}$, and
$P_{\textrm{pred}}^{\textrm{r}}$ or $P_{\textrm{FPEOS}}^{\textrm{r}}$, to
denote respectively when the full or restricted sets are used.

| $P_{\textrm{FPEOS}}^{\textrm{r}}$ | $P_{\textrm{FPEOS}}^{\textrm{f}}$ | $P_{\textrm{nn}}^{\textrm{AA}}$
---|---|---|---
MAPE | 1.6 (0.7) | 1300 (2000) | 1.9 (0.2)
MALE | 0.64 (0.2) | 1.9 (0.7) | 0.77 (0.1)
adMALE | 2.2 (0.3) | 3.5 (0.5) | 2.5 (0.4)
$f_{20}$ | 99.7 (0.3) | 98 (0.3) | 99.6 (0.3)
$f_{5}$ | 95 (2) | 93 (1) | 94 (2)
Table 10: Comparison of error metrics between FPEOS direct interpolation with
AA neural network model. The MAPE values for $P_{\textrm{FPEOS}}^{\textrm{f}}$
are somewhat meaningless, due to the presence of large outliers.

To evaluate the performance of the FPEOS networks, we first consider Fig. 16.
In the top panels of Fig. 16, the FPEOS interpolation is evaluated on the full
test sets; meanwhile, in the bottom panel, it is evaluated on the restricted
test sets, which exclude data falling outside the bounds of the training data.
We see that the bottom panels have far fewer large errors than the top panels;
in particular, the $P_{\textrm{pred}}^{\textrm{f}}$ predictions have some
extremely large errors of over $100\%$. However, in general, the vast majority
of the data is closely aligned with the reference (as expected).

Next, we consider Table 10, in which we compare error metrics for
$P_{\textrm{FPEOS}}^{\textrm{r}}$, $P_{\textrm{FPEOS}}^{\textrm{f}}$, and the
AA neural network model $P_{\textrm{nn}}^{\textrm{AA}}$. It is clear from this
table that the AA model is slightly less accurate than the FPEOS interpolation
on the restricted test sets, $P_{\textrm{FPEOS}}^{\textrm{r}}$. However, it is
more accurate than the FPEOS interpolation on the full datasets,
$P_{\textrm{FPEOS}}^{\textrm{f}}$, and in particular has fewer large outliers.
Overall, we can conclude that direct interpolation of the FPEOS dataset is
highly accurate; however, even when going slightly beyond the limits of the
training data, the interpolation quickly becomes out of control. This is a key
advantage of our method, which does not rely on an element-by-element
interpolation, and therefore can be used to make predictions over a much wider
range of temperatures and densities.

## References

  * [1] M Koenig, A Benuzzi-Mounaix, A Ravasio, T Vinci, N Ozaki, S Lepape, D Batani, G Huser, T Hall, D Hicks, A MacKinnon, P Patel, H S Park, T Boehly, M Borghesi, S Kar, and L Romagnani.  Progress in the study of warm dense matter.  Plasma Physics and Controlled Fusion, 47(12B):B441–B449, nov 2005\.
  * [2] Basic Research Needs for High Energy Density Laboratory Physics.  U.S. DOE, 2009.
  * [3] E. Brown, M.A. Morales, C. Pierleoni, and D. Ceperley.  Quantum Monte Carlo Techniques and Applications for Warm Dense Matter, volume 96 of Lecture Notes in Computational Science and Engineering, pages 123–149.  Springer International Publishing, 2014.
  * [4] Katerina Falk.  Experimental methods for warm dense matter research.  High Power Laser Science and Engineering, 6:e59, 2018.
  * [5] M. Bonitz, T. Dornheim, Zh. A. Moldabekov, S. Zhang, P. Hamann, H. Kählert, A. Filinov, K. Ramakrishna, and J. Vorberger.  Ab initio simulation of warm dense matter.  Physics of Plasmas, 27(4):042710, 2020.
  * [6] Tobias Dornheim, Simon Groth, and Michael Bonitz.  The uniform electron gas at warm dense matter conditions.  Physics Reports, 744:1–86, 2018.
  * [7] A. L. Kritcher, T. Döppner, C. Fortmann, T. Ma, O. L. Landen, R. Wallace, and S. H. Glenzer.  In-flight measurements of capsule shell adiabats in laser-driven implosions.  Phys. Rev. Lett., 107:015002, Jul 2011.
  * [8] John D. Lindl, Peter Amendt, Richard L. Berger, S. Gail Glendinning, Siegfried H. Glenzer, Steven W. Haan, Robert L. Kauffman, Otto L. Landen, and Laurence J. Suter.  The physics basis for ignition using indirect-drive targets on the National Ignition Facility.  Phys. Plasmas, 11(2):339–491, 2004.
  * [9] R. Betti and O. A. Hurricane.  Inertial-confinement fusion with lasers.  Nature Physics, 12(5):435–448, May 2016.
  * [10] R. S. Craxton, K. S. Anderson, T. R. Boehly, V. N. Goncharov, D. R. Harding, J. P. Knauer, R. L. McCrory, P. W. McKenty, D. D. Meyerhofer, J. F. Myatt, A. J. Schmitt, J. D. Sethian, R. W. Short, S. Skupsky, W. Theobald, W. L. Kruer, K. Tanaka, R. Betti, T. J. B. Collins, J. A. Delettrez, S. X. Hu, J. A. Marozas, A. V. Maximov, D. T. Michel, P. B. Radha, S. P. Regan, T. C. Sangster, W. Seka, A. A. Solodov, J. M. Soures, C. Stoeckl, and J. D. Zuegel.  Direct-drive inertial confinement fusion: A review.  Physics of Plasmas, 22(11):110501, 2015.
  * [11] W. B. Hubbard, T. Guillot, J. I. Lunine, A. Burrows, D. Saumon, M. S. Marley, and R. S. Freedman.  Liquid metallic hydrogen and the structure of brown dwarfs and giant planets.  Physics of Plasmas, 4(5):2011–2015, 1997.
  * [12] Andrea L. Kritcher, Damian C. Swift, Tilo Döppner, Benjamin Bachmann, Lorin X. Benedict, Gilbert W. Collins, Jonathan L. DuBois, Fred Elsner, Gilles Fontaine, Jim A. Gaffney, Sebastien Hamel, Amy Lazicki, Walter R. Johnson, Natalie Kostinski, Dominik Kraus, Michael J. MacDonald, Brian Maddox, Madison E. Martin, Paul Neumayer, Abbas Nikroo, Joseph Nilsen, Bruce A. Remington, Didier Saumon, Phillip A. Sterne, Wendi Sweet, Alfredo A. Correa, Heather D. Whitley, Roger W. Falcone, and Siegfried H. Glenzer.  A measurement of the equation of state of carbon envelopes of white dwarfs.  Nature, 584(7819):51–54, Aug 2020.
  * [13] J Daligault and S Gupta.  Electron–ion scattering in dense multi-component plasmas: Application to the outer crust of an accreting neutron star.  The Astrophysical Journal, 703(1):994, 2009.
  * [14] J. J. Fortney, S. H. Glenzer, M. Koenig, B. Militzer, D. Saumon, and D. Valencia.  Frontiers of the physics of dense plasmas and planetary interiors: Experiments, theory, and applications.  Physics of Plasmas, 16(4):041003, 2009.
  * [15] Jeffrey H. Nguyen and Neil C. Holmes.  Melting of iron at the physical conditions of the earth’s core.  Nature, 427(6972):339–342, 2004.
  * [16] Bruce A Remington.  High energy density laboratory astrophysics.  Plasma Physics and Controlled Fusion, 47(5A):A191, apr 2005.
  * [17] J. Clérouin, P. Noiret, P. Blottiau, V. Recoules, B. Siberchicot, P. Renaudin, C. Blancard, G. Faussurier, B. Holst, and C. E. Starrett.  A database for equations of state and resistivities measurements in the warm dense matter regime.  Physics of Plasmas, 19(8):082702, 2012.
  * [18] Siegfried H. Glenzer and Ronald Redmer.  X-ray Thomson scattering in high energy density plasmas.  Rev. Mod. Phys., 81:1625–1663, Dec 2009.
  * [19] O. S. Humphries, R. S. Marjoribanks, Q. Y. van den Berg, E. C. Galtier, M. F. Kasim, H. J. Lee, A. J. F. Miscampbell, B. Nagler, R. Royle, J. S. Wark, and S. M. Vinko.  Probing the electronic structure of warm dense nickel via resonant inelastic x-ray scattering.  Phys. Rev. Lett., 125:195001, Nov 2020.
  * [20] O. Ciricosta, S. M. Vinko, H.-K. Chung, B.-I. Cho, C. R. D. Brown, T. Burian, J. Chalupský, K. Engelhorn, R. W. Falcone, C. Graves, V. Hájková, A. Higginbotham, L. Juha, J. Krzywinski, H. J. Lee, M. Messerschmidt, C. D. Murphy, Y. Ping, D. S. Rackstraw, A. Scherz, W. Schlotter, S. Toleikis, J. J. Turner, L. Vysin, T. Wang, B. Wu, U. Zastrau, D. Zhu, R. W. Lee, P. Heimann, B. Nagler, and J. S. Wark.  Direct measurements of the ionization potential depression in a dense plasma.  Phys. Rev. Lett., 109:065002, Aug 2012.
  * [21] D. Kraus, J. Vorberger, J. Helfrich, D. O. Gericke, B. Bachmann, V. Bagnoud, B. Barbrel, A. Blažević, D. C. Carroll, W. Cayzac, T. Döppner, L. B. Fletcher, A. Frank, S. Frydrych, E. J. Gamboa, M. Gauthier, S. Göde, E. Granados, G. Gregori, N. J. Hartley, B. Kettle, H. J. Lee, B. Nagler, P. Neumayer, M. M. Notley, A. Ortner, A. Otten, A. Ravasio, D. Riley, F. Roth, G. Schaumann, D. Schumacher, W. Schumaker, K. Siegenthaler, C. Spindloe, F. Wagner, K. Wünsch, S. H. Glenzer, M. Roth, and R. W. Falcone.  The complex ion structure of warm dense carbon measured by spectrally resolved x-ray scattering.  Physics of Plasmas, 22(5):056307, 2015.
  * [22] E. I. Moses, R. N. Boyd, B. A. Remington, C. J. Keane, and R. Al-Ayat.  The National Ignition Facility: Ushering in a new age for high energy density science.  Phys. Plasmas, 16(4):041006, 2009.
  * [23] Christoph Bostedt, Sébastien Boutet, David M. Fritz, Zhirong Huang, Hae Ja Lee, Henrik T. Lemke, Aymeric Robert, William F. Schlotter, Joshua J. Turner, and Garth J. Williams.  Linac coherent light source: The first five years.  Rev. Mod. Phys., 88:015007, Mar 2016.
  * [24] David Pile.  First light from sacla.  Nature Photonics, 5(8):456–457, Aug 2011.
  * [25] Thomas Tschentscher, Christian Bressler, Jan Grünert, Anders Madsen, Adrian P. Mancuso, Michael Meyer, Andreas Scherz, Harald Sinn, and Ulf Zastrau.  Photon beam transport and scientific instruments at the european xfel.  Applied Sciences, 7(6), 2017.
  * [26] T. J. Callow, S. B. Hansen, E. Kraisler, and A. Cangi.  First-principles derivation and properties of density-functional average-atom models.  Phys. Rev. Research, 4:023055, Apr 2022.
  * [27] J.A. Gaffney, S.X. Hu, P. Arnault, A. Becker, L.X. Benedict, T.R. Boehly, P.M. Celliers, D.M. Ceperley, O. Čertík, J. Clérouin, G.W. Collins, L.A. Collins, J.-F. Danel, N. Desbiens, M.W.C. Dharma-wardana, Y.H. Ding, A. Fernandez-Pañella, M.C. Gregor, P.E. Grabowski, S. Hamel, S.B. Hansen, L. Harbour, X.T. He, D.D. Johnson, W. Kang, V.V. Karasiev, L. Kazandjian, M.D. Knudson, T. Ogitsu, C. Pierleoni, R. Piron, R. Redmer, G. Robert, D. Saumon, A. Shamp, T. Sjostrom, A.V. Smirnov, C.E. Starrett, P.A. Sterne, A. Wardlow, H.D. Whitley, B. Wilson, P. Zhang, and E. Zurek.  A review of equation-of-state models for inertial confinement fusion materials.  High Energy Density Physics, 28:7–24, 2018.
  * [28] Burkhard Militzer, Felipe González-Cataldo, Shuai Zhang, Kevin P. Driver, and Fran çois Soubiran.  First-principles equation of state database for warm dense matter computation.  Phys. Rev. E, 103:013203, Jan 2021.
  * [29] Y. H. Ding and S. X. Hu.  First-principles equation-of-state table of beryllium based on density-functional theory calculations.  Physics of Plasmas, 24(6):062702, 2017.
  * [30] S. X. Hu, B. Militzer, V. N. Goncharov, and S. Skupsky.  First-principles equation-of-state table of deuterium for inertial confinement fusion applications.  Phys. Rev. B, 84:224109, Dec 2011.
  * [31] David A. Liberman.  Self-consistent field model for condensed matter.  Phys. Rev. B, 20:4981–4989, Dec 1979.
  * [32] C. E. Starrett, J. Clérouin, V. Recoules, J. D. Kress, L. A. Collins, and D. E. Hanson.  Average atom transport properties for pure and mixed species in the hot and warm dense matter regimes.  Physics of Plasmas, 19(10):102709, 2012.
  * [33] Gérald Faussurier and Christophe Blancard.  Pressure in warm and hot dense matter using the average-atom model.  Phys. Rev. E, 99:053201, May 2019.
  * [34] C. E. Starrett, J. Daligault, and D. Saumon.  Pseudoatom molecular dynamics.  Phys. Rev. E, 91:013104, Jan 2015.
  * [35] R. Piron and T. Blenski.  Variational-average-atom-in-quantum-plasmas (vaaqp) code and virial theorem: Equation-of-state and shock-hugoniot calculations for warm dense al, fe, cu, and pb.  Phys. Rev. E, 83:026403, Feb 2011.
  * [36] James D. McHardy.  An introduction to the theory and use of sesame equations of state.  12 2018.
  * [37] D Kraus, B Bachmann, B Barbrel, R W Falcone, L B Fletcher, S Frydrych, E J Gamboa, M Gauthier, D O Gericke, S H Glenzer, S Göde, E Granados, N J Hartley, J Helfrich, H J Lee, B Nagler, A Ravasio, W Schumaker, J Vorberger, and T Döppner.  Characterizing the ionization potential depression in dense carbon plasmas with high-precision spectrally resolved x-ray scattering.  Plasma Physics and Controlled Fusion, 61(1):014015, nov 2018.
  * [38] Tobias Dornheim, Maximilian Böhme, Dominik Kraus, Tilo Döppner, Thomas R. Preston, Zhandos A. Moldabekov, and Jan Vorberger.  Accurate temperature diagnostics for matter under extreme conditions.  Nature Communications, 13:7911, 2022.
  * [39] Zhandos A. Moldabekov, Tobias Dornheim, and Attila Cangi.  Thermal excitation signals in the inhomogeneous warm dense electron gas.  Scientific Reports, 12(1), jan 2022.
  * [40] Tobias Dornheim, Maximilian P. Böhme, David A. Chapman, Dominik Kraus, Thomas R. Preston, Zhandos A. Moldabekov, Niclas Schlünzen, Attila Cangi, Tilo Döppner, and Jan Vorberger.  Imaginary-time correlation function thermometry: A new, high-accuracy and model-free temperature analysis technique for x-ray Thomson scattering data.  Physics of Plasmas, 30(4), 04 2023.  042707\.
  * [41] W. Kohn and L. J. Sham.  Self-consistent equations including exchange and correlation effects.  Phys. Rev., 140:A1133–A1138, Nov 1965.
  * [42] Yan Alexander Wang and Emily A. Carter.  Orbital-Free Kinetic-Energy Density Functional Theory, pages 117–184.  Springer Netherlands, Dordrecht, 2002.
  * [43] John P. Perdew and Karla Schmidt.  Jacob’s ladder of density functional approximations for the exchange-correlation energy.  AIP Conference Proceedings, 577(1):1–20, 2001.
  * [44] Aron J. Cohen, Paula Mori-Sánchez, and Weitao Yang.  Challenges for density functional theory.  Chem. Rev., 112(1):289–320, 2012.
  * [45] Neil Qiang Su and Xin Xu.  Development of new density functional approximations.  Annual Review of Physical Chemistry, 68(1):155–182, 2017.  PMID: 28226221.
  * [46] Valentin V. Karasiev, Lázaro Calderín, and S. B. Trickey.  Importance of finite-temperature exchange correlation for warm dense matter calculations.  Phys. Rev. E, 93:063207, Jun 2016.
  * [47] Kushal Ramakrishna, Tobias Dornheim, and Jan Vorberger.  Influence of finite temperature exchange-correlation effects in hydrogen.  Phys. Rev. B, 101:195129, May 2020.
  * [48] J. C. Smith, A. Pribram-Jones, and K. Burke.  Exact thermal density functional theory for a model system: Correlation components and accuracy of the zero-temperature exchange-correlation approximation.  Phys. Rev. B, 93:245131, Jun 2016.
  * [49] Brittany P. Harding, Zachary Mauri, and Aurora Pribram-Jones.  Approximate bounds and temperature dependence of adiabatic connection integrands for the uniform electron gas.  The Journal of Chemical Physics, 156(13):134104, 2022.
  * [50] Zhandos Moldabekov, Tobias Dornheim, Maximilian Böhme, Jan Vorberger, and Attila Cangi.  The relevance of electronic perturbations in the warm dense electron gas.  The Journal of Chemical Physics, 155(12):124116, 2021.
  * [51] Zhandos Moldabekov, Tobias Dornheim, Jan Vorberger, and Attila Cangi.  Benchmarking exchange-correlation functionals in the spin-polarized inhomogeneous electron gas under warm dense conditions.  Phys. Rev. B, 105:035134, Jan 2022.
  * [52] Zhandos A. Moldabekov, Mani Lokamani, Jan Vorberger, Attila Cangi, and Tobias Dornheim.  Assessing the accuracy of hybrid exchange-correlation functionals for the density response of warm dense electrons.  The Journal of Chemical Physics, 158(9), 03 2023.  094105\.
  * [53] Zhandos A. Moldabekov, Mani Lokamani, Jan Vorberger, Attila Cangi, and Tobias Dornheim.  Non-empirical mixing coefficient for hybrid xc functionals from analysis of the xc kernel.  The Journal of Physical Chemistry Letters, 14(5):1326–1333, 2023\.  PMID: 36724891.
  * [54] V.V. Karasiev and S.B. Trickey.  Issues and challenges in orbital-free density functional calculations.  Computer Physics Communications, 183(12):2519–2527, 2012.
  * [55] Valentin V. Karasiev, Debajit Chakraborty, Olga A. Shukruto, and S. B. Trickey.  Nonempirical generalized gradient approximation free-energy functional for orbital-free simulations.  Phys. Rev. B, 88:161108, Oct 2013.
  * [56] Mel Levy and Hui Ou-Yang.  Exact properties of the pauli potential for the square root of the electron density and the kinetic energy functional.  Phys. Rev. A, 38:625–629, Jul 1988.
  * [57] K. Finzel.  The first order atomic fragment approach—An orbital-free implementation of density functional theory.  The Journal of Chemical Physics, 151(2), 07 2019.  024109\.
  * [58] Eduardo V. Ludeña, Edison X. Salazar, Mauricio H. Cornejo, Darío E. Arroyo, and Valentin V. Karasiev.  The liu-parr power series expansion of the pauli kinetic energy functional with the incorporation of shell-inducing traits: Atoms.  International Journal of Quantum Chemistry, 118(14):e25601, 2018\.
  * [59] D. M. Ceperley.  Fermion nodes.  Journal of Statistical Physics, 63(5):1237–1267, Jun 1991.
  * [60] Matthias Troyer and Uwe-Jens Wiese.  Computational complexity and fundamental limitations to fermionic quantum monte carlo simulations.  Phys. Rev. Lett., 94:170201, May 2005.
  * [61] Maximilian Böhme, Zhandos A. Moldabekov, Jan Vorberger, and Tobias Dornheim.  Static electronic density response of warm dense hydrogen: Ab initio path integral monte carlo simulations.  Phys. Rev. Lett., 129:066402, Aug 2022.
  * [62] Kamal Choudhary, Brian DeCost, Chi Chen, Anubhav Jain, Francesca Tavazza, Ryan Cohn, Cheol Woo Park, Alok Choudhary, Ankit Agrawal, Simon J. L. Billinge, Elizabeth Holm, Shyue Ping Ong, and Chris Wolverton.  Recent advances and applications of deep learning methods in materials science.  npj Computational Materials, 8(1):59, 2022.
  * [63] L. Fiedler, K. Shah, M. Bussmann, and A. Cangi.  Deep dive into machine learning density functional theory for materials science and chemistry.  Phys. Rev. Mater., 6:040301, Apr 2022.
  * [64] Ryan Pederson, Bhupalee Kalita, and Kieron Burke.  Machine learning and density functional theory.  Nature Reviews Physics, 4(6):357–358, 2022.
  * [65] Mihail Bogojeski, Leslie Vogt-Maranto, Mark E. Tuckerman, Klaus-Robert Müller, and Kieron Burke.  Quantum chemical accuracy from density functional approximations via machine learning.  Nature Communications, 11(1):5223, 2020.
  * [66] James Kirkpatrick, Brendan McMorrow, David H. P. Turban, Alexander L. Gaunt, James S. Spencer, Alexander G. D. G. Matthews, Annette Obika, Louis Thiry, Meire Fortunato, David Pfau, Lara Román Castellanos, Stig Petersen, Alexander W. R. Nelson, Pushmeet Kohli, Paula Mori-Sánchez, Demis Hassabis, and Aron J. Cohen.  Pushing the frontiers of density functionals by solving the fractional electron problem.  Science, 374(6573):1385–1389, 2021.
  * [67] David Pfau, James S. Spencer, Alexander G. D. G. Matthews, and W. M. C. Foulkes.  Ab initio solution of the many-electron schrödinger equation with deep neural networks.  Phys. Rev. Res., 2:033429, Sep 2020.
  * [68] Giuseppe Carleo and Matthias Troyer.  Solving the quantum many-body problem with artificial neural networks.  Science, 355(6325):602–606, 2017.
  * [69] Peter W. Hatfield, Jim A. Gaffney, Gemma J. Anderson, Suzanne Ali, Luca Antonelli, Suzan Başeğmez du Pree, Jonathan Citrin, Marta Fajardo, Patrick Knapp, Brendan Kettle, Bogdan Kustowski, Michael J. MacDonald, Derek Mariscal, Madison E. Martin, Taisuke Nagayama, Charlotte A. J. Palmer, J. Luc Peterson, Steven Rose, J J Ruby, Carl Shneider, Matt J. V. Streeter, Will Trickey, and Ben Williams.  The data-driven future of high-energy-density physics.  Nature, 593(7859):351–361, may 2021.
  * [70] Tobias Dornheim, Attila Cangi, Kushal Ramakrishna, Maximilian Böhme, Shigenori Tanaka, and Jan Vorberger.  Effective static approximation: A fast and reliable tool for warm-dense matter theory.  Phys. Rev. Lett., 125:235001, Dec 2020.
  * [71] J. A. Ellis, L. Fiedler, G. A. Popoola, N. A. Modine, J. A. Stephens, A. P. Thompson, A. Cangi, and S. Rajamanickam.  Accelerating finite-temperature kohn-sham density functional theory with deep neural networks.  Phys. Rev. B, 104:035120, Jul 2021.
  * [72] Lenz Fiedler, Nils Hoffmann, Parvez Mohammed, Gabriel A Popoola, Tamar Yovell, Vladyslav Oles, J Austin Ellis, Sivasankaran Rajamanickam, and Attila Cangi.  Training-free hyperparameter optimization of neural networks for electronic structures in matter.  Machine Learning: Science and Technology, 3(4):045008, oct 2022\.
  * [73] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [74] Kurt Hornik.  Approximation capabilities of multilayer feedforward networks.  Neural Networks, 4(2):251–257, 1991.
  * [75] S. J. Ali, D. C. Swift, C. J. Wu, and R. G. Kraus.  Development of uncertainty-aware equation-of-state models: Application to copper.  Journal of Applied Physics, 128(18):185902, 2020.
  * [76] Jim A Gaffney, Lin Yang, and Suzanne Ali.  Constraining model uncertainty in plasma equation-of-state models with a physics-constrained gaussian process, 2022.
  * [77] Beth A. Lindquist and Ryan B. Jadrich.  Uncertainty quantification for a multi-phase carbon equation of state model.  Journal of Applied Physics, 131(15):155104, 2022.
  * [78] Katherine L. Mentzer and J. Luc Peterson.  Neural network surrogate models for equations of state.  Physics of Plasmas, 30(3):032704, 2023.
  * [79] Of course, the performance of the network is evaluated on data that was not seen during training; the point is that the training and test data is drawn from the same database.
  * [80] E. Wigner and F. Seitz.  On the constitution of metallic sodium.  Phys. Rev., 43:804–810, May 1933.
  * [81] R. P. Feynman, N. Metropolis, and E. Teller.  Equations of state of elements based on the generalized fermi-thomas theory.  Phys. Rev., 75:1561–1573, May 1949.
  * [82] M. W. C. Dharma-wardana and Fran çois Perrot.  Density-functional theory of hydrogen plasmas.  Phys. Rev. A, 26:2096–2104, Oct 1982.
  * [83] Balazs F. Rozsnyai.  Photoabsorption in hot plasmas based on the ion-sphere and ion-correlation models.  Phys. Rev. A, 43:3035–3042, Mar 1991.
  * [84] C. Blancard and G. Faussurier.  Equation of state and transport coefficients for dense plasmas.  Phys. Rev. E, 69:016409, Jan 2004.
  * [85] P.A. Sterne, S.B. Hansen, B.G. Wilson, and W.A. Isaacs.  Equation of state, occupation probabilities and conductivities in the average atom purgatorio code.  High Energy Density Physics, 3(1):278–282, 2007.  Radiative Properties of Hot Dense Matter.
  * [86] G. Massacrier, M. Böhme, J. Vorberger, F. Soubiran, and B. Militzer.  Reconciling ionization energies and band gaps of warm dense matter derived with ab initio simulations and average atom models.  Phys. Rev. Research, 3:023026, Apr 2021.
  * [87] John P. Perdew and Yue Wang.  Accurate and simple analytic representation of the electron-gas correlation energy.  Phys. Rev. B, 45:13244–13249, Jun 1992.
  * [88] Timothy J. Callow, Eli Kraisler, and Attila Cangi.  Improved calculations of mean ionization states with an average-atom model.  Phys. Rev. Res., 5:013049, Jan 2023.
  * [89] Timothy Callow, Daniel Kotik, Ekaterina Tsvetoslavova Stankulova, Eli Kraisler, and Attila Cangi.  atoMEC, August 2021.  https://doi.org/10.5281/zenodo.5205718.
  * [90] Timothy J. Callow, Daniel Kotik, Eli Kraisler, and Attila Cangi.  atoMEC: An open-source average-atom Python code.  In Meghann Agarwal, Chris Calloway, Dillon Niederhut, and David Shupe, editors, Proceedings of the 21st Python in Science Conference, pages 31 – 39, 2022.
  * [91] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [92] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [93] Susi Lehtola, Conrad Steigemann, Micael J.T. Oliveira, and Miguel A.L. Marques.  Recent developments in libxc — a comprehensive library of functionals for density functional theory.  SoftwareX, 7:1–5, 2018.
  * [94] mendeleev – a python resource for properties of chemical elements, ions and isotopes, ver. 0.9.0.  https://github.com/lmmentel/mendeleev, 2014–.
  * [95] Joblib Development Team.  Joblib: running python functions as pipeline jobs.  https://joblib.readthedocs.io/, 2020.
  * [96] Richard M. More.  Quantum-statistical model for high-density matter.  Phys. Rev. A, 19:1234–1246, Mar 1979.
  * [97] R.M. More.  Pressure ionization, resonances, and the continuity of bound and free states.  volume 21 of Advances in Atomic and Molecular Physics, pages 305–356. Academic Press, 1985.
  * [98] J C Pain.  A model of dense-plasma atomic structure for equation-of-state calculations.  Journal of Physics B: Atomic, Molecular and Optical Physics, 40(8):1553, mar 2007.
  * [99] The ideal approximation, as implied by the name, is a known approximation, and cannot be derived from the functional derivative of the free energy.
  * [100] We note that, in Eq. (12), we have written the xc-free energy $F_{\textrm{xc}}[n]$ as part of the internal energy. However, in principle, it also contains an entropic contribution, because the entropy $S[n]$ is approximated by the non-interacting entropy functional in KS-DFT. We direct readers to Ref. [26] for a more detailed discussion of the xc energy term in ground-state and finite-temperature KS-DFT.
  * [101] Thomas Blenski and Kenichi Ishikawa.  Pressure ionization in the spherical ion-cell model of dense plasmas and a pressure formula in the relativistic pauli approximation.  Phys. Rev. E, 51:4869–4881, May 1995.
  * [102] Nadine Wetta, Jean-Christophe Pain, and Olivier Heuzé.  D’yakov-kontorovitch instability of shock waves in hot plasmas.  Phys. Rev. E, 98:033205, Sep 2018.
  * [103] Gérald Faussurier, Christophe Blancard, and Mandy Bethkenhagen.  Carbon ionization from a quantum average-atom model up to gigabar pressures.  Phys. Rev. E, 104:025209, Aug 2021.
  * [104] John C. Slater.  Hellmann‐feynman and virial theorems in the $x\alpha$ method.  The Journal of Chemical Physics, 57(6):2389–2396, 1972.
  * [105] A. G. McLellan.  Virial theorem generalized.  American Journal of Physics, 42(3):239–243, 1974.
  * [106] P. Ziesche, J. Gräfenstein, and O. H. Nielsen.  Quantum-mechanical stress and a generalized virial theorem for clusters and solids.  Phys. Rev. B, 37:8167–8178, May 1988.
  * [107] P Legrand and F Perrot.  Virial theorem and pressure calculations in the gga.  Journal of Physics: Condensed Matter, 13(2):287, jan 2001.
  * [108] Paul W. Ayers, Robert G. Parr, and Agnes Nagy.  Local kinetic energy and local temperature in the density-functional theory of electronic structure.  International Journal of Quantum Chemistry, 90(1):309–326, 2002\.
  * [109] Leon Cohen.  Local kinetic energy in quantum mechanics.  The Journal of Chemical Physics, 70(2):788–789, 1979.
  * [110] M. Kohout and A. Savin.  Atomic shell structure and electron numbers.  International Journal of Quantum Chemistry, 60(4):875–882, 1996\.
  * [111] Walter Johnson.  Atoms at finite temperatures, August 2000.  [Online; last accessed 26.02.23].
  * [112] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [113] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton.  Imagenet classification with deep convolutional neural networks.  Commun. ACM, 60(6):84–90, may 2017.
  * [114] Tian Xie and Jeffrey C. Grossman.  Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties.  Phys. Rev. Lett., 120:145301, Apr 2018.
  * [115] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural Computation, 9(8):1735–1780, 1997.
  * [116] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  Commun. ACM, 63(11):139–144, oct 2020.
  * [117] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning representations by back-propagating errors.  Nature, 323(6088):533–536, 1986.
  * [118] Trevor Hastie, Robert Tibshirani, and Jerome Friedman.  Elements of Statistical Learning.  Springer, Berlin, 2nd edition, 2009.  12th printing.
  * [119] The reason we use two different notations is because $P_{\textrm{ref}}$ is used to denote the reference pressure in general, and is used later when evaluating the raw AA results, as well as the neural network results. $Y^{0}$ is only used in the context of the neural network training procedure, and is used to denote a specific reference pressure for the given stage / subset of the training workflow. Furthermore, as will be discussed in Section IV.1, $Y^{0}$ is a transformation of the original pressure via a scaling relation.
  * [120] J. Kittler.  Combining classifiers: A theoretical framework.  Pattern Analysis and Applications, 1(1):18–27, 1998.
  * [121] M. G. Kendall.  A new measure of rank correlation.  Biometrika, 30(1/2):81–93, 1938.
  * [122] P. Hohenberg and W. Kohn.  Inhomogeneous electron gas.  Phys. Rev., 136:B864–B871, Nov 1964.
  * [123] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2014.
  * [124] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [125] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019.
  * [126] H.A. Scott and S.B. Hansen.  Advances in NLTE modeling for integrated simulations.  High Energy Density Physics, 6(1):39–47, 2010.
  * [127] Daniel A. Holladay, Christopher J. Fontes, Wesley P. Even, and Ryan G. McClarren.  An accelerated approach to inline non-lte modeling.  High Energy Density Physics, 34:100746, 2020.