[1]Mgcini Keith Phuthi

[1]Venkatasubramanian Viswanathan

[1]Mechanical Engineering, University of Michigan, 500 State St., Ann Arbor,
48105, Michigan, USA

2]Department of Physics, Carnegie Mellon University, 5000 Forbes Ave.,
Pittsburgh, 15213, Pennsylvania, USA

3] University of Science and Technology of China, Hefei, 230026, China 4]
Suzhou Institute for Advanced Research, University of Science and Technology
of China, Suzhou, 215213, China

# Vibrational Entropy and Free Energy of Solid Lithium using Covariance of
Atomic Displacements Enabled by Machine Learning

<EMAIL_ADDRESS>Yang Huang<EMAIL_ADDRESS>Michael Widom
<EMAIL_ADDRESS><EMAIL_ADDRESS>*  [  [  [

###### Abstract

Vibrational properties of solids are key to determining stability, response
and functionality. However, they are challenging to computationally predict at
Ab-Initio accuracy, even for elemental systems. Ab-Initio methods for modeling
atomic interactions are limited in the system sizes and simulation times that
can be achieved. Due to these limitations, Machine Learning Interatomic
Potentials (MLIPs) are gaining popularity and success as a faster, more
scalable approach for modeling atomic interactions, potentially at Ab-Initio
accuracy. Even with faster potentials, methodologies for predicting entropy,
free energy and vibrational properties vary in accuracy, cost and difficulty
to implement. Using the Covariance of Atomic Displacements (CAD) to predict
entropy, free energy and finite-temperature phonon dispersions is a promising
approach but thorough benchmarking has been hampered by the cost of Ab-Initio
methods for sampling. In this work, we use a MLIP and the CAD to characterize
the convergence of the predicted properties and determine optimal sampling
strategies. We focus on solid lithium at zero pressure, showing that the MLIP-
CAD approach reproduces experimental entropy, phonon dispersions and the
martensitic transition while also comparing to more established methods.

###### keywords:

Vibrational Entropy, Free Energy, Machine Learning Interatomic Potential,
Lithium

##  1 Introduction

Computational materials discovery and characterization are powerful tools for
engineering materials for critical applications such as sustainable energy [1,
2]. One of the most important properties of a material is the free energy
landscape as a function of chemical composition, temperature and pressure. The
free energy landscape determines the (meta-)stability of a material and is a
strong indicator for synthesizability [3]. After establishing stability, the
material can be characterized for various other functional properties. There
exist a number of methods of determining vibrational free energies depending
on the material class. For solids, some of the most commonly used are
Thermodynamic Integration (TI) [4], Temperature-Dependent Effective Potentials
(TDEP) [5, 6] and the Quasiharmonic Approximation (QHA) [7]. These methods
employ varying degrees of approximation and can sometimes be difficult to
implement for some material classes such as alloys [8]. Their feasibility is
also determined by the implementation details or computational cost of the
interatomic potential used to model the atomic interactions [9].

A host of methods lie on the Pareto frontier of computational cost vs accuracy
for modeling interatomic interactions, with some shown in Fig. 1. Quantum
chemistry methods are generally the most accurate but infeasible for most
materials of interest [10]. Density Functional Theory (DFT) is more feasible
and extensively used due to explicitly including important quantum mechanical
effects and for its generalizability across the entire periodic table with
acceptable accuracy for most applications[10]. DFT however suffers from
finite-size effects in simulations where long-range correlations matter since
it scales poorly with the number of electrons ($\mathcal{O}(N^{3})$) and is
usually too slow for simulation protocols that require long timescale dynamics
for the statistics of the calculated property to converge [11].

In cases where statistics and/or long-range correlations are significant,
empirical potentials and Machine Learning Interatomic Potentials (MLIPs)
fitted to experimental and/or ab-initio data are more suitable at the cost of
some accuracy [4]. The better scaling with the number of atoms
($\mathcal{O}(N)$) and orders of magnitude speed-up in property evaluation
allows for collection of better statistics. Empirical potentials are however
sometimes only qualitative in their predictions as they do not explicitly
incorporate complicated quantum mechanical interactions. MLIPs on the other
hand can be trained on ab-initio data to reproduce ab-initio predictions to
remarkable accuracy [12, 9]. They however suffer from poor generalizability in
extrapolating outside the domain of the training data [13, 14]. MLIPs are
therefore most reliable when predicting in-distribution properties in cases
where collecting statistics (in size or simulation time) is the limiting
factor for DFT or higher-order methods. The methodology described in this work
takes advantage of the this ability to perform accurate in-distribution
sampling.

Once an interatomic potential is chosen, a method needs to be selected to
perform the vibrational free energy calculations. Thermodynamic Integration
methods are often considered the gold standard as they are formally exact in
predicting free energy differences [7]. However, TI is generally challenging
due to convergence issues and the need to identify an integration path for
each system of interest. Some forms of TI require one to sequentially go from
one ($T,P$) condition starting from the reference state making it impossible
to fully parallelize across ($T,P$) space and thus computationally inefficient
[15]. Hamiltonian Thermodynamic Integration can be parallelized across $T,P$
but requires long simulations where at each step contributions to the full
Hamiltonian are turned on. Wang et al. performed HTI to compare different
phases of lithium using a Deep Potential MLIP [16] where at each step in the
integration path, 300,000 timesteps at 1fs intervals were used. 16 integration
steps were used leading to 4,800,000 timesteps at each condition of
temperature and pressure for each phase.

Alternatively, the Harmonic and Quasi-Harmonic Approximations are commonly
used for solids. This is because the QHA is the simplest non-trivial approach
to including effects of interacting phonons [17, 18]. Coupled with the finite-
displacement approach for predicting phonon spectra, the QHA is independent of
the implementation of interatomic potential making it much easier to implement
with MLIPs. Due to relatively low computational cost and feasibility of
Density Functional Perturbation Theory, it is also used in Ab-Initio studies
but is still costly for most cases [7, 19]. When an MLIP is used, it only
needs to be accurate in predicting the phonon spectrum of relevant phases at
different pressures. The QHA explicitly ignores anharmonic contributions,
except that of thermal expansion and its success hinges on the validity of the
harmonic approximation in the material of interest or cancellation of errors.
For highly anharmonic phases such as $\beta$-Ti or Zr, the QHA fails because
the 0K phonon spectrum from which finite-temperature vibrational properties
are derived has imaginary frequencies [20, 21, 5].

The Temperature Dependent Effective Potential (TDEP) and Covariance of Atomic
Displacements (CAD) methods are effective phonon theories at finite
temperatures [5, 20]. They both use statistics collected from MD simulations
to construct “effective” force constant matrices giving phonon dispersion
relations. The dispersion relations obtained are then used in standard
harmonic formulae for entropy, free energy, heat capacity etc [6, 20]. The
difference between the two methods is in how the force constant matrix is
obtained. In TDEP, corrections to the force constant matrix of increasing
order are obtained via an optimization procedure that minimizes the difference
between a model Hamiltonian and the MD lattice dynamics. The TDEP method has
been used in studies using AIMD, machine learning and empirical force fields
[6, 22, 8]. In the CAD method, the force constant matrix is obtained by
relating it to a matrix constructed using the covariance of atomic
displacements from their ideal lattice sites. The derivation is briefly
described in Section 4.1 and in more detail by Huang and Widom [20].

There also exist other approaches based on the Kirkwood approximation for
entropy from MD simulations using pair correlation functions [15]. However
some of these methods work best with pair-potentials, making it difficult to
apply MLIPs.

In this work, we benchmark and discuss the combination of MLIPs with the CAD
to provide computationally feasible and accurate predictions of vibrational
properties of solids. In particular, we investigate the entropy, phonon
dispersion and free energy of lithium metal at finite temperatures for which a
suitably accurate MLIP is available. By performing convergence tests and
benchmarks, we discuss the choice of simulation protocols, parameters and
potential pitfalls. These tests are much easier to perform with the cheaper
interatomic potential. We show that the CAD method requires relatively small
systems and significantly shorter simulation times than TI. Looking at the
particular case of solid lithium, we compare CAD predictions with the QHA,
TDEP and experimental results of thermodynamic properties and the lithium
martensitic transition.

##  2 Results

A workflow for computing vibrational properties of solids using CAD is given
in Fig. 1. The first step is to find or develop an appropriate atomistic
potential for the system. In the case of MLIPs, a number of different methods
exist such as generating datasets from scratch using active learning [23, 24]
to fine-tuning pre-trained universal potentials [25, 26]. The choice of MLIP
depends on the accuracy and speed of performing MD. In Fig. 1 , we also show
that only MLIPs and and empirical potentials fall into the desirable window
where we can perform converged CAD calculations within an hour per NVT
simulation on a single high-performance computing node. Among the MLIPs and
empirical potentials considered, the NequIP32 model is the most accurate. It
is important to note that the Root Mean Square Error (RMSE) is not a complete
indicator of accuracy. One should perform a number of different benchmarks on
multiple properties to make sure the MLIP reproduces the DFT results it was
trained on. The physics of the architecture is also important. The NequIP
message passing neural netwrokd architecture implicitly incorporates very
long-range interactions up to dozens of Angstroms generally resulting in
better force predictions which govern the dynamics [27].

Once the potential is determined, NVT simulations can be performed followed by
CAD calculations as post-processing of MD snapshots. The MD and CAD portions
of the workflow can be performed at any desired conditions in parallel.

By studying the phases of solid lithium, we quantify the convergence and
uncertainties in the CAD calculated properties. Lithium is well-described by
the QHA, having a Grunëisen parameter of 0.98 [28] and it has a rich phase
diagram with a temperature-driven martensitic transition from FCC to BCC at
lower pressures [19]. We therefore compare the calculated properties to
available experimental measurements and QHA predictions of thermodynamic
properties.

Figure 1: Schematic workflow for CAD calculations with MLIPs. The inset plot
demonstrates that to run a converged CAD calculation in less than an hour, it
is necessary to use an MLIP such as NequIP or faster potential.

###  2.1 Convergence with system size and sampling time

Figure 2: Convergence of CAD entropy and vibrational free energy with and
without symmeterization for BCC lithium at 300K with respect to a) number of
atoms in unit cell using 10,000 timesteps b) number of MD steps using 432
atoms. Standard errors (too small to be visible) over five different sets of
initial velocities are included as shaded regions. The black dashed line is
the experimental value. Electronic entropy is not included for clarity and is
very small at ${\sim}$0.04$k_{B}$.

Due to the stochastic nature of molecular dynamics simulations, it is
important to quantify the convergence and statistical uncertainties of
calculated properties with system size and simulation time. We determine this
convergence by predicting entropy and vibrational free energy of BCC lithium
at 300K. In Fig. 2, we demonstrate the convergence criterion of the predicted
entropy and free energy to be ${\sim}$20,000 timesteps, equivalent to 40ps
sampled with a timestep of 2fs. To quantify uncertainty, we repeat each
calculation 4 times with different random seeds for the initial velocities and
include the standard deviation of results. The converged simulation time is
significantly shorter than the simulation times needed in methods such as TI.
We also discuss the effect of different sampling strategies such as using a
sampling stride to reduce correlations, sampling later in the trajectory and
varying the timestep in the Supplementary Information. All the strategies
converge similarly to the same values as long as the simulation remains stable
except in the case of phonon dispersion discussed in Section 2.2. We also note
that, for integrated properties such as entropy and free energy, using crystal
symmetry to constrain CAD frequencies has little effect for sufficiently long
simulations but drastically affects the phonon dispersion. The predictions
using symmetry tend to converge slightly faster.

In Fig. 2, we observe that entropy and free energy values converge at
${\sim}$1,000 atoms. We further note that the free energy contributions $Ts$
and $f_{vib}$ at ${\sim}$500 atoms are 0.19meV/atom and 0.04meV/atom from the
converged value respectively. These are small energy differences, well within
typical error bounds of even typical DFT calculations. Therefore, in cases
where high-throughput predictions are performed, it is acceptable to use
${\sim}$500 atoms to save computational effort. This also makes CAD more
accessible to Ab-Initio Molecular Dynamics albeit still at significant cost.

For stable simulations, the predicted entropy and free energy are also robust
as can be seen from the small standard deviation. The standard deviation with
the random seed used to generate the initial velocities is small and barely
visible in Fig. 2. One therefore does not need to perform a large number of
simulations with different initial velocities. Performing more than one
simulation is however still useful as simulations in which instabilities or
phase change occurs can be identified with large variations in entropy values.
An example of diverging entropy and free energy values is shown in Fig. 4,
near the experimental melting point $T_{m}$=453.69K.

The determined convergence criteria mean that the vibrational entropy and free
energy are easily accessible to MLIPs and empirical potential simulations.
Each MD simulation at a given temperature and volume can be completed within
40min on a single NVIDIA A100 GPU using NequIP, one of the most costly MLIPs
and orders of magnitude less for others. We also observe that the convergence
criteria for system size and simulation times correspond similarly to those
observed in the TDEP method [6].

###  2.2 Finite-Temperature Phonon Dispersion Relations

Figure 3: Phonon dispersions and DOS for BCC lithium at 300K predicted using
the CAD and the harmonic approximation at the corresponding equilibrium volume
determined by NPT simulations. Both approaches reproduce experimental results
well.

The formulation of the CAD method in terms of effective harmonic frequencies
naturally allows the prediction of phonon dispersion relations and density of
states (DOS). In Fig. 3, we compare an example of a CAD dispersion and DOS for
BCC lithium at 300K to experimental dispersion measured at the same
temperature. Apparently, the CAD dispersion reproduces the experimental
results very well. It particularly captures the softening of phonon modes as
temperature increases, for example at the $H$ point resulting in a decreasing
maximum frequency with temperature compared corresponding to a smaller lattice
parameter. A traditional harmonic phonon calculation at the predicted
equilibrium volume corresponding to $T$=300K is an alternative approach to
obtaining a “finite-temperature” dispersion. The harmonic dispersion also
corresponds very well with experiments. This highlights the consistency
between the two different approaches to obtaining the dispersion relations
albeit the harmonic calculation being significantly cheaper since it does not
require MD simulations.

When predicting the CAD dispersion and DOS, there are a few notes of caution.
The CAD dispersion converges after slightly longer simulation time than
entropy, but not necessarily needing more timesteps. It is usually better to
run a simulation with larger timesteps, even as much as 5fs. Additionally,
symmeterization of the covariance matrix is critically important for
predicting a reasonable CAD dispersion even if the entropy and free energy are
not as affected. Instabilities in the simulation can lead to “oversoftening”
of some phonon modes and even imaginary frequencies. The physical significance
of these imaginary frequencies should be considered with care as the
frequencies in CAD are empirical. Multiple repeated runs with different
initial velocities and larger supercells give a clearer picture of the degree
of variation in the CAD dispersion and DOS.

###  2.3 Thermodynamics at Zero Pressure

Figure 4: Thermodynamic properties as a function of temperature, namely a)
Lattice constant b) Electronic entropy c) Total Entropy of FCC and BCC,
including electronic entropy, d) Residual of entropy with respect to NIST
experimental data. All methods are in good agreement with experimental data
and each other below 100K, capturing the small discontinuity $\Delta
s$=0.067$k_{B}$ within 5%. Above 100K, CAD, TDEP and QHA do well and are well
within experimental error but diverge in different directions. e) Gibbs free
energy referenced at $T=$298.15K. CAD and TDEP predictions are slightly
different from QHA and experiment which agree well. f) Difference between the
FCC and BCC Gibbs free energy. All methods predict $T_{c}{\sim}$150K compared
to 80K from experiment. The shaded area is the standard error for each phase
in CAD.

We show some thermodynamic properties of lithium metal as a function of
temperature at zero pressure, effectively corresponding to room pressure in
Fig. 4. For CAD, each calculation is performed at the corresponding
equilibrium volume at the given temperature found by averaging the volume in
an NPT simulation. Therefore thermal expansion due to classical dynamics is
accounted for.

In Fig. 4, we show that the lattice constant as predicted by the volume that
minimizes the CAD free energy is most consistent with experimental data. The
equilibrium volume as predicted by NPT simulations in MD also have matching
qualitative behavior. The CAD lattice constant is closer to experimental
values and always higher than the MD lattice constant as it accounts for the
zero-point volume correction. The difference between them however is small
when considering it’s effect on predicted entropies and free energies. The QHA
predicted lattice constant deviates significantly more from the trend in
experimental data despite matching the experimental room temperature lattice
constant very well.

In Fig. 4 we show the entropy as predicted by CAD, TDEP and QHA compared to
experiment including the small linear correction of electronic entropy shown
in the Supplementary Information. At low temperatures, it is expected that the
QHA performs well as it is a low temperature approximation. CAD and TDEP
perform just as well at low temperatures as the dynamics converge to becoming
more harmonic in the low-temperature regime. Figure 4 also shows the entropy
of BCC lithium, including electronic entropy up to the melting point,
$T_{m}$=453.69K. Both CAD and QHA describe the system within experimental
error bounds (${\sim}0.05k_{B}$) of the NIST data up to the melting point
[31]. We also highlight that there is significant variation in experimental
heat capacity ($C_{p}$) measurements near the melting point as shown in Fig.
4. The NIST data is fit to the most recent Douglas et al. [32] experiments.
The Douglas et al. results are the lowest values measured, all other
experiments measured higher heat capacity by up to 30% which would result in
higher entropy above 300K. More experiments with modern techniques would be
useful to investigate the discrepancy between different experiments. CAD and
TDEP generally estimate higher entropy values at higher temperatures than QHA
which could explain the higher experimental values of heat capacity and can
potentially be attributed to anharmonicity beyond the QHA. In all properties
shown, there is a clear correspondence between the TDEP and CAD predicted
quantities since they both employ classical (Born-Oppenheimer) molecular
dynamics to sample the configuration space.

In Fig. 4, we show the Gibbs free energy of BCC lithium as a function of
temperature. The CAD and TDEP predictions decrease faster with temperature
whereas the QHA prediction matches very well with experiment. This difference
is attributed to the slightly higher entropy in Fig. 4. In this case, it
appears the QHA performs significantly better. The CAD and TDEP free energy
shown is based on eq. 5 which is a harmonic equation but using CAD
frequencies. Korotaev et al. note that using the harmonic expression for the
free energy is questionable even in the TDEP approach and argues that it is
better to use

$f=u-Ts,$ (1)

where $u$ is the average internal energy from the molecular dynamics
simulation and $s$ is the TDEP entropy. In the case of lithium at low
temperature however, such an approximation leads to a gross violation of the
second law of thermodynamics as shown by the non-monotonic free energy as a
function of temperature in Fig. 4. One potential cause of this issue is in the
inability of the MLIP to perfectly reproduce DFT. Furthermore, the underlying
DFT used to train the MLIP also has an intrinsic error from parameter choices,
choice of exchange-correlation functional and other approximations [33]. These
differences however may not account for why the QHA performs so well for the
free energy.

We believe the most likely source source of error is the assumption of Born-
Oppenheimer dynamics in the low temperature limit for the very light lithium
atoms. The Born-Oppenheimer approximation ignores nuclear quantum effects such
as coherence [34]. The assumption of Boltzmann distributed velocities and
classical dynamics is itself flawed at low temperatures which are better
described with a full quantum mechanical treatment. For instance, one can
calculate the thermal De Broglie wavelength for lithium at 10K as
${\sim}2\AA{}$, significantly larger than the atomic displacements. The effect
is to significantly overestimate the enthalpy at low temperatures which
amounts by assuming that the equipartition of energy holds which is invalid at
low temperatures [4]. Feng et al. found that nuclear quantum effects reduced
the free energy of the high pressure cI16 phase of lithium by up to 15meV/atom
at nearly 300K using Path-Integral Molecular Dynamics [34]. It would therefore
be useful to either improve or rigorously validate the use of harmonic
expressions using CAD or TDEP frequencies to identify the source of error.
More discussion on this issue is included in the Supplementary Information.

An alternative approach might be to employ Wigner’s quantum correction to
classical statistical mechanics [35]. This is a rapidly convergent
perturbative expansion of the free energy in the ratio of $\hbar\omega/k_{B}T$
[36].

We also note that using eq. 13 apparently gives better agreement with
experiment at higher temperature explaining how it has been successfully used
to model other transitions such as that of Ti which occurs at 1155.15K where a
classical approach might be more appropriate [37].

Free energy calculations are most useful if used to determine phase stability
and phase transition critical temperatures ($T_{c}$). Experimentally, at
$T_{c}$=80($\pm 10$)K, a martensitic transition from the face-centered cubic
(FCC) structure to the body-centered cubic (BCC) structure occurs with
increasing temperature [38].

In Fig. 4 we show the Gibbs free energy difference between FCC and BCC lithium
as a function of temperature. We observe that CAD, TDEP and QHA methods
predict a phase transition at about 150K where the free energies cross,
significantly higher than the experimental value. Predicting $T_{c}$ for this
transition accurately has been a long standing challenge for computational
methods. The free energy differences associated with the transition are of the
order of 1meV/atom. An error in the cohesive energy of 0.5meV/atom can change
$T_{c}$ by up to 100K. The NequIP MLIP has a reported RMSE of 1meV/atom for
energy. The effect of the zero-point correction to the volume, a purely
quantum mechanical effect plays a small role. We estimate that this effect
contributes at most a 0.1meV/atom difference. An additional source of error is
the choice of Exchange-Correlation Functional. We found that while PBE and
PBESol gave very similar results, SCAN gave a cohesive energy about
0.3meV/atom higher.

In Fig. 4, it appears as though the FCC and BCC phases have very similar Gibbs
free energy above $T_{c}$ in the CAD prediction, unlike in the QHA before a
large divergence. The large divergence is due to melting which can be observed
in the atomic configurations and happens for FCC before it happens for BCC.
The divergence and deviation of the CAD entropy may be correlated with the
(meta-)stability of the phase being simulated, even if it remains the dominant
phase in the simulation. This would mean that if one is not careful, it is
possible to overestimate an unstable structure’s entropy. An overestimated
entropy leads to a more negative contribution of the $Ts$ term to the free
energy. Ab-initio calculations by Jerabek et al.[19] and Xie et al. [39]
argued that the FCC and BCC phase are quasi-degenerate with BCC being
metastable below the $T_{c}$ with barriers less than 1meV/atom corresponding
to ${\sim}$12K. The Zero Point Energy (ZPE) and vibrational entropy are
important in stabilizing BCC as it is a temperature-driven transition [40].
With increasing temperature, the transition from FCC to BCC becomes “almost
barrierless” hence FCC is not metastable above $T_{c}$. This conclusion aligns
with the observation in Fig. 4.

These subtle caveats imply that confident use of the CAD method requires one
to already know the stable phase. There are however some simple checks one can
make to test the entropy prediction such as checking if the magnitude of the
entropy is reasonable compared to other phases at the same temperature, if the
dependence of free energy with temperature is strictly decreasing and more
weakly, if the entropy decreases with increasing pressure if the material does
not have negative thermal expansion [41]. These constraints are placed by the
second law of thermodynamics and are discussed in Appendix A. This also means
the CAD entropy is a potential descriptor for meta-stability since divergence
of the entropy suggests there exists a possible phase transition that may not
be apparent in structural analysis or requires too long a simulation time to
occur.

In the entropy, free energy and phase transition predictions, the CAD and TDEP
approaches have very similar results. This suggests that the limitations
observed are likely limitations to either the grounding of the methods in the
harmonic theory, or the use of classical molecular dynamics simulations for a
strongly quantum system. And so far, the two methods appear interchangeable
for the case of lithium metal.

##  3 Discussion

The CAD method overcomes some of the computational challenges of TI and in
combination with MLIPs, is significantly computationally cheaper. Firstly, CAD
properties (vibrational entropy, free energy etc.) can be spot-evaluated at a
particular temperature and pressure for a particular phase. This means the
evaluations are completely parallelizable in $T,P$ space as with HTI and TDEP.
However, one only needs ${\sim}20,000$ samples from a single NVT simulation,
${{\sim}}240$ times less timesteps at each $T,P$ evaluation compared to HTI.
Perhaps most importantly, since CAD is implemented as a post-processing step
to an NVT/NPT simulation it is significantly easier to use. The researcher
does not need to spend significant effort determining parameters like
integration steps, and reference states. Unlike TI however, the CAD method
only accounts for vibrational properties. Liquid, configurational entropy and
free energy would have to be treated within a different framework.

Compared to CAD, the QHA is significantly computationally cheaper when using a
MLIP. This is because in the QHA, only the calculation of the phonon
dispersion requires the interatomic potential at each pressure. Subsequent
steps are mostly post-processing using standard analytical formulae for the
temperature dependence. The CAD method, like TDEP requires MD simulations at
each temperature. The CAD method however does not fail for anharmonic
materials at 0K unlike the QHA. Huang and Widom showed that the entropy and
phase transition of $\beta$-Ti can be predicted using the CAD method.
Additionally, CAD implicitly accounts for anharmonicity as it is an effective
harmonic theory, with the expectation of better results at higher temperatures
compared to QHA. As lithium does not have strong anharmonicity, studies on
other materials would have to be carried out to validate this assertion.

##  4 Methods

In this section, we describe the basics of the theory behind the CAD
calculations, model architectures, data generation and simulation details.

###  4.1 Covariance of Atomic Displacements

The CAD method introduced by Huang and Widom[20] allows the prediction of
vibrational properties such as entropy and free energy of solids at finite
temperatures, and phonon dispersion including anharmonicity. The vibrational
properties are derived from a covariance matrix
($\mathbf{\Sigma}_{\mathcal{U}}$) constructed from molecular dynamics or
Monte-Carlo simulations. Each matrix element of
$\mathbf{\Sigma}_{\mathcal{U}}$ is given by the covariance between the
displacement ($x_{i\mu}$) of atom $i$ in each coordinate $\mu$ and the
displacement ($x_{j\nu}$) of every other atom $j$ in each coordinate $\nu$
relative to their ideal lattice positions over the course of a trajectory.
$\mathbf{\Sigma}_{\mathcal{U}}$ is therefore a $3N\times 3N$ matrix with
elements $\Sigma_{i\mu j\nu}=\langle x_{i\mu}x_{j\nu}\rangle$. One can relate
the mass-weighted covariance matrix,
$\tilde{\Sigma}_{i\mu,j\nu}=\Sigma_{i\mu,j\nu}\sqrt{m_{i}m_{j}}$ to the mass-
reduced force constant matrix
$\tilde{C}_{i\mu,j\nu}=C_{i\mu,j\nu}/\sqrt{m_{i}m_{j}}$ using the relation

$\tilde{C}=\frac{1}{\beta}\tilde{\Sigma}^{-1}.$ (2)

Equation 2 identifies the eigenvalues of $\tilde{C}$ as inverse frequencies
$\omega_{k\mu}^{-1}$ which share eigenvectors with $\tilde{\Sigma}$. For a
harmonic potential eq. 2 is exact whereas for a general potential, it is an
effective force constant matrix from which to extract effective vibrational
frequencies that include anharmonic effects. The frequencies are effective in
that they implicitly account for anharmonicity beyond the quasiharmonic
approximation which only accounts for anharmonicity due to thermal expansion.
The expectation is therefore that the CAD frequencies represent higher
temperature properties more accurately.

Given the set of CAD frequencies, it is therefore possible to use the standard
harmonic expressions for thermodynamic properties. The vibrational entropy at
temperature $T$ is given by

$s_{vib}=k_{B}\sum\limits_{k\mu}\left[\frac{\beta\hbar\omega_{k\mu}}{e^{\beta\hbar\omega_{k\mu}}-1}-\ln(1-e^{-\beta\hbar\omega_{k\mu}})\right],$
(3)

where $k_{B}$ is the Boltzmann constant, $\beta=1/k_{B}T$ and $\hbar$ is the
reduced Planck’s constant, $e$ is the exponential factor and $\Lambda$ is the
thermal de Broglie wavelength of the material.

An alternative expression for the entropy of a single atomic species can be
written succinctly in terms of the covariance matrix as

$s^{classical}_{vib}=\frac{1}{2}k_{B}\ln\left[(2\pi
e/\Lambda^{2})^{3N}\det{(\mathbf{\Sigma}_{\mathcal{U}})}\right]+\frac{3}{2}k_{B},$
(4)

a form that does not explicitly invoke vibrational frequencies or the harmonic
approximation. The expression in eq. 4 is termed the classical entropy as it
is derived from information-theoretic principles but diverges to $-\infty$ as
$T\rightarrow 0$. The third law of thermodynamics demands that the entropy
should vanish instead. We therefore use eq. 3 unless otherwise stated as it
more faithfully reproduces experimental results over the entire temperature
range as shown in Fig. 2 and does not erroneously diverge.

The vibrational Helmholtz free energy can be expressed as

$f_{vib}=\frac{1}{N}\sum\limits_{k\mu}\left[\frac{1}{2}\hbar\omega_{k\mu}+\frac{1}{\beta}\ln\left(1-e^{-\beta\hbar\omega_{k\mu}}\right)\right].$
(5)

Other properties such as vibrational heat capacity, thermal conductivity, can
also be derived using standard formulae. Additionally, given the frequencies
and eigenvectors of the force constant matrix, one can also construct finite-
temperature phonon dispersion curves and the density of states (DOS).

The total free energy per atom ($f(v,T)$) of a single element system depends
on specific volume ($v$) and is given by

$f(v,T)=u(v,T=0)+f_{vib}(v,T)+f_{el}(v,T).$ (6)

$u(v,T=0)$ is the potential energy per atom of the ideal crystal structure. We
reference all potential energies relative to that of FCC lithium in
equilibrium at 0K. $u$ is obtained from the chosen interatomic potential which
in this case is the NequIP MLIP. It is also reasonable to use the DFT
calculated $u$ for the static cell since it is cheap to calculate and more
reliable. $f_{el}$ is the electronic free energy which can be approximated
from DFT using Fermi-Dirac smearing at the corresponding temperature.

We are usually interested in the stability of materials under conditions of
pressure and temperature hence the Gibbs free energy ($g(p,T)$) or similarly
the NPT ensemble is the more appropriate thermodynamic potential. $g$ can be
calculated from $f$ in the canonical NVT ensemble or enthalpy ($h$) in the NPS
ensemble via Legendre transforms given by

$g(P,T)=f(v,T)+Pv=h(P,s)-Ts$ (7)

We work under the assumption that we can perform simulations in the more
convenient canonical NVT ensemble and use eq. 7 to get $g$. This assumption is
exact in the thermodynamic limit and is already implicitly assumed in MD when
we approximate the true thermodynamic quantities as the mean values from
simulations [42]. We demonstrate the consistency of using the NVT simulation
in place of the NPT simulation in Fig. 4 which shows that the equilibrium
volume obtained using an NPT simulation corresponds well with the volume that
minimizes $f(v,T)$ obtained from multiple NVT simulations for a given $T$.
Expressed in terms of values we extract from simulation,

$g(P,T)=u(v)+f_{vib}(v,T)+f_{el}(T).$ (8)

The most stable structure out of a set of candidates at given temperature and
pressure is the one that minimizes eq. 8.

###  4.2 NequIP Potential

The MLIP used in this work was developed by Phuthi et al. and details of the
dataset and architecture are available [23]. This MLIP was shown to have high
accuracy for describing a number of properties of lithium metal. The NequIP
architecture which uses equivariant features of rotation order up to 2 is
known to be very data efficient and performs with state-of-the-art accuracy on
a number of benchmarks [23, 27].

###  4.3 Simulation Parameters

All DFT, MD and CAD simulations are managed within our open-source simulation
framework, Atomic Simulation Tools (ASIMTools). The framework allows the
building of simulation protocols incorporating multiple different software
tools, allowing reproducible simulations.

For DFT calculations, parameters were converged to less than 1meV/atom in
energy. DFT calculations were done using Quantum Espresso [43] within the
Generalized Gradient Approximation using the Perdew-Burke-Eizenhoff exchange
correlation functional [44] and the Projector Augmented Wave approach [45]
with a plane wave cutoff energy of 1360eV. We used the pseudopotential Li.pbe-
s-kjpaw_psl.1.0.0.UPF from http://www.quantum-espresso.org. A uniform
Brillouin Zone spacing of 0.02$\AA{}^{-1}$ with a Monkhorst-Pack [46] sampling
procedure was used. To help with convergence of the the Fermi surface,
Methfessel-Paxton [47] smearing using a smearing width of 0.27eV was chosen.

For TDEP calculations, the same trajectories were used as in the CAD
calculations and only up to second-order force constants were included.

###  4.4 Molecular Dynamics details

Molecular Dynamics simulations were performed in the Large-scale
Atomic/Molecular Massively Parallel Simulator (LAMMPS) [48]. Unless otherwise
stated, simulations were performed using a timestep of 2fs. Isobaric (NPT) or
isothermal (NVT) simulations were performed using Nosé-Hoover style barostat
and thermostat equations of motion by Shinoda et al. [49]. For CAD runs,
20,000 timesteps were used unless stated otherwise. For FCC, supercells with
500 atoms and for BCC, supercells with 432 atoms were used unless stated
otherwise.

Supplementary information The provided supplementary information provides
additional benchmarks for the CAD method and sampling strategies as well as
more details of the simulations performed.

Acknowledgements This work was supported in part by Oracle Cloud credits and
related resources provided by Oracle for Research. MW and YH were supported by
DOE grant DE-SC0014506.

## Declarations

###  4.5 Data and Code Availability

The MLIP and training data used in this work is available on Zenodo [50]. The
code for the CAD post-processing is available from
https://alloy.phys.cmu.edu/software/CAD.html. An interface for running
benchmarks, MD (using LAMMPS or ASE) and CAD altogether is available as an
extension to ASIMTools and will be made available on github. For all TDEP
calculations, the open-source TDEP software by Knoop et al. was used [6, 51].

###  4.6 Conflicts of Interest

The authors declare no conflicts of interest.

###  4.7 Author Contribution

Conceptualization, M.K.P; Methodology - Modeling - M.K.P., Y.H., M.W.;
Investigation, M.K.P.; Formal Analysis, M.K.P., M.W.; Data Curation, M.K.P. ;
Writing - Original Draft, M.K.P., M.W. and V.V.; Writing - Review & Editing -
M.K.P.; Supervision. M.W. and V.V.

## Appendix A

We briefly discuss why the second law of thermodynamics imposes that the Gibbs
free energy at constant pressure must be a decreasing function of temperature.
One can write out the differential form of the Gibbs free energy

$dG=-SdT+Vdp+\mu dN$ (9)

where

$\frac{\partial G}{\partial T}\big{|}_{p,N}=-S<0.$ (10)

Since the entropy $S$ is positive, $G$ can only be a decreasing function of
$T$.

Now consider the dependence of entropy on pressure at constant temperature
i.e.

$\frac{\partial S}{\partial p}\bigg{|}_{T,N}=-\frac{\partial V}{\partial
T}\bigg{|}_{p,N}$ (11)

The right-hand side is derived using a Maxwell relation and is related to the
coefficient of thermal expansion ($\alpha$) by

$\alpha=\frac{1}{V}\left(\frac{\partial V}{\partial T}\bigg{|}_{p}\right).$
(12)

For most materials, the thermal expansion coefficient is positive but there
are a few exceptions such as pure silicon. One can therefore not rule out that
a materials with increasing entropy as a function of pressure is unstable but
it is likely.

## References

  *   * Jain et al. [2016] Jain, A., Shin, Y., Persson, K.A.: Computational predictions of energy materials using density functional theory. Nature Reviews Materials 1(1), 1–13 (2016) https://doi.org/10.1038/natrevmats.2015.4 . Publisher: Nature Publishing Group. Accessed 2024-06-03
  * Chen et al. [2020] Chen, C., Zuo, Y., Ye, W., Li, X., Deng, Z., Ong, S.P.: A Critical Review of Machine Learning of Energy Materials. Advanced Energy Materials 10(8), 1903242 (2020) https://doi.org/10.1002/aenm.201903242 . Accessed 2024-06-03
  * T. Szczypiński et al. [2021] T. Szczypiński, F., Bennett, S., E. Jelfs, K.: Can we predict materials that can be synthesised? Chemical Science 12(3), 830–840 (2021) https://doi.org/10.1039/D0SC04321D . Publisher: Royal Society of Chemistry. Accessed 2024-06-03
  * Frenkel and Smit [2001] Frenkel, D., Smit, B.: Understanding Molecular Simulation: From Algorithms to Applications, 2nd edition edn. Academic Press, San Diego (2001)
  * Hellman et al. [2011] Hellman, O., Abrikosov, I.A., Simak, S.I.: Lattice dynamics of anharmonic solids from first principles. Physical Review B 84(18), 180301 (2011) https://doi.org/10.1103/PhysRevB.84.180301 . Publisher: American Physical Society. Accessed 2024-05-15
  * Hellman et al. [2013] Hellman, O., Steneteg, P., Abrikosov, I.A., Simak, S.I.: Temperature dependent effective potential method for accurate free energy calculations of solids. Physical Review B 87(10), 104111 (2013) https://doi.org/10.1103/PhysRevB.87.104111 . Accessed 2024-03-27
  * Korotaev et al. [2018] Korotaev, P., Belov, M., Yanilkin, A.: Reproducibility of vibrational free energy by different methods. Computational Materials Science 150, 47–53 (2018) https://doi.org/10.1016/j.commatsci.2018.03.057 . Accessed 2024-05-04
  * Grabowski et al. [2019] Grabowski, B., Ikeda, Y., Srinivasan, P., Körmann, F., Freysoldt, C., Duff, A.I., Shapeev, A., Neugebauer, J.: Ab initio vibrational free energies including anharmonicity for multicomponent alloys. npj Computational Materials 5(1), 1–6 (2019) https://doi.org/10.1038/s41524-019-0218-8 . Number: 1 Publisher: Nature Publishing Group. Accessed 2023-11-06
  * Behler and Csányi [2021] Behler, J., Csányi, G.: Machine learning potentials for extended systems: a perspective. The European Physical Journal B 94(7), 142 (2021) https://doi.org/10.1140/epjb/s10051-021-00156-1 . Accessed 2024-05-13
  * Martin [2008] Martin, R.M.: Electronic Structure: Basic Theory and Practical Methods, 1st edition edn. Cambridge University Press, Cambridge, UK ; New York (2008)
  * Kresse and Furthmüller [1996] Kresse, G., Furthmüller, J.: Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Computational Materials Science 6(1), 15–50 (1996) https://doi.org/10.1016/0927-0256(96)00008-0 . Accessed 2023-04-03
  * Kocer et al. [2022] Kocer, E., Ko, T.W., Behler, J.: Neural Network Potentials: A Concise Overview of Methods. Annual Review of Physical Chemistry 73(Volume 73, 2022), 163–186 (2022) https://doi.org/10.1146/annurev-physchem-082720-034254 . Publisher: Annual Reviews. Accessed 2024-05-13
  * Liu et al. [2023] Liu, Y., He, X., Mo, Y.: Discrepancies and error evaluation metrics for machine learning interatomic potentials. npj Computational Materials 9(1), 1–13 (2023) https://doi.org/10.1038/s41524-023-01123-3 . Publisher: Nature Publishing Group. Accessed 2024-06-03
  * Vita and Schwalbe-Koda [2023] Vita, J.A., Schwalbe-Koda, D.: Data efficiency and extrapolation trends in neural network interatomic potentials. Machine Learning: Science and Technology 4(3), 035031 (2023) https://doi.org/10.1088/2632-2153/acf115 . Publisher: IOP Publishing. Accessed 2024-06-03
  * Sluss et al. [2022] Sluss, C.C., Pittman, J., Nicholson, D.M., Keffer, D.J.: Exploration of Entropy Pair Functional Theory. Entropy 24(5), 603 (2022) https://doi.org/10.3390/e24050603 . Number: 5 Publisher: Multidisciplinary Digital Publishing Institute. Accessed 2024-06-12
  * Wang et al. [2023] Wang, X., Wang, Z., Gao, P., Zhang, C., Lv, J., Wang, H., Liu, H., Wang, Y., Ma, Y.: Data-driven prediction of complex crystal structures of dense lithium. Nature Communications 14(1), 2924 (2023) https://doi.org/10.1038/s41467-023-38650-y . Number: 1 Publisher: Nature Publishing Group. Accessed 2023-06-20
  * Fultz [2010] Fultz, B.: Vibrational thermodynamics of materials. Progress in Materials Science 55(4), 247–352 (2010) https://doi.org/10.1016/j.pmatsci.2009.05.002 . Accessed 2024-05-18
  * Allen [2020] Allen, P.B.: Theory of thermal expansion: Quasi-harmonic approximation and corrections from quasi-particle renormalization. Modern Physics Letters B 34(02), 2050025 (2020) https://doi.org/10.1142/S0217984920500256 . Publisher: World Scientific Publishing Co. Accessed 2024-06-03
  * Jerabek et al. [2022] Jerabek, P., Burrows, A., Schwerdtfeger, P.: Solving a problem with a single parameter: A smooth bcc to fcc phase transition for metallic lithium. Chemical Communications 58(96), 13369–13372 (2022) https://doi.org/10.1039/D2CC04928G . arXiv:2211.06317 [cond-mat]. Accessed 2024-04-19
  * Huang and Widom [2022] Huang, Y., Widom, M.: Vibrational Entropy of Crystalline Solids from Covariance of Atomic Displacements. Entropy 24(5), 618 (2022) https://doi.org/10.3390/e24050618 . Number: 5 Publisher: Multidisciplinary Digital Publishing Institute. Accessed 2023-04-25
  * Kadkhodaei et al. [2017] Kadkhodaei, S., Hong, Q.-J., Van De Walle, A.: Free energy calculation of mechanically unstable but dynamically stabilized bcc titanium. Physical Review B 95(6), 064101 (2017) https://doi.org/10.1103/PhysRevB.95.064101 . Accessed 2024-06-11
  * Romero et al. [2015] Romero, A.H., Gross, E.K.U., Verstraete, M.J., Hellman, O.: Thermal conductivity in PbTe from first principles. Physical Review B 91(21), 214310 (2015) https://doi.org/10.1103/PhysRevB.91.214310 . Accessed 2024-05-22
  * Phuthi et al. [2024] Phuthi, M.K., Yao, A.M., Batzner, S., Musaelian, A., Guan, P., Kozinsky, B., Cubuk, E.D., Viswanathan, V.: Accurate Surface and Finite-Temperature Bulk Properties of Lithium Metal at Large Scales Using Machine Learning Interaction Potentials. ACS Omega 9(9), 10904–10912 (2024) https://doi.org/10.1021/acsomega.3c10014 . Publisher: American Chemical Society. Accessed 2024-03-15
  * Vandermause et al. [2020] Vandermause, J., Torrisi, S.B., Batzner, S., Xie, Y., Sun, L., Kolpak, A.M., Kozinsky, B.: On-the-fly active learning of interpretable Bayesian force fields for atomistic rare events. npj Computational Materials 6(1), 1–11 (2020) https://doi.org/10.1038/s41524-020-0283-z . Publisher: Nature Publishing Group. Accessed 2024-06-08
  * Deng et al. [2023] Deng, B., Zhong, P., Jun, K., Riebesell, J., Han, K., Bartel, C.J., Ceder, G.: CHGNet as a pretrained universal neural network potential for charge-informed atomistic modelling. Nature Machine Intelligence 5(9), 1031–1041 (2023) https://doi.org/10.1038/s42256-023-00716-3 . Publisher: Nature Publishing Group. Accessed 2024-06-08
  * Batatia et al. [2022] Batatia, I., Kovacs, D.P., Simm, G.N.C., Ortner, C., Csanyi, G.: MACE: Higher Order Equivariant Message Passing Neural Networks for Fast and Accurate Force Fields. (2022). https://openreview.net/forum?id=YPpSngE-ZU Accessed 2024-06-08
  * Batzner et al. [2022] Batzner, S., Musaelian, A., Sun, L., Geiger, M., Mailoa, J.P., Kornbluth, M., Molinari, N., Smidt, T.E., Kozinsky, B.: E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nature Communications 13(1), 2453 (2022) https://doi.org/10.1038/s41467-022-29939-5 . Number: 1 Publisher: Nature Publishing Group. Accessed 2022-06-20
  * Chuan-Hui and Xiao-Bao [2001] Chuan-Hui, N., Xiao-Bao, S.: Grüneisen parameter of lithium, sodium and potassium at T=298 K. Journal of Physics and Chemistry of Solids 62(7), 1359–1361 (2001) https://doi.org/10.1016/S0022-3697(01)00007-5 . Accessed 2024-05-18
  * Beg and Nielsen [1976] Beg, M.M., Nielsen, M.: Temperature dependence of lattice dynamics of lithium 7. Physical Review B 14(10), 4266–4273 (1976) https://doi.org/10.1103/PhysRevB.14.4266 . Publisher: American Physical Society. Accessed 2022-05-10
  * Owen and Williams [1954] Owen, E.A., Williams, G.I.: X-Ray Measurements on Lithium at Low Temperatures. Proceedings of the Physical Society. Section A 67(10), 895–900 (1954) https://doi.org/10.1088/0370-1298/67/10/306 . Accessed 2022-11-10
  * Chase [1998] Chase, M.: Nist-janaf thermochemical tables, 4th edition (1998)
  * Douglas et al. [1955] Douglas, T.B., Epstein, L.F., Dever, J.L., Howland, W.H.: Lithium: Heat Content from 0 to 900°, Triple Point and Heat of Fusion, and Thermodynamic Properties of the Solid and Liquid1. Journal of the American Chemical Society 77(8), 2144–2150 (1955) https://doi.org/10.1021/ja01613a031 . Publisher: American Chemical Society. Accessed 2024-06-07
  * Sholl [2009] Sholl, D.S.: Density Functional Theory a Practical Introduction. Wiley, Hoboken, N.J (2009)
  * Markland and Ceriotti [2018] Markland, T.E., Ceriotti, M.: Nuclear quantum effects enter the mainstream. Nature Reviews Chemistry 2(3), 1–14 (2018) https://doi.org/10.1038/s41570-017-0109 . Publisher: Nature Publishing Group. Accessed 2024-05-21
  * Wigner [1932] Wigner, E.: On the Quantum Correction For Thermodynamic Equilibrium. Physical Review 40(5), 749–759 (1932) https://doi.org/10.1103/PhysRev.40.749 . Publisher: American Physical Society. Accessed 2024-06-03
  * Deserno and Turgut [2023] Deserno, M., Turgut, O.T.: Leading quantum correction to the classical free energy. American Journal of Physics 91(11), 923–931 (2023) https://doi.org/10.1119/5.0106687 . Accessed 2024-06-03
  * Huang et al. [2022] Huang, Y., Widom, M., Gao, M.C.: Ab initio free energies of liquid metal alloys: Application to the phase diagrams of Li-Na and K-Na. Physical Review Materials 6(1), 013802 (2022) https://doi.org/10.1103/PhysRevMaterials.6.013802 . Accessed 2023-11-16
  * Hultgren [1963] Hultgren, R.: Selected Values of Thermodynamic Properties of Metals And Alloys. Wiley, ??? (1963). Google-Books-ID: AjFRAAAAMAAJ
  * Xie et al. [2008] Xie, Y., Ma, Y.M., Cui, T., Li, Y., Qiu, J., Zou, G.T.: Origin of bcc to fcc phase transition under pressure in alkali metals. New Journal of Physics 10(6), 063022 (2008) https://doi.org/10.1088/1367-2630/10/6/063022 . Accessed 2024-04-19
  * Schaeffer et al. [2015] Schaeffer, A.M., Cai, W., Olejnik, E., Molaison, J.J., Sinogeikin, S., Santos, A.M., Deemyad, S.: Boundaries for martensitic transition of 7Li under pressure. Nature Communications 6(1), 8030 (2015) https://doi.org/10.1038/ncomms9030 . Number: 1 Publisher: Nature Publishing Group. Accessed 2023-09-27
  * Liu et al. [2022] Liu, Z.-K., Wang, Y., Shang, S.-L.: Zentropy Theory for Positive and Negative Thermal Expansion. Journal of Phase Equilibria and Diffusion 43(6), 598–605 (2022) https://doi.org/10.1007/s11669-022-00942-z . Accessed 2024-06-12
  * Frenkel and Smit [2002] Appendix H - Statistical Mechanics of the Gibbs “Ensemble”. In: Frenkel, D., Smit, B. (eds.) Understanding Molecular Simulation (Second Edition), pp. 563–571. Academic Press, San Diego (2002). https://doi.org/10.1016/B978-012267351-1/50027-4 . https://www.sciencedirect.com/science/article/pii/B9780122673511500274 Accessed 2024-05-18
  * Giannozzi et al. [2009] Giannozzi, P., Baroni, S., Bonini, N., Calandra, M., Car, R., Cavazzoni, C., Ceresoli, D., Chiarotti, G.L., Cococcioni, M., Dabo, I., Corso, A.D., Gironcoli, S.d., Fabris, S., Fratesi, G., Gebauer, R., Gerstmann, U., Gougoussis, C., Kokalj, A., Lazzeri, M., Martin-Samos, L., Marzari, N., Mauri, F., Mazzarello, R., Paolini, S., Pasquarello, A., Paulatto, L., Sbraccia, C., Scandolo, S., Sclauzero, G., Seitsonen, A.P., Smogunov, A., Umari, P., Wentzcovitch, R.M.: QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. Journal of Physics: Condensed Matter 21(39), 395502 (2009) https://doi.org/10.1088/0953-8984/21/39/395502 . Accessed 2023-02-27
  * Perdew et al. [1996] Perdew, J.P., Burke, K., Ernzerhof, M.: Generalized Gradient Approximation Made Simple. Physical Review Letters 77(18), 3865–3868 (1996) https://doi.org/10.1103/PhysRevLett.77.3865 . Accessed 2019-12-03
  * Blöchl [1994] Blöchl, P.E.: Projector augmented-wave method. Physical Review B 50(24), 17953–17979 (1994) https://doi.org/10.1103/PhysRevB.50.17953 . Publisher: American Physical Society. Accessed 2023-02-27
  * Monkhorst and Pack [1976] Monkhorst, H.J., Pack, J.D.: Special points for Brillouin-zone integrations. Physical Review B 13(12), 5188–5192 (1976) https://doi.org/10.1103/PhysRevB.13.5188 . Accessed 2021-02-16
  * Methfessel and Paxton [1989] Methfessel, M., Paxton, A.T.: High-precision sampling for Brillouin-zone integration in metals. Physical Review B 40(6), 3616–3621 (1989) https://doi.org/10.1103/PhysRevB.40.3616 . Publisher: American Physical Society. Accessed 2023-02-27
  * Thompson et al. [2022] Thompson, A.P., Aktulga, H.M., Berger, R., Bolintineanu, D.S., Brown, W.M., Crozier, P.S., Veld, P.J.i.t., Kohlmeyer, A., Moore, S.G., Nguyen, T.D., Shan, R., Stevens, M.J., Tranchida, J., Trott, C., Plimpton, S.J.: LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Comp. Phys. Comm. 271, 108171 (2022) https://doi.org/%****␣sn-article.bbl␣Line␣875␣****10.1016/j.cpc.2021.108171
  * Shinoda et al. [2004] Shinoda, W., Shiga, M., Mikami, M.: Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Physical Review B 69(13), 134103 (2004) https://doi.org/10.1103/PhysRevB.69.134103 . Publisher: American Physical Society. Accessed 2024-01-10
  * Mgcini [2024] Mgcini, P.: Dataset and Potentials for Lithium Metal, Phuthi M. K. et al. (2024). Zenodo (2024). https://doi.org/10.5281/zenodo.10470793 . https://zenodo.org/records/10470793 Accessed 2024-05-18
  * Knoop et al. [2024] Knoop, F., Shulumba, N., Castellano, A., Batista, J.P.A., Farris, R., Verstraete, M.J., Heine, M., Broido, D., Kim, D.S., Klarbring, J., Abrikosov, I.A., Simak, S.I., Hellman, O.: TDEP: Temperature Dependent Effective Potentials. Journal of Open Source Software 9(94), 6150 (2024) https://doi.org/10.21105/joss.06150 . Accessed 2024-06-01

## Supplementary Information for: Vibrational Entropy and Free Energy of Solid
Lithium using Covariance of Atomic Displacements Enabled by Machine Learning

### Consistency of CAD thermodynamics with Molecular Dynamics

In Fig. 5, we show the energy-volume curves for the static potential energy
($u$) which is typically used as the lattice constant prediction at 0K using
DFT. The lattice constant is usually determined by finding the minimum of the
curve and is shown as an orange vertical line. In addition we show the
Helmholtz free energy predicted by CAD as a function of volume and the
corresponding minimum which includes the effect of thermal expansion and the
zero-point volume correction. The lattice constant predicted using the minimum
of the Helmholtz free energy is therefore always higher than that of the
internal energy. One can also predict the lattice constant using NPT
simulations in MD and that gives a result in between the other cases as it
accounts for expansion due to kinetic energy but not the zero-point volume
correction.

Figure 5: Energy-volume curves for the static internal energy ($u(v,T=0)$) and
CAD predicted Helmholtz free energy ($f(v,T)$). The CAD predicted equilibrium
volume is slightly higher but consistent with the equilibrium volume found
using NPT simulations in MD ($v0$)

We can also do a comparison of the individual contribution of the $pV$ and
$Ts$ terms in the free energies as shown in Fig. 6. Since the equilibrium
volumes from NPT simulations were used for the NVT simulations in CAD, the
pressure and thus $pV$ are very close to zero, much smaller than the
contribution of $Ts$ term from the entropy.

Figure 6: Relative contribution of $pV$ and $Ts$ term at zero pressure as a
function of temperature.

The electronic entropy adds a small linear addition to the vibrational entropy
and is shown in Fig. 7. The entropy of each phase is calculated by applying
Fermi-Dirac smearing to static DFT calculations at the corresponding
equilibrium volume. There is very little difference between the entropy values
hence one would expext that this does not affect phase transitions as much.
More thorough consideration might be necessary for cases where interesting
electronic ground states emerge such as the formation of electrides in lithium
at extreme gigapascal pressures.

Figure 7: Electronic entropy of BCC and FCC lithium as a function of
temperature.

### More details on sampling strategies

We tested a number of different sampling strategies and quantify their effects
for a few examples by predicting the entropy of BCC lithium at 300 K.

In Fig. 8, we compare the convergence of sampling from the initial ideal
crystal in the MD simulation with the case where we start sampling after the
5,000th timestep. This could be of significance if the initial configuration
or velocities fluctuate at the beginning of the simulation or if computer
storage for the trajectories is limited. We see that there is no significant
difference in the observed results after sufficient sampling as long as the
BCC structure remains stable throughout the simulation.

Figure 8: Starting sampling later in the simulation leads to faster
convergence. If storage is a limiting factor, then one can only start storing
the coordinates of the trajectory after some burn-in time.

Another important consideration in MD simulations is the timestep chosen in
the simulation. The timestep is typically chosen to match with the fastest
mode expected in the simulation [42] on the order of 1fs. In Fig. 9 we show
that the timestep has little effect for a small range of investigated
timesteps with differences on the order of 0.02$k_{B}$. It therefore seems
advantageous to use larger timesteps to sample more effectively and we found
that this can sometimes lead to faster convergence of sensitive properties
like the phonon dispersion. The effect of timestep at different temperatures
would have to be considered further.

Figure 9: Interestingly, the choice of timestep has very little effect on the
predicted entropy with all predictions being within the small standard
deviation of results. Even for relatively large timesteps of 0.005fs. For
faster convergence one can therefore use a larger timestep.

A common strategy for improving sampling is to only take samples for post-
processing after every N steps where N is the stride. This makes the samples
less correlated therefore leading to better statistical averages. We compare
the effect of different N with the 1,000 of snapshots used in the CAD
calculation. We see that the stride also does not have a significant effect by
noting that the poor convergence at lower strides is actually due to not
sampling for long enough. The convergence of overall length of the simulation
is more significant which occurs at about 40ps.

Figure 10: Convergence of the entropy with different choice of sampling stride
using 1,000 timesteps for each point.

### Inconsistency of Classical Free Energy with the Second Law of
Thermodynamics at Low Temperature

As stated in the text, one might naïvely expect that calculating the Gibbs
free energy using

$g=u-Ts+pV$ (13)

is valid at all temperatures where $u$ is the internal energy from MD
averaging. This however leads to a gross violation of the second law of
thermodynamics. This is attributed to the failure of classical Born-
Oppenheimer dynamics at lower temperatures, more so for light elements where
there is expected to be significant coherence and nuclear quantum effects. It
is there only valid to use the quantum mechanical free energy contribution in
the low temperature regime. The cause is due to the overstimation of the
enthalpy as shown in Fig. 12. The Boltzmann velocity distribution from a
classical simulation does not correspond to the more coherent quantum
mechanical dynamics which is dominated by the long wavelength phonon modes. In
other words, the equipartition theorem does not hold.

Figure 11: Using eq. 13 leads to a violation of the second law of
thermodynamics as g is supposed to be a strictly decreasing function of T.
Figure 12: Plot of the energy contributions to the FCC enthalpy at low
temperatures. The experimental enthalpy derived from integrating the heat
capacity is significantly lower than the MD predicted contributions, largely
due to quantum effects.


# Obituary, Professor John D. Barrow 1952–2020111Published in Astronomy &
Geophysics, Volume 61, Issue 6, December 2020, Pages 6.11–6.12,
https://doi.org/10.1093/astrogeo/ataa081
The Sharpest of Minds

João Magueijo1 and John Webb2

1Theoretical Physics Group, The Blackett Laboratory, Imperial College, Prince
Consort Rd., London, SW7 2BZ, UK.
2Clare Hall, University of Cambridge, Herschel Rd, Cambridge CB3 9AL, UK.

On Saturday 26 September, around 4am, John Barrow died aged 67, with his wife
Elizabeth and son Roger at his side. From a scientific perspective, it is hard
to conceive a more premature end. During lockdown alone, whilst undergoing
chemotherapy and in the full knowledge that his cancer was inoperable, John
managed to co-author 11 scientific papers and write a new book (“One Plus
One”). Even by his own standards of productivity, this is staggeringly
impressive, an achievement he was openly proud of. From a broader perspective,
with a wife, 3 children and 5 young grandchildren, many strong friendships,
and so much more to offer the world, he departed far too soon.

Figure 1: John David Barrow (born 1952), British cosmologist, theoretical
physicist, and mathematician. Barrow was a Professor of Mathematical Sciences
at the University of Cambridge since 1999. He was elected Fellow of the Royal
Society in 2003 and was awarded the Faraday Prize in 2008 for his work on
science popularisation. Photographed in 2007 in Barrow’s office at the Centre
for Mathematical Sciences, Cambridge. Photograph provided by Simon Stone,
Science Photo Library, London.

John David Barrow was born on 29 November 1952 to parents Walter and Lois. He
had one sibling, an elder sister, Brenda, who pre-deceased John earlier this
year. He grew up in Wembley, London, attending Barham Primary School and then
Ealing Grammar School for Boys from 1964-71. Whilst his interest in science
was sparked at the age of 12 when he was given a chemistry set, it must have
seemed to the young John Barrow that a sporting career was far more likely. He
was an outstanding runner, beating Steve Ovett at Withdean Stadium in Brighton
in 1971 to become the English Schools 800m champion (over the same distance
for which Ovett went on to win Gold at the 1980 Moscow Olympics). One of the
authors, unaware of this precedent, was once flabbergasted to see the middle-
aged John sprinting into the distance to catch a train. “One never loses the
fast-twitch muscles”, he commented, laughing about this incident on his
deathbed. But John’s mind was even sharper than his muscles and in the end
prevailed in defining his career.

Figure 2: John beating Steve Ovett in the English Schools 800m, Withdean
Stadium, Brighton, 1971.

Even before finishing school, John had already met his future wife, Elizabeth.
Like many of the baby boomer generation, John and his family were regular
Sunday School attendees, where his teacher was one Mr. East, Elizabeth’s
father. The two children became close long before romance blossomed, although
that had happened by the time university beckoned. Equipped with A-level
Physics and Maths, John went up to Durham University to study Mathematics,
whereas Elizabeth won a place at the Hammersmith Hospital Nursing School.
There ensued regular train journeys over the next few years and when visits
were on hold (at exam times John went into a self-imposed “lockdown”), home-
made food parcels were posted by Elizabeth, usually her delicious cakes. In
1974, John gained 1st Class Honours and was accepted at Magdalen College,
Oxford, considering himself “very lucky” to be taken on as a graduate student
by Dennis Sciama. This was the same year that Brandon Carter coined the term
“Anthropic Principle”, a highly controversial set of ideas that subsequently
inspired John’s remarkable 1986 book, co-authored by Frank Tipler, “The
Cosmological Anthropic Principle”. As influential as Carter’s paper was, it
was clearly Barrow and Tipler’s detailed study that garnered more attention
for anthropic theories. But even though John was a pioneer in the field, he
always showed a disarming lack of dogmatism about it.

John’s brilliance had already emerged as a Durham undergraduate and it
flourished at Oxford. He and Elizabeth were fortunate to find very convenient
College accommodation and the couple married in 1975. Their metamorphosis into
a powerful partnership had begun, with Elizabeth “carrying out the tasks that
John could not and John doing what he did best”, as she put it herself. In
1977 John was awarded his PhD, titled “Non-Uniform Cosmological Models”, in
the same year being awarded no less than Oxford’s Wallace Research Prize, a
Junior Research Lectureship at Christ Church College, Oxford, and a Lindemann
Trust Fellowship. After brief negotiations with Christ Church, John and
Elizabeth spent the first year of John’s Lectureship in the Astronomy
Department at the University of California Berkeley, returning to Oxford for
the remaining 2 years. In 1979 he won (jointly with Bernard Carr) a Gravity
Research Foundation award for the awesomely titled essay “Shear Hell Holes and
Anisotropic Universes”. John then returned to the Physics Department at
Berkeley as a Miller Fellow for a further year.

During these four years at Berkeley and Oxford, John produced more than 30
papers, covering a bewilderingly diverse range of topics, including but not
limited to, anisotropic cosmologies, black holes, big bang nucleosynthesis,
and grand unification. The Berkeley years also saw the birth of Elizabeth and
John’s two sons, David and Roger. In 1981 John was offered a Lectureship in
the Astronomy Centre at the University of Sussex. He remained at Sussex for
the next 18 years, becoming a full professor in 1989. It was also there that
John and Elizabeth completed their family, with the birth of their third
child, Louise, in 1984.

The Sussex years saw John’s writing career skyrocket, turning him into one of
the world’s foremost popularisers of science. His first book, “The Left Hand
of Creation” (in collaboration with Joe Silk) appeared in 1983. John’s
seemingly endless energy resulted in 22 books over the next 37 years, covering
the interplay between maths, cosmology, and the arts. In 1999, John’s
outstanding capacity for science communication attracted him to become the the
Director of the Millennium Mathematics Project at Cambridge University, a
Professor of Mathematical Physics, and a Fellow of Clare Hall. He remained at
Cambridge for the rest of his life. Under John’s guidance, the MMP project
developed into a BBC collaboration, reaching over 1.5 million UK
schoolchildren, parents and teachers.

Figure 3: Some of John’s many books.

John was made a Fellow of the Royal Society in 2003. He wrote over 550
scientific articles, becoming an acknowledged world expert in topics as
disparate as anisotropic models of the universe (the connection with chaos
theory, for example), inflation (intermediate inflation models and the
definition of slow-roll parameters), and the theory of cosmological
singularities (the concept of sudden singularities). He was the undisputed
world leader in theories of “varying constants”, pioneered by Paul Dirac in
the 1930s. In close relation to this, it is possible that John’s greatest
scientific passion was understanding how such a small number of natural laws
and parameters could result in such complex outcomes: biology, consciousness,
even nosy physicists laying claim to understanding the origin and evolution of
the universe. He eagerly and objectively sought to know whether the laws of
nature were programmed into the universe and immutable, or if instead the
universe just “made it up as it went along”. If he saw religious overtones in
this quest, he never voiced it to either of us. John received the Templeton
Prize in 2006 and in 2019 he was appointed to the Pontifical Academy of
Sciences by Pope Francis.

But above all we will miss John Barrow, the person. John had a unique sense of
humour, sometimes mischievously so. He was very open-minded scientifically,
combining a deceivingly conservative demeanour with an irrepressible gleam in
the eye in the face of the wildest speculations. Everything was fair game to
him. And although he was never dogmatic about his science, he was also firm in
the face of irrational rejection. One of the authors keeps a trove of peer-
review correspondence in which John sometimes reveals hilarious nonchalance
and lack of emotion in response to the basest insults from irascible referees.
He was the ultimate democrat. He treated students and senior scientists with
equal respect and interest. In some ways he shone even more impressively in
the company of non-scientists, because of his remarkable ability to converse
logically at any level on almost any topic. He even lectured the British
cabinet on chaos theory, at the request of Mrs. Thatcher.

We will also remember John as a bon-vivant. He was fond of haute cuisine and
was a member of a chocolate tasting club (his sons complain that he used to
hide the best ones). He was an enthusiast of music and theatre, with a
definite artistic vein, so much so that in 2002 he penned his own theatre
script “Infinities”. His play was performed that same year in Milan (in
Italian) and in Valencia (in Spanish), receiving multiple accolades (for
example, the 2002 Italian Play of the Year Award). He loved beautiful Italy
and charming Italy loved him. It is telling that his last holiday, a few weeks
before his demise, was in that country.

We were privileged to know John Barrow, the polymath, not only as an admired
professional collaborator, but also as a steadfast friend. Quite apart from
intellectual excellence and creativity, few people combine generosity,
courtesy, sincerity, and humour to the same degree. His memory retention was
second to none and he was frequently referred to as a “walking encyclopaedia”.
One of us tried numerous times to entice him to one of the local village pub
quizzes, albeit unfairly on the opposition, although he never did. John leaves
behind his wife Elizabeth, to whom he was married for 45 years, as well as 3
children and 5 grand-children. On the day of his death, four year-old Poppy
was seated at the family table in order for Elizabeth to explain that Grandpa
had left forever and gone up to the stars. As the news was broken, she stared
at the table top looking concerned and quickly said “Oh no, Grandma! He’s
forgotten to take his mobile phone with him.”

With or without wireless coverage, John will always be with us.