# atoMEC: An open-source average-atom Python code

Timothy J. Callow‡§∗, Daniel Kotik‡§, Eli Kraisler¶, Attila Cangi‡§ *
Corresponding author<EMAIL_ADDRESS>Center for Advanced Systems
Understanding (CASUS), D-02826 Görlitz, Germany§ Helmholtz-Zentrum Dresden-
Rossendorf, D-01328 Dresden, Germany¶ Fritz Haber Center for Molecular
Dynamics and Institute of Chemistry, The Hebrew University of Jerusalem,
9091401 Jerusalem, IsraelCopyright © 2022 Timothy J. Callow et al. This is an
open-access article distributed under the terms of the Creative Commons
Attribution License, which permits unrestricted use, distribution, and
reproduction in any medium, provided the original author and source are
credited.

###### Abstract

Average-atom models are an important tool in studying matter under extreme
conditions, such as those conditions experienced in planetary cores, brown and
white dwarfs, and during inertial confinement fusion. In the right context,
average-atom models can yield results with similar accuracy to simulations
which require orders of magnitude more computing time, and thus can greatly
reduce financial and environmental costs. Unfortunately, due to the wide range
of possible models and approximations, and the lack of open-source codes,
average-atom models can at times appear inaccessible. In this paper, we
present our open-source average-atom code, atoMEC. We explain the aims and
structure of atoMEC to illuminate the different stages and options in an
average-atom calculation, and to facilitate community contributions. We also
discuss the use of various open-source Python packages in atoMEC, which have
expedited its development.

###### Index Terms:

computational physics, plasma physics, atomic physics, materials science

## Introduction

The study of matter under extreme conditions — materials exposed to high
temperatures, high pressures, or strong electromagnetic fields — is critical
to our understanding of many important scientific and technological processes,
such as nuclear fusion and various astrophysical and planetary physics
phenomena [GFG+16]. Of particular interest within this broad field is the warm
dense matter (WDM) regime, which is typically characterized by temperatures in
the range of $10^{3}-10^{6}$ degrees (Kelvin), and densities ranging from
dense gases to highly compressed solids ($\sim 0.01-1000\ \textrm{g cm}^{-3}$)
[BDM+20]. In this regime, it is important to account for the quantum
mechanical nature of the electrons (and in some cases, also the nuclei).
Therefore conventional methods from plasma physics, which either neglect
quantum effects or treat them coarsely, are usually not sufficiently accurate.
On the other hand, methods from condensed-matter physics and quantum
chemistry, which account fully for quantum interactions, typically target the
ground-state only, and become computationally intractable for systems at high
temperatures.

Nevertheless, there are methods which can, in principle, be applied to study
materials at any given temperature and density whilst formally accounting for
quantum interactions. These methods are often denoted as "first-principles"
because, formally speaking, they yield the exact properties of the system,
under certain well-founded theoretical approximations. Density-functional
theory (DFT), initially developed as a ground-state theory [HK64, KS65] but
later extended to non-zero temperatures [Mer65, PPF+11], is one such theory
and has been used extensively to study materials under WDM conditions
[GDRT14]. Even though DFT reformulates the Schrödinger equation in a
computationally efficient manner [Koh99], the cost of running calculations
becomes prohibitively expensive at higher temperatures. Formally, it scales as
$\mathcal{O}(N^{3}\tau^{3})$, with $N$ the particle number (which usually also
increases with temperature) and $\tau$ the temperature [CRNB18]. This poses a
serious computational challenge in the WDM regime. Furthermore, although DFT
is a formally exact theory, in practice it relies on approximations for the
so-called "exchange-correlation" energy, which is, roughly speaking,
responsible for simulating all the quantum interactions between electrons.
Existing exchange-correlation approximations have not been rigorously tested
under WDM conditions. An alternative method used in the WDM community is path-
integral Monte–Carlo [DGB18], which yields essentially exact properties;
however, it is even more limited by computational cost than DFT, and becomes
unfeasibly expensive at lower temperatures due to the fermion sign problem.

It is therefore of great interest to reduce the computational complexity of
the aforementioned methods. The use of graphics processing units in DFT
calculations is becomingly increasingly common, and has been shown to offer
significant speed-ups relative to conventional calculations using central
processing units [MED11, JFC+13]. Some other examples of promising
developments to reduce the cost of DFT calculations include machine-learning-
based solutions [SRH+12, BVL+17, EFP+21] and stochastic DFT [CRNB18, BNR13].
However, in this paper, we focus on an alternative class of models known as
"average-atom" models. Average-atom models have a long history in plasma
physics [CHKC22]: they account for quantum effects, typically using DFT, but
reduce the complex system of interacting electrons and nuclei to a single atom
immersed in a plasma (the "average" atom). An illustration of this principle
(reduced to two dimensions for visual purposes) is shown in Fig. 1. This
significantly reduces the cost relative to a full DFT simulation, because the
particle number is restricted to the number of electrons per nucleus, and
spherical symmetry is exploited to reduce the three-dimensional problem to one
dimension.

Naturally, to reduce the complexity of the problem as described, various
approximations must be introduced. It is important to understand these
approximations and their limitations for average-atom models to have genuine
predictive capabilities. Unfortunately, this is not always the case: although
average-atom models share common concepts, there is no unique formal theory
underpinning them. Therefore a variety of models and codes exist, and it is
not typically clear which models can be expected to perform most accurately
under which conditions. In a previous paper [CHKC22], we addressed this issue
by deriving an average-atom model from first principles, and comparing the
impact of different approximations within this model on some common
properties.

In this paper, we focus on computational aspects of average-atom models for
WDM. We introduce atoMEC [CKTS+21]: an open-source average-atom code for
studying Matter under Extreme Conditions. One of the main aims of atoMEC is to
improve the accessibility and understanding of average-atom models. To the
best of our knowledge, open-source average-atom codes are in scarce supply:
with atoMEC, we aim to provide a tool that people can use to run average-atom
simulations and also to add their own models, which should facilitate
comparisons of different approximations. The relative simplicity of average-
atom codes means that they are not only efficient to run, but also efficient
to develop: this means, for example, that they can be used as a test-bed for
new ideas that could be later implemented in full DFT codes, and are also
accessible to those without extensive prior expertise, such as students.
atoMEC aims to facilitate development by following good practice in software
engineering (for example extensive documentation), a careful design structure,
and of course through the choice of Python and its widely used scientific
stack, in particular the NumPy [HMvdW+20] and SciPy [VGO+20] libraries.

This paper is structured as follows: in the next section, we briefly review
the key theoretical points which are important to understand the functionality
of atoMEC, assuming no prior physical knowledge of the reader. Following that,
we present the key functionality of atoMEC, discuss the code structure and
algorithms, and explain how these relate to the theoretical aspects
introduced. Finally, we present an example case study: we consider helium
under the conditions often experienced in the outer layers of a white dwarf
star, and probe the behavior of a few important properties, namely the band-
gap, pressure, and ionization degree.

Figure 1: Illustration of the average-atom concept. The many-body and fully-
interacting system of electron density (shaded blue) and nuclei (red points)
on the left is mapped into the much simpler system of independent atoms on the
right. Any of these identical atoms represents the "average-atom". The effects
of interaction from neighboring atoms are implicitly accounted for in an
approximate manner through the choice of boundary conditions.

## Theoretical background

Properties of interest in the warm dense matter regime include the equation-
of-state data, which is the relation between the density, energy, temperature
and pressure of a material [HRD08]; the mean ionization state and the electron
ionization energies, which tell us about how tightly bound the electrons are
to the nuclei; and the electrical and thermal conductivities. These properties
yield information pertinent to our understanding of stellar and planetary
physics, the Earth’s core, inertial confinement fusion, and more besides. To
exactly obtain these properties, one needs (in theory) to determine the
thermodynamic ensemble of the quantum states (the so-called _wave-functions_)
representing the electrons and nuclei. Fortunately, they can be obtained with
reasonable accuracy using models such as average-atom models; in this section,
we elaborate on how this is done.

We shall briefly review the key theory underpinning the type of average-atom
model implemented in atoMEC. This is intended for readers without a background
in quantum mechanics, to give some context to the purposes and mechanisms of
the code. For a comprehensive derivation of this average-atom model, we direct
readers to Ref. [CHKC22]. The average-atom model we shall describe falls into
a class of models known as _ion-sphere_ models, which are the simplest (and
still most widely used) class of average-atom model. There are alternative
(more advanced) classes of model such as _ion-correlation_ [Roz91] and
_neutral pseudo-atom_ models [SS14] which we have not yet implemented in
atoMEC, and thus we do not elaborate on them here.

As demonstrated in Fig. 1, the idea of the ion-sphere model is to map a fully-
interacting system of many electrons and nuclei into a set of independent
atoms which do not interact explicitly with any of the other spheres.
Naturally, this depends on several assumptions and approximations, but there
is formal justification for such a mapping [CHKC22]. Furthermore, there are
many examples in which average-atom models have shown good agreement with more
accurate simulations and experimental data [FB19], which further justifies
this mapping.

Although the average-atom picture is significantly simplified relative to the
full many-body problem, even determining the wave-functions and their ensemble
weights for an atom at finite temperature is a complex problem. Fortunately,
DFT reduces this complexity further, by establishing that the electron
_density_ — a far less complex entity than the wave-functions — is sufficient
to determine all physical observables. The most popular formulation of DFT,
known as Kohn–Sham DFT (KS-DFT) [KS65], allows us to construct the _fully-
interacting_ density from a _non-interacting_ system of electrons, simplifying
the problem further still. Due to the spherical symmetry of the atom, the non-
interacting electrons — known as KS electrons (or KS orbitals) — can be
represented as a wave-function that is a product of radial and angular
components,

$\phi_{nlm}(\mathbf{r})=X_{nl}(r)Y_{l}^{m}(\theta,\phi)\,,$ (1)

where $n,\ l,\ \textrm{and}\ m$ are the _quantum numbers_ of the orbitals,
which come from the fact that the wave-function is an eigenfunction of the
Hamiltonian operator, and $Y_{l}^{m}(\theta,\phi)$ are the spherical harmonic
functions.1 The radial coordinate $r$ represents the absolute distance from
the
nucleus.\raisebox{10.00002pt}{\hypertarget{f2}{}}\hyperlink{id1}{1}\raisebox{10.00002pt}{\hypertarget{f2}{}}\hyperlink{id1}{1}footnotetext:
Please note that the notation in Eq. (1) does not imply Einstein summation
notation. All summations in this paper are written explicitly; Einstein
summation notation is not used.

We therefore only need to determine the radial KS orbitals $X_{nl}(r)$. These
are determined by solving the radial KS equation, which is similar to the
Schrödinger equation for a non-interacting system, with an additional term in
the potential to mimic the effects of electron-electron interaction (within
the single atom). The radial KS equation is given by:

$\left[-\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)+v_{\textrm{s}}[n](r)\right]X_{nl}(r)=\epsilon_{nl}X_{nl}(r).$
(2)

We have written the above equation in a way that emphasizes that it is an
eigenvalue equation, with the eigenvalues $\epsilon_{nl}$ being the energies
of the KS orbitals.

On the left-hand side, the terms in the round brackets come from the kinetic
energy operator acting on the orbitals. The $v_{\textrm{s}}[n](r)$ term is the
KS potential, which itself is composed of three different terms,

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

where $R_{\textrm{WS}}$ is the radius of the atomic sphere, $n(r)$ is the
electron density, $Z$ the nuclear charge, and $F_{\textrm{xc}}[n]$ the
exchange-correlation free energy functional. Thus the three terms in the
potential are respectively the electron-nuclear attraction, the classical
Hartree repulsion, and the exchange-correlation (xc) potential.

We note that the KS potential and its constituents are functionals of the
electron density $n(r)$. Were it not for this dependence on the density,
solving Eq. 2 just amounts to solving an ordinary linear differential equation
(ODE). However, the electron density is in fact constructed from the orbitals
in the following way,

$n(r)=2\sum_{nl}(2l+1)f_{nl}(\epsilon_{nl},\mu,\tau)|X_{nl}(r)|^{2}\,,$ (4)

where $f_{nl}(\epsilon_{nl},\mu,\tau)$ is the Fermi–Dirac distribution, given
by

$f_{nl}(\epsilon_{nl},\mu,\tau)=\frac{1}{1+e^{(\epsilon_{nl}-\mu)/\tau}}\,,$
(5)

where $\tau$ is the temperature, and $\mu$ is the chemical potential, which is
determined by fixing the number of electrons to be equal to a pre-determined
value $N_{\textrm{e}}$ (typically equal to the nuclear charge $Z$). The
Fermi–Dirac distribution therefore assigns weights to the KS orbitals in the
construction of the density, with the weight depending on their energy.

Therefore, the KS potential that determines the KS orbitals via the ODE (2),
is itself dependent on the KS orbitals. Consequently, the KS orbitals and
their dependent quantities (the density and KS potential) must be determined
via a so-called self-consistent field (SCF) procedure. An initial guess for
the orbitals, $X_{nl}^{0}(r)$, is used to construct the initial density
$n^{0}(r)$ and potential $v_{\textrm{s}}^{0}(r)$. The ODE (2) is then solved
to update the orbitals. This process is iterated until some appropriately
chosen quantities — in atoMEC the total free energy, density and KS potential
— are converged, i.e. $n^{i+1}(r)=n^{i}(r),\
v_{\textrm{s}}^{i+1}(r)=v_{\textrm{s}}^{i}(r),\ F^{i+1}=F^{i}$, within some
reasonable numerical tolerance. In Fig. 2, we illustrate the life-cycle of the
average-atom model described so far, including the SCF procedure. On the left-
hand side of this figure, we show the physical choices and mathematical
operations, and on the right-hand side, the representative classes and
functions in atoMEC. In the following section, we shall discuss some aspects
of this figure in more detail.

Some quantities obtained from the completion of the SCF procedure are directly
of interest. For example, the energy eigenvalues $\epsilon_{nl}$ are related
to the electron ionization energies, i.e. the amount of energy required to
excite an electron bound to the nucleus to being a free (conducting) electron.
These predicted ionization energies can be used, for example, to help
understand ionization potential depression, an important but somewhat
controversial effect in WDM [STJ+14]. Another property that can be
straightforwardly obtained from the energy levels and their occupation numbers
is the mean ionization state $\bar{Z}$2,

$\bar{Z}=\sum_{n,l}(2l+1)f_{nl}(\epsilon_{nl},\mu,\tau)$ (6)

which is an important input parameter for various models, such as adiabats
which are used to model inertial confinement fusion
[KDF+11].\raisebox{10.00002pt}{\hypertarget{f1}{}}\hyperlink{id2}{2}\raisebox{10.00002pt}{\hypertarget{f1}{}}\hyperlink{id2}{2}footnotetext:
The summation in Eq. (6) is often shown as an integral because the energies
above a certain threshold form a continuous distribution (in most models).

Various other interesting properties can also be calculated following some
post-processing of the output of an SCF calculation, for example the pressure
exerted by the electrons and ions. Furthermore, response properties, i.e.
those resulting from an external perturbation like a laser pulse, can also be
obtained from the output of an SCF cycle. These properties include, for
example, electrical conductivities [Sta16] and dynamical structure factors
[SPS+14].

Figure 2: Schematic of the average-atom model set-up and the self-consistent
field (SCF) cycle. On the left-hand side, the physical choices and
mathematical operations that define the model and SCF cycle are shown. On the
right-hand side, the (higher-order) functions and classes in atoMEC
corresponding to the items on the left-hand side are shown. Some liberties are
taken with the code snippets in the right-hand column of the figure to improve
readability; more precisely, some non-crucial intermediate steps are not
shown, and some parameters are also not shown or simplified. The dotted lines
represent operations that are taken care of within the models.CalcEnergy
function, but are shown nevertheless to improve understanding.

## Code structure and details

In the following sections, we describe the structure of the code in relation
to the physical problem being modeled. Average-atom models typically rely on
various parameters and approximations. In atoMEC, we have tried to structure
the code in a way that makes clear which parameters come from the physical
problem studied compared to choices of the model and numerical or algorithmic
choices.

###  atoMEC.Atom: Physical parameters

The first step of any simulation in WDM (which also applies to simulations in
science more generally) is to define the physical parameters of the problem.
These parameters are unique in the sense that, if we had an exact method to
simulate the real system, then for each combination of these parameters there
would be a unique solution. In other words, regardless of the model — be it
average atom or a different technique — these parameters are always required
and are independent of the model.

In average-atom models, there are typically three parameters defining the
physical problem, which are:

  * •

the atomic species;

  * •

the temperature of the material, $\tau$;

  * •

the mass density of the material, $\rho_{\textrm{m}}$.

The mass density also directly corresponds to the mean distance between two
nuclei (atomic centers), which in the average-atom model is equal to twice the
radius of the atomic sphere, $R_{\textrm{WS}}$. An additional physical
parameter not mentioned above is the net charge of the material being
considered, i.e. the difference between the nuclear charge $Z$ and the
electron number $N_{\textrm{e}}$. However, we usually assume zero net charge
in average-atom simulations (i.e. the number of electrons is equal to the
atomic charge).

In atoMEC, these physical parameters are controlled by the Atom object. As an
example, we consider aluminum under ambient conditions, i.e. at room
temperature, $\tau=300\ \textrm{K}$, and normal metallic density,
$\rho_{\textrm{m}}=2.7\ \textrm{g cm}^{-3}$. We set this up as:

from atoMEC import AtomAl = Atom("Al", 300, density=2.7, units_temp="K")

Figure 3: Auto-generated print statement from calling the atoMEC.Atom object.

By default, the above code automatically prints the output seen in Fig. 3. We
see that the first two arguments of the Atom object are the chemical symbol of
the element being studied, and the temperature. In addition, at least one of
"density" or "radius" must be specified. In atoMEC, the default (and only
permitted) units for the mass density are $\textrm{g cm}^{-3}$; _all_ other
input and output units in atoMEC are by default Hartree atomic units, and
hence we specify "K" for Kelvin.

The information in Fig. 3 displays the chosen parameters in units commonly
used in the plasma and condensed-matter physics communities, as well as some
other information directly obtained from these parameters. The chemical symbol
("Al" in this case) is passed to the mendeleev library [men14] to generate
this data, which is used later in the calculation.

This initial stage of the average-atom calculation, i.e. the specification of
physical parameters and initialization of the Atom object, is shown in the top
row at the top of Fig. 2.

###  atoMEC.models: Model parameters

After the physical parameters are set, the next stage of the average-atom
calculation is to choose the model and approximations within that class of
model. As discussed, so far the only class of model implemented in atoMEC is
the ion-sphere model. Within this model, there are still various choices to be
made by the user. In some cases, these choices make little difference to the
results, but in other cases they have significant impact. The user might have
some physical intuition as to which is most important, or alternatively may
want to run the same physical parameters with several different model
parameters to examine the effects. Some choices available in atoMEC, listed
approximately in decreasing order of impact (but this can depend strongly on
the system under consideration), are:

  * •

the boundary conditions used to solve the KS equations;

  * •

the treatment of the unbound electrons, which means those electrons not
tightly bound to the nucleus, but rather delocalized over the whole atomic
sphere;

  * •

the choice of exchange and correlation functionals, the central approximations
of DFT [CMSY12];

  * •

the spin polarization and magnetization.

We do not discuss the theory and impact of these different choices in this
paper. Rather, we direct readers to Refs. [CHKC22] and [CKC22] in which all of
these choices are discussed.

In atoMEC, the ion-sphere model is controlled by the models.ISModel object.
Continuing with our aluminum example, we choose the so-called "neumann"
boundary condition, with a "quantum" treatment of the unbound electrons, and
choose the LDA exchange functional (which is also the default). This model is
set up as:

from atoMEC import modelsmodel = models.ISModel(Al, bc="neumann",
xfunc_id="lda_x", unbound="quantum") By default, the above code prints the
output shown in Fig. 4. The first (and only mandatory) input parameter to the
models.ISModel object is the Atom object that we generated earlier. Together
with the optional spinpol and spinmag parameters in the models.ISModel object,
this sets either the total number of electrons (spinpol=False) or the number
of electrons in each spin channel (spinpol=True).

The remaining information displayed in Fig. 4 shows directly the chosen model
parameters, or the default values where these parameters are not specified.
The exchange and correlation functionals - set by the parameters xfunc_ id and
cfunc_ id \- are passed to the LIBXC library [LSOM18] for processing. So far,
only the "local density" family of approximations is available in atoMEC, and
thus the default values are usually a sensible choice. For more information on
exchange and correlation functionals, there are many reviews in the
literature, for example Ref. [CMSY12].

This stage of the average-atom calculation, i.e. the specification of the
model and the choices of approximation within that, is shown in the second row
of Fig. 2.

Figure 4: Auto-generated print statement from calling the models.ISModel
object.

###  ISModel.CalcEnergy: SCF calculation and numerical parameters

Once the physical parameters and model have been defined, the next stage in
the average-atom calculation (or indeed any DFT calculation) is the SCF
procedure. In atoMEC, this is invoked by the ISModel.CalcEnergy function. This
function is called CalcEnergy because it finds the KS orbitals (and associated
KS density) which minimize the total free energy.

Clearly, there are various mathematical and algorithmic choices in this
calculation. These include, for example: the basis in which the KS orbitals
and potential are represented, the algorithm used to solve the KS equations
(2), and how to ensure smooth convergence of the SCF cycle. In atoMEC, the SCF
procedure currently follows a single pre-determined algorithm, which we
briefly review below.

In atoMEC, we represent the radial KS quantities (orbitals, density and
potential) on a logarithmic grid, i.e. $x=\log(r)$. Furthermore, we make a
transformation of the orbitals $P_{nl}(x)=X_{nl}(x)e^{x/2}$. Then the
equations to be solved become:

$\displaystyle\frac{\textrm{d}^{2}P_{nl}(x)}{\textrm{d}x^{2}}-2e^{2x}(W(x)-\epsilon_{nl})P_{nl}(x)=0$
(7) $\displaystyle
W(x)=v_{\textrm{s}}[n](x)+\frac{1}{2}\left(l+\frac{1}{2}\right)^{2}e^{-2x}\,.$
(8)

In atoMEC, we solve the KS equations using a matrix implementation of
Numerov’s algorithm [PGW12]. This means we diagonalize the following equation:

$\displaystyle\hat{H}\vec{P}$ $\displaystyle=\vec{\epsilon}\hat{B}\vec{P}\,,\
\textrm{where}$ (9) $\displaystyle\hat{H}$
$\displaystyle=\hat{T}+\hat{B}+W_{\textrm{s}}(\vec{x})\,,$ (10)
$\displaystyle\hat{T}$ $\displaystyle=-\frac{1}{2}e^{-2\vec{x}}\hat{A}\,,$
(11) $\displaystyle\hat{A}$
$\displaystyle=\frac{\hat{I}_{-1}-2\hat{I}_{0}+\hat{I}_{1}}{\textrm{d}x^{2}}\,,\textrm{and}$
(12) $\displaystyle\hat{B}$
$\displaystyle=\frac{\hat{I}_{-1}+10\hat{I}_{0}+\hat{I}_{1}}{12}\,,$ (13)

In the above, $\hat{I}_{-1/0/1}$ are lower shift, identify, and upper shift
matrices.

The Hamiltonian matrix $\hat{H}$ is sparse and we only seek a subset of
eigenstates with lower energies: therefore there is no need to perform a full
diagonalization, which scales as $\mathcal{O}(N^{3})$, with $N$ being the size
of the radial grid. Instead, we use SciPy’s sparse matrix diagonalization
function scipy.sparse.linalg.eigs, which scales more efficiently and allows us
to go to larger grid sizes.

After each step in the SCF cycle, the relative changes in the free energy $F$,
density $n(r)$ and potential $v_{\textrm{s}}(r)$ are computed. Specifically,
the quantities computed are

$\displaystyle\Delta F$
$\displaystyle=\left|\frac{F^{i}-F^{i-1}}{F^{i}}\right|$ (14)
$\displaystyle\Delta n$
$\displaystyle=\frac{\int\mathrm{d}r|n^{i}(r)-n^{i-1}(r)|}{\int\mathrm{d}rn^{i}(r)}$
(15) $\displaystyle\Delta v$
$\displaystyle=\frac{\int\mathrm{d}r|v^{i}_{\textrm{s}}(r)-v_{\textrm{s}}^{i-1}(r)|}{\int\mathrm{d}rv_{\textrm{s}}^{i}(r)}\,.$
(16)

Once all three of these metrics fall below a certain threshold, the SCF cycle
is considered converged and the calculation finishes.

The SCF cycle is an example of a non-linear system and thus is prone to
chaotic (non-convergent) behavior. Consequently a range of techniques have
been developed to ensure convergence [SM91]. Fortunately, the tendency for
calculations not to converge becomes less likely for temperatures above zero
(and especially as temperatures increase). Therefore we have implemented only
a simple linear mixing scheme in atoMEC. The potential used in each
diagonalization step of the SCF cycle is not simply the one generated from the
most recent density, but a mix of that potential and the previous one,

$v_{\textrm{s}}^{(i)}(r)=\alpha
v_{\textrm{s}}^{i}(r)+(1-\alpha)v_{\textrm{s}}^{i-1}(r)\,.$ (17)

In general, a lower value of the mixing fraction $\alpha$ makes the SCF cycle
more stable, but requires more iterations to converge. Typically a choice of
$\alpha\approx 0.5$ gives a reasonable balance between speed and stability.

We can thus summarize the key parameters in an SCF calculation as follows:

  * •

the maximum number of eigenstates to compute, in terms of both the principal
and angular quantum numbers;

  * •

the numerical grid parameters, in particular the grid size;

  * •

the convergence tolerances, Eqs. (14) to (16);

  * •

the SCF parameters, i.e. the mixing fraction and the maximum number of
iterations.

The first three items in this list essentially control the accuracy of the
calculation. In principle, for each SCF calculation — i.e. a unique set of
physical and model inputs — these parameters should be independently varied
until some property (such as the total free energy) is considered suitably
converged with respect to that parameter. Changing the SCF parameters should
not affect the final results (within the convergence tolerances), only the
number of iterations in the SCF cycle.

Let us now consider an example SCF calculation, using the Atom and model
objects we have already defined:

from atoMEC import configconfig.numcores = -1 # parallelizenmax = 3 # max
value of principal quantum numberlmax = 3 # max value of angular quantum
number# run SCF calculationscf_out = model.CalcEnergy( nmax, lmax,
grid_params={"ngrid": 1500}, scf_params={"mixfrac": 0.7}, ) We see that the
first two parameters passed to the CalcEnergy function are the nmax and lmax
quantum numbers, which specify the number of eigenstates to compute. Precisely
speaking, there is a unique Hamiltonian for each value of the angular quantum
number $l$ (and in a spin-polarized calculation, also for each spin quantum
number). The sparse diagonalization routine then computes the first nmax
eigenvalues for each Hamiltonian. In atoMEC, these diagonalizations can be run
in parallel since they are independent for each value of $l$. This is done by
setting the config.numcores variable to the number of cores desired
(config.numcores=-1 uses all the available cores) and handled via the joblib
library [Job20].

The remaining parameters passed to the CalcEnergy function are optional; in
the above, we have specified a grid size of 1500 points and a mixing fraction
$\alpha=0.7$. The above code automatically prints the output seen in Fig. 5.
This output shows the SCF cycle and, upon completion, the breakdown of the
total free energy into its various components, as well as other useful
information such as the KS energy levels and their occupations.

Figure 5: Auto-generated print statement from calling the ISModel.CalcEnergy
function

Additionally, the output of the SCF function is a dictionary containing the
staticKS.Orbitals, staticKS.Density, staticKS.Potential and staticKS.Density
objects. For example, one could extract the eigenfunctions as follows:

orbs = scf_out["orbitals"] # orbs objectks_eigfuncs = orbs.eigfuncs #
eigenfunctions The initialization of the SCF procedure is shown in the third
and fourth rows of Fig. 2, with the SCF procedure itself shown in the
remaining rows.

This completes the section on the code structure and algorithmic details. As
discussed, with the output of an SCF calculation, there are various kinds of
post-processing one can perform to obtain other properties of interest. So far
in atoMEC, these are limited to the computation of the pressure
(ISModel.CalcPressure), the electron localization function
(atoMEC.postprocess.ELFTools) and the Kubo–Greenwood conductivity
(atoMEC.postprocess.conductivity). We refer readers to our pre-print [CKC22]
for details on how the electron localization function and the Kubo–Greenwood
conductivity can be used to improve predictions of the mean ionization state.

## Case-study: Helium

In this section, we consider an application of atoMEC in the WDM regime.
Helium is the second most abundant element in the universe (after hydrogen)
and therefore understanding its behavior under a wide range of conditions is
important for our understanding of many astrophysical processes. Of particular
interest are the conditions under which helium is expected to undergo a
transition from insulating to metallic behavior in the outer layers of white
dwarfs, which are characterized by densities of around $1-20\textrm{ g
cm}^{-3}$ and temperatures of $10-50$ kK [PR20]. These conditions are a
typical example of the WDM regime. Besides predicting the point at which the
insulator-to-metallic transition occurs in the density-temperature spectrum,
other properties of interest include equation-of-state data (relating
pressure, density, and temperature) and electrical conductivity.

To calculate the insulator-to-metallic transition point, the key quantity is
the electronic _band-gap_. The concept of band-structures is a complicated
topic, which we try to briefly describe in layman’s terms. In solids,
electrons can occupy certain energy ranges — we call these the energy bands.
In insulating materials, there is a gap between these energy ranges that
electrons are forbidden from occupying — this is the so-called band-gap. In
conducting materials, there is no such gap, and therefore electrons can
conduct electricity because they can be excited into any part of the energy
spectrum. Therefore, a simple method to determine the insulator-to-metallic
transition is to determine the density at which the band-gap becomes zero.

In Fig. 6, we plot the density-of-states (DOS) as a function of energy, for
different densities and at fixed temperature $\tau=50$ kK. The DOS shows the
energy ranges that the electrons are allowed to occupy; we also show the
actual energies occupied by the electrons (according to Fermi–Dirac
statistics) with the black dots. We can clearly see in this figure that the
band-gap (the region where the DOS is zero) becomes smaller as a function of
density. From this figure, it seems the transition from insulating to metallic
state happens somewhere between 5 and 6 $\textrm{g cm}^{-3}.$

Figure 6: Helium density-of-states (DOS) as a function of energy, for
different mass densities $\rho_{\textrm{m}}$, and at temperature $\tau=50$ kK.
Black dots indicate the occupations of the electrons in the permitted energy
ranges. Dashed black lines indicate the band-gap (the energy gap between the
insulating and conducting bands). Between 5 and 6 $\textrm{g cm}^{-3}$, the
band-gap disappears.

In Fig. 7, we plot the band-gap as a function of density, for a fixed
temperature $\tau=50$ kK. Visually, it appears that the relationship between
band-gap and density is linear at this temperature. This is confirmed using a
linear fit, which has a coefficient of determination value of almost exactly
one, $R^{2}=0.9997$. Using this fit, the band-gap is predicted to close at
$5.5\ \textrm{g cm}^{-3}$. Also in this figure, we show the fraction of
ionized electrons, which is given by $\bar{Z}/N_{\textrm{e}}$, using Eq. (6)
to calculate $\bar{Z}$, and $N_{\textrm{e}}$ being the total electron number.
The ionization fraction also relates to the conductivity of the material,
because ionized electrons are not bound to any nuclei and therefore free to
conduct electricity. We see that the ionization fraction mostly increases with
density (excepting some strange behavior around $\rho_{\textrm{m}}=1\
\textrm{g cm}^{-3}$), which is further evidence of the transition from
insulating to conducting behaviour with increasing density.

Figure 7: Band-gap (red circles) and ionization fraction (blue squares) for
helium as a function of mass density, at temperature $\tau=50$ kK. The
relationship between the band-gap and the density appears to be linear.

As a final analysis, we plot the pressure as a function of mass density and
temperature in Fig. 8. The pressure is given by the sum of two terms: (i) the
electronic pressure, calculated using the method described in Ref. [FB19], and
(ii) the ionic pressure, calculated using the ideal gas law. We observe that
the pressure increases with both density and temperature, which is the
expected behavior. Under these conditions, the density dependence is much
stronger, especially for higher densities.

The code required to generate the above results and plots can be found in this
repository.

Figure 8: Helium pressure (logarithmic scale) as a function of mass density
and temperature. The pressure increases with density and temperature (as
expected), with a stronger dependence on density.

## Conclusions and future work

In this paper, we have presented atoMEC: an average-atom Python code for
studying materials under extreme conditions. The open-source nature of atoMEC,
and the choice to use (pure) Python as the programming language, is designed
to improve the accessibility of average-atom models.

We gave significant attention to the code structure in this paper, and tried
as much as possible to connect the functions and objects in the code with the
underlying theory. We hope that this not only improves atoMEC from a user
perspective, but also facilitates new contributions from the wider average-
atom, WDM and scientific Python communities. Another aim of the paper was to
communicate how atoMEC benefits from a strong ecosystem of open-source
scientific libraries — especially the Python libraries NumPy, SciPy, joblib
and mendeleev, as well as LIBXC.

We finish this paper by emphasizing that atoMEC is still in the early stages
of development, and there are many opportunities to improve and extend the
code. These include, for example:

  * •

adding new average-atom models, and different approximations to the existing
models.ISModel model;

  * •

optimizing the code, in particular the routines in the numerov module;

  * •

adding new postprocessing functionality, for example to compute structure
factors;

  * •

improving the structure and design choices of the code.

Of course, these are just a snapshot of the avenues for future development in
atoMEC. We are open to contributions in these areas and many more besides.

## Acknowledgements

This work was partly funded 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 Ministry for Science, Culture and Tourism
(SMWK) with tax funds on the basis of the budget approved by the Saxon State
Parliament.

## References

  * [BDM+20] 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.  Phys. Plasmas, 27(4):042710, 2020.  doi:10.1063/1.5143225.
  * [BNR13] Roi Baer, Daniel Neuhauser, and Eran Rabani.  Self-averaging stochastic Kohn-Sham density-functional theory.  Phys. Rev. Lett., 111:106402, Sep 2013.  doi:10.1103/PhysRevLett.111.106402.
  * [BVL+17] Felix Brockherde, Leslie Vogt, Li Li, Mark E. Tuckerman, Kieron Burke, and Klaus-Robert Müller.  Bypassing the Kohn-Sham equations with machine learning.  Nature Communications, 8(1):872, Oct 2017.  doi:10.1038/s41467-017-00839-3.
  * [CHKC22] 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.  doi:10.1103/PhysRevResearch.4.023055.
  * [CKC22] Timothy J. Callow, Eli Kraisler, and Attila Cangi.  Accurate and efficient computation of mean ionization states with an average-atom Kubo-Greenwood approach, 2022.  doi:10.48550/ARXIV.2203.05863.
  * [CKTS+21] Timothy Callow, Daniel Kotik, Ekaterina Tsvetoslavova Stankulova, Eli Kraisler, and Attila Cangi.  atomec, August 2021.  If you use this software, please cite it using these metadata.  doi:10.5281/zenodo.5205719.
  * [CMSY12] Aron J. Cohen, Paula Mori-Sánchez, and Weitao Yang.  Challenges for density functional theory.  Chemical Reviews, 112(1):289–320, 2012.  doi:10.1021/cr200107z.
  * [CRNB18] Yael Cytter, Eran Rabani, Daniel Neuhauser, and Roi Baer.  Stochastic density functional theory at finite temperatures.  Phys. Rev. B, 97:115207, Mar 2018.  doi:10.1103/PhysRevB.97.115207.
  * [DGB18] Tobias Dornheim, Simon Groth, and Michael Bonitz.  The uniform electron gas at warm dense matter conditions.  Phys. Rep., 744:1 – 86, 2018.  doi:10.1016/j.physrep.2018.04.001.
  * [EFP+21] 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.  doi:10.1103/PhysRevB.104.035120.
  * [FB19] 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.  doi:10.1103/PhysRevE.99.053201.
  * [GDRT14] Frank Graziani, Michael P Desjarlais, Ronald Redmer, and Samuel B Trickey.  Frontiers and challenges in warm dense matter, volume 96.  Springer Science & Business, 2014.  doi:10.1007/978-3-319-04912-0.
  * [GFG+16] S H Glenzer, L B Fletcher, E Galtier, B Nagler, R Alonso-Mori, B Barbrel, S B Brown, D A Chapman, Z Chen, C B Curry, F Fiuza, E Gamboa, M Gauthier, D O Gericke, A Gleason, S Goede, E Granados, P Heimann, J Kim, D Kraus, M J MacDonald, A J Mackinnon, R Mishra, A Ravasio, C Roedel, P Sperling, W Schumaker, Y Y Tsui, J Vorberger, U Zastrau, A Fry, W E White, J B Hasting, and H J Lee.  Matter under extreme conditions experiments at the Linac Coherent Light Source.  J. Phys. B, 49(9):092001, apr 2016.  doi:10.1088/0953-4075/49/9/092001.
  * [HK64] P. Hohenberg and W. Kohn.  Inhomogeneous electron gas.  Phys. Rev., 136(3B):B864–B871, Nov 1964.  doi:10.1103/PhysRev.136.B864.
  * [HMvdW+20] 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.  doi:10.1038/s41586-020-2649-2.
  * [HRD08] Bastian Holst, Ronald Redmer, and Michael P. Desjarlais.  Thermophysical properties of warm dense hydrogen using quantum molecular dynamics simulations.  Phys. Rev. B, 77:184201, May 2008.  doi:10.1103/PhysRevB.77.184201.
  * [JFC+13] Weile Jia, Jiyun Fu, Zongyan Cao, Long Wang, Xuebin Chi, Weiguo Gao, and Lin-Wang Wang.  Fast plane wave density functional theory molecular dynamics calculations on multi-GPU machines.  Journal of Computational Physics, 251:102–115, 2013.  doi:10.1016/j.jcp.2013.05.005.
  * [Job20] Joblib Development Team.  Joblib: running Python functions as pipeline jobs.  https://joblib.readthedocs.io/, 2020.
  * [KDF+11] 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.  doi:10.1103/PhysRevLett.107.015002.
  * [Koh99] W. Kohn.  Nobel lecture: Electronic structure of matter—wave functions and density functionals.  Rev. Mod. Phys., 71:1253–1266, 10 1999.  doi:10.1103/RevModPhys.71.1253.
  * [KS65] W. Kohn and L. J. Sham.  Self-consistent equations including exchange and correlation effects.  Phys. Rev., 140(4A):A1133–A1138, Nov 1965.  doi:10.1103/PhysRev.140.A1133.
  * [LSOM18] 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.  doi:10.1016/j.softx.2017.11.002.
  * [MED11] Stefan Maintz, Bernhard Eck, and Richard Dronskowski.  Speeding up plane-wave electronic-structure calculations using graphics-processing units.  Computer Physics Communications, 182(7):1421–1427, 2011.  doi:10.1016/j.cpc.2011.03.010.
  * [men14] mendeleev – A Python resource for properties of chemical elements, ions and isotopes, ver. 0.9.0.  https://github.com/lmmentel/mendeleev, 2014.
  * [Mer65] N. David Mermin.  Thermal properties of the inhomogeneous electron gas.  Phys. Rev., 137:A1441–A1443, Mar 1965.  doi:10.1103/PhysRev.137.A1441.
  * [PGW12] Mohandas Pillai, Joshua Goglio, and Thad G. Walker.  Matrix numerov method for solving schrödinger’s equation.  American Journal of Physics, 80(11):1017–1019, 2012.  doi:10.1119/1.4748813.
  * [PPF+11] S. Pittalis, C. R. Proetto, A. Floris, A. Sanna, C. Bersier, K. Burke, and E. K. U. Gross.  Exact conditions in finite-temperature density-functional theory.  Phys. Rev. Lett., 107:163001, Oct 2011.  doi:10.1103/PhysRevLett.107.163001.
  * [PR20] Martin Preising and Ronald Redmer.  Metallization of dense fluid helium from ab initio simulations.  Phys. Rev. B, 102:224107, Dec 2020.  doi:10.1103/PhysRevB.102.224107.
  * [Roz91] Balazs F. Rozsnyai.  Photoabsorption in hot plasmas based on the ion-sphere and ion-correlation models.  Phys. Rev. A, 43:3035–3042, Mar 1991.  doi:10.1103/PhysRevA.43.3035.
  * [SM91] H. B. Schlegel and J. J. W. McDouall.  Do You Have SCF Stability and Convergence Problems?, pages 167–185.  Springer Netherlands, Dordrecht, 1991.  doi:10.1007/978-94-011-3262-6_2.
  * [SPS+14] A. N. Souza, D. J. Perkins, C. E. Starrett, D. Saumon, and S. B. Hansen.  Predictions of x-ray scattering spectra for warm dense matter.  Phys. Rev. E, 89:023108, Feb 2014.  doi:10.1103/PhysRevE.89.023108.
  * [SRH+12] John C. Snyder, Matthias Rupp, Katja Hansen, Klaus-Robert Müller, and Kieron Burke.  Finding density functionals with machine learning.  Phys. Rev. Lett., 108:253002, Jun 2012.  doi:10.1103/PhysRevLett.108.253002.
  * [SS14] C.E. Starrett and D. Saumon.  A simple method for determining the ionic structure of warm dense matter.  High Energy Density Physics, 10:35–42, 2014.  doi:10.1016/j.hedp.2013.12.001.
  * [Sta16] C.E. Starrett.  Kubo–Greenwood approach to conductivity in dense plasmas with average atom models.  High Energy Density Physics, 19:58–64, 2016.  doi:10.1016/j.hedp.2016.04.001.
  * [STJ+14] Sang-Kil Son, Robert Thiele, Zoltan Jurek, Beata Ziaja, and Robin Santra.  Quantum-mechanical calculation of ionization-potential lowering in dense plasmas.  Phys. Rev. X, 4:031004, Jul 2014.  doi:10.1103/PhysRevX.4.031004.
  * [VGO+20] 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.  doi:10.1038/s41592-019-0686-2.