# Vibrational entropy of crystalline solids from covariance of atomic
displacements

Yang Huang  Physics Department, Carnegie Mellon University.  Michael Widom
Physics Department, Carnegie Mellon University.

###### Abstract

The vibrational entropy of a solid at finite temperature is investigated from
the perspective of information theory. Ab initio molecular dynamics (AIMD)
simulations generate ensembles of atomic configurations at finite temperature
from which we obtain the $N$-body distribution of atomic displacements,
$\rho_{N}$. We calculate the information-theoretic entropy from the
expectation value of $\ln{\rho_{N}}$. At a first level of approximation,
treating individual atomic displacements independently, our method may be
applied using Debye-Waller B-factors, allowing diffraction experiments to
obtain an upper bound on the thermodynamic entropy. At the next level of
approximation we correct the overestimation through inclusion of displacement
covariances. We apply this approach to elemental body-centered cubic sodium
and face-centered cubic aluminum, showing good agreement with experimental
values above the Debye temperatures of the metals. Below the Debye
temperatures we extract an effective vibrational density of states from
eigenvalues of the covariance matrix, and then evaluate the entropy quantum
mechanically, again yielding good agreement with experiment down to low
temperatures. Our method readily generalizes to complex solids, as we
demonstrate for a high entropy alloy. Further, our method applies in cases
where the quasiharmonic approximation fails, as we demonstrate by calculating
the HCP/BCC transition in Ti.

##  I Introduction

The importance of entropy as a component of thermodynamic free energy,
together with the difficulty of its calculation, motivates continuing efforts
seeking improved computational approaches [1, 2, 3, 4, 5, 6]. The entropy is a
function of the state of the system, and is in principle determined by the
instantaneous values of every degree of freedom. Most computational approaches
to entropy calculation do not make explicit use of these values, and instead
apply some form of thermodynamic integration to relate the entropy in the
state of interest to some reference point of known entropy [7, 8, 9]. Our
approach recognizes that the entropy equals, in suitable units, the
information required to fully specify the state of the system. We capture this
information in the form of many-body correlation functions obtained from ab
initio molecular dynamics (AIMD).

Multiple types of excitation contribute to the entropy of a solid. Neglecting
correlations among these, we may approximate the entropy as a sum

$S\approx S^{v}+S^{e}+S^{c}+\cdots$ (1)

where $S^{v}$ arises from atomic vibrations [10], $S^{e}$ includes electronic
excitations, the non-vibrational configurational term $S^{c}$ incorporates
vacancies and chemical species substitutions [11]. The additional terms may
include magnetism and other effects [12]. The present paper primarily
addresses the vibrational contribution, but for comparison with experiment we
must include the electronic entropy. While our initial approach is classical,
and intended for applicability at elevated temperatures close to melting, we
show how quantum effects can be incorporated to obtain accurate results below
the Debye temperature. Additionally, the electronic entropy is intrinsically a
quantum mechanical property.

The following section presents our computational methods. The heart of our
approach rests on approximating the many-body displacement correlation
function as a Gaussian distribution matching the simulated covariance of
atomic displacements. We then apply the method to two test cases, face-
centered cubic Al and body-centered cubic Na. In each case we compare with
standard thermodynamic data. We also show the applicability of a simple
approximation based on experimental Debye-Waller factors (thermal B-factors)
that allow experimental diffraction measurements to obtain approximate
thermodynamic entropies.

Our principal results for these test cases are illustrated in Fig. 1 parts (a)
and (b). Notice first that the Debye-Waller factors yield good qualitative
results, lying within 1$k_{\rm B}$ of the experimental values, but remaining
consistently high. This is because the Debye-Waller factors treat the
individual atomic vibrations independently, and neglect the mutual information
contained in displacement correlation functions that must reduce the
vibrational entropy [1, 13, 14, 15, 6]. Including the covariances of
displacements and electronic entropies (curves labeled classical) improves the
agreement, but with negative deviations at low temperatures due to the
$\log{T}$ divergence of the classical vibrational entropy.

Figure 1:  Entropies of (a) Al, and (b) Na. Black circles show experimental
values from the NIST JANAF Tables [16, 17]. Orange triangles are calculated
from Eq. (16) using B-factors obtained from [18, 19]. Red squares add the
classical vibrational entropy calculated from Eq. (10) to the electronic
entropy calculated from Eq. (15). Parts (c) and (d) show residuals after
subtracting the experimental values. In addition, we show our vibrational
quantum model calculated from Eq. (14) using effective vibrational frequencies
calculated by Eq. (13), and the quasiharmonic prediction with vibrational
frequencies calculated by Phonopy. All calculations are performed at the
experimentally determined volumes for each temperature. All curves, except for
Debye-Waller and single-site, include electronic entropy.

To overcome the deficiency of classical statistical mechanics at low
temperatures, we introduce a quantum version of our method where we interpret
eigenvalues of the covariance matrix as effective vibrational frequencies
renormalized by anharmonic forces. This reveals a relationship between our
method and a different approach based on velocity autocorrelation functions
[20, 21, 22], although the two approaches differ in important respects.
Residuals of several versions of our calculated entropies relative to
experiment are plotted in Fig. 1 (c) and (d), and are further discussed in
Section III.1.

We then apply our method to two examples that are scientifically interesting
and technically challenging. First we examine the high entropy alloy MoNbTaW
[23, 24]. Here the chemically disordered structure makes the conventional
phonon-based approach time consuming. Unfortunately it also increases the
demands on AIMD run times and limits our ability to improve statistics through
symmetrization. Next, we address the temperature-driven HCP to BCC transition
of Ti. Owing to the presence of imaginary frequency modes in the BCC state,
the usual harmonic and quasiharmonic approaches cannot be applied, while our
method succeeds.

##  II Methods

###  II.1 Probability density function

Our approach focuses on the $N$-body probability density function
$\rho_{N}({\bm{\mathcal{U}}},{\bm{\mathcal{P}}})$ of a classical $N$-atom
system in Cartesian phase space. The displacement variable
${\bm{\mathcal{U}}}=({\bf u}_{1},{\bf u}_{2},...,{\bf u}_{N})$, where ${\bf
u}_{i}\equiv{\bf r}_{i}-{\bf R}_{i}$ defines atomic displacement of the
position ${\bf r}_{i}$ of atom $i$ from its mean position (ideal lattice site)
${\bf R}_{i}$, and ${\bm{\mathcal{P}}}=({\bf p}_{1},{\bf p}_{2},...,{\bf
p}_{N})$ incorporates the momenta $\\{{\bf p}_{i}\\}$. Owing to the additivity
of kinetic and potential energy, the phase space probability factors into a
product of density functions $f_{u}({\bm{\mathcal{U}}})$ and
$f_{p}({\bm{\mathcal{P}}})$

$\rho_{N}({\bm{\mathcal{U}}},{\bm{\mathcal{P}}})=h^{3N}f_{u}({\bm{\mathcal{U}}})f_{p}({\bm{\mathcal{P}}}).$
(2)

The factor $h^{3N}$ comes from the constraint that the probability density
integrates to 1,

$\frac{1}{h^{3N}}\iint\cdots\int_{-\infty}^{\infty}\rho_{N}({\bm{\mathcal{U}}},{\bm{\mathcal{P}}})d^{3N}\,{\bm{\mathcal{U}}}\,d^{3N}{\bm{\mathcal{P}}}=1.$
(3)

The entropy according to Gibbs [25] is

$S=-\frac{1}{h^{3N}}\iint\cdots\int_{-\infty}^{\infty}\rho_{N}({\bm{\mathcal{U}}},{\bm{\mathcal{P}}})\ln\rho_{N}({\bm{\mathcal{U}}},{\bm{\mathcal{P}}})\,d^{3N}{\bm{\mathcal{U}}}\,d^{3N}{\bm{\mathcal{P}}}.$
(4)

This is identical to the Shannon [26] information-theoretic entropy, in
suitable units.

According to classical Maxwell-Boltzmann statistics, the momentum distribution
function is Gaussian,

$f_{\bm{\mathcal{P}}}({\bm{\mathcal{P}}})=\frac{\exp(-\frac{1}{2}{\bm{\mathcal{P}}}^{\rm
T}\bm{\Sigma_{p}^{-1}}{\bm{\mathcal{P}}})}{\sqrt{(2\pi)^{3N}\det(\bm{\Sigma}_{p})}},$
(5)

with $\Sigma_{p}$ a diagonal matrix of entries $m_{i}/\beta$ where $m_{i}$ is
the mass of atom $i$ and $\beta=1/k_{\rm B}T$. Formally, we set $\bm{M}={\rm
diag}\left(m_{1},m_{1},m_{1},m_{2}\cdots,m_{N}\right)$, so that
$\Sigma_{p}=\bm{M}/\beta$.

In contrast to the simplicity of the momentum distribution, the density
function $f_{u}({\bm{\mathcal{U}}})$ is difficult to describe precisely,
considering the many-body and anharmonic interactions among atoms. We choose
to approximate it as a Gaussian with suitable covariance. Hence we write

$f({\bm{\mathcal{U}}})=\frac{\exp(-\frac{1}{2}{\bm{\mathcal{U}}}^{\rm
T}\bm{\Sigma_{u}^{-1}}{\bm{\mathcal{U}}})}{\sqrt{(2\pi)^{3N}\det(\bm{\Sigma_{u}})}},$
(6)

where $\bm{\Sigma_{u}}$ is the covariance matrix

$\bm{\Sigma_{u}}=\left(\begin{array}[]{cccc}\bm{\sigma}_{1,1}&\bm{\sigma}_{1,2}&\cdots&\bm{\sigma}_{1,N}\\\
\bm{\sigma}_{2,1}&\bm{\sigma}_{2,2}&\cdots&\bm{\sigma}_{2,N}\\\
\vdots&\vdots&\ddots&\vdots\\\
\bm{\sigma}_{N,1}&\bm{\sigma}_{N,2}&\cdots&\bm{\sigma}_{,NN}\\\
\end{array}\right).$ (7)

The $\bm{\sigma}_{i,j}$ element of $\bm{\Sigma_{u}}$ is the $3\times 3$
covariance matrix of the displacements ${\bf u}_{i}$ and ${\bf u}_{j}$ of the
$i^{\rm th}$ and $j^{\rm th}$ atoms,

$\bm{\sigma}_{i,j}=\left(\begin{array}[]{ccc}{\langle
x_{i}x_{j}\rangle}&{\langle x_{i}y_{j}\rangle}&{\langle x_{i}z_{j}\rangle}\\\
{\langle y_{i}x_{j}\rangle}&{\langle y_{i}y_{j}\rangle}&{\langle
y_{i}z_{j}\rangle}\\\ {\langle z_{i}x_{j}\rangle}&{\langle
z_{i}y_{j}\rangle}&{\langle z_{i}z_{j}\rangle}\end{array}\right),$ (8)

with $x$, $y$, and $z$ the Cartesian coordinates of the displacement ${\bf
u}$. Diagonal elements of the covariance matrix yield the variances, e.g. for
our cubic lattices $\bm{\sigma}_{i,i}={\langle x^{2}\rangle}\bm{1}$. Due to
the Gaussian approximation, the many-body density $f_{u}({\bm{\mathcal{U}}})$
factors into a product of two-body correlations. Note that these two-body
terms include anharmonic effects through the values of their covariances.

Within these approximations, the entropy $S$ of N atoms becomes

$S=\frac{1}{2}\ln(\det{(\bm{\Sigma_{u}})})+\frac{3}{2}\sum_{i=1}^{N}\ln(m_{i}/\beta\hbar^{2})+3N.$
(9)

If all masses are equal, $S$ simplifies to

$S=\frac{1}{2}\ln\left((2\pi e\Lambda)^{3N}\det{(\bm{\Sigma_{u}})}\right)$
(10)

where $\Lambda=\sqrt{2\pi\hbar^{2}/mk_{B}T}$ is the thermal de Broglie
wavelength for mass $m$ at temperature $T$. Subject to the Gaussian
approximation, our method resembles the approach of Morris and Ho [1], who
applied it to a one-dimensional model system. However, the formalism of Eq.
(4) applies generally, and we will examine corrections to the Gaussian
approximation in Section III.3.2.

Fig. 2 illustrates the covariance matrix $\bm{\Sigma_{u}}$ for FCC Al at
T=900K. Repeating patterns reflect the symmetries of the FCC structure.
Translational symmetry requires that the covariance submatrix
$\bm{\sigma}_{i,j}$ depends only on the relative position
$\bm{R}_{i,j}=\bm{R}_{j}-\bm{R}_{i}=h\bm{a}+k\bm{b}+l\bm{c}$, of the $i^{th}$
and $j^{th}$ atoms. Consequently, covariance matrices $\bm{\sigma}_{i,j}$
sharing the same Miller indices $hkl$ share the same value,
$\bm{\sigma}_{hkl}$. All $3\times 3$ matrices along the diagonal are
equivalent and share the form $\bm{\sigma}_{000}$ shown in part (d), whose off
diagonal elements vanish due to mirror symmetries. Three-fold rotational
symmetry can be seen in the covariance matrices $\bm{\sigma}_{011}$,
$\bm{\sigma}_{101}$, $\bm{\sigma}_{110}$ (parts (e)-(g)) whose non-zero off-
diagonal elements are $yz$, $xz$, and $xy$ components.

Figure 2:  (a) Complete AIMD simulated covariance matrix of FCC Al at 900K in
a $4\times 4\times 4$ supercell of 256 atoms. (b) Submatrix of a 4-atom
tetrahedron. (c) Unit cell of FCC Al illustrating tetrahedron of four nearest
neighbors. (d)-(g) Single site variance matrix $\bm{\sigma}_{000}$ and three
nearest-neighbor covariance matrices $\bm{\sigma}_{011}$, $\bm{\sigma}_{101}$
and $\bm{\sigma}_{110}$. Red color indicates positive covariance while blue
color indicates negative covariance. Color bar indicates ${\rm
sgn}(\sigma)\ln{(|\sigma/\sigma_{\rm min}|)}$.

###  II.2 Relation to force constant matrix

The probability density $\rho(x)$ of a classical oscillator in the harmonic
potential $U=\frac{1}{2}m\omega^{2}x^{2}$ is

$\rho(x)=\sqrt{\frac{\beta m\omega^{2}}{2\pi}}e^{-\frac{1}{2}\beta
m\omega^{2}x^{2}},$ (11)

and the variance of its displacement is
$\sigma^{2}=\left<x^{2}\right>=1/(\beta m\omega^{2})$. The force constant
$C=U^{\prime\prime}=m\omega^{2}$ is related to the variance by
$C=1/\beta\sigma^{2}$. For an $N$-particle system, the force constant matrix
$C$ is defined in term of the second derivative of the potential $U$,

$C_{i\mu,j\nu}=\frac{\partial^{2}U}{\partial u_{i\mu}\partial u_{j\nu}}.$ (12)

where $u_{i\mu},u_{j\nu}$ are elements of displacement ${\bm{\mathcal{U}}}$ in
which $i,j$ denote atoms and $\mu,\nu$ denote $x,y,z$ Cartesian coordinates.
The mass-weighted covariance matrix,
$\tilde{\Sigma}_{i\mu,j\nu}=\sqrt{m_{i}m_{j}}\,\Sigma_{i\mu,j\nu}$, relates to
the mass-reduced force constant matrix
$\tilde{C}_{i\mu,j\nu}=C_{i\mu,j\nu}/\sqrt{m_{i}m_{j}}$, by

$\tilde{C}=\frac{1}{\beta}\tilde{\Sigma}^{-1},$ (13)

hence measurement of the covariance matrix yields the complete set of force
constants. The matrices $\tilde{C}$ and $\tilde{\Sigma}$ are singular because
of center of mass translation invariance. To invert the singular matrix, we
represent $\tilde{\Sigma}=\sum_{k\mu}\lambda_{k\mu}|k\mu\rangle\langle k\mu|$
where $\\{(\lambda_{k\mu}\equiv\beta\omega_{k\mu}^{2})^{-1},|k\mu\rangle\\}$
is the set of eigenvalues and eigenvectors of $\tilde{\Sigma}$. Then, noting
that $\tilde{C}$ and $\tilde{\Sigma}$ share common eigenvectors, we invert the
nonvanishing eigenvalues to obtain
$\tilde{C}=\sum_{k\mu}\omega_{k\mu}^{2}|k\mu\rangle\langle k\mu|$.

For a harmonic potential $U$, the relationship Eq. (13) is exact; for an
anharmonic system we may take Eq. (13) as defining temperature-dependent
effective force constants.

###  II.3 Quantum harmonic entropy

The entropies predicted by our classical theory agree quite well with the
experimental values at high temperatures, but they fall below experiment at
temperatures below the Debye temperatures $\Theta_{D}$, as seen in Fig. 1. The
negative deviation is a consequence of the negative divergence of
$\log{(u^{2}/\Lambda^{2})}\sim 2\log{T}$ as $T\rightarrow 0$. Experimentally
$S\rightarrow 0$ for all materials, by the third law of thermodynamics,
because quantum mechanics inhibits the excitation of vibrational modes with
frequencies greater than $k_{\rm B}T/\hbar$.

To overcome the singularity of classical entropy, we adopt entropy of the
quantum harmonic oscillator, using harmonic frequencies $\omega_{k\mu}$
obtained from eigenvalues of our covariance matrix as discussed in Sec. II.2.
Summing over the nonzero vibrational frequencies, the entropy with quantum
corrections is

$S=k_{B}\sum_{k\mu}\left[-\ln(1-e^{-\beta\hbar\omega_{k\mu}})+\frac{\beta\hbar\omega_{k\mu}}{e^{\beta\hbar\omega_{k\mu}}-1}\right].$
(14)

This yields better agreement when temperature is below the Debye temperature
as shown in Fig. 1. In particular, the limit $S\rightarrow 0$ as $T\rightarrow
0$ is obeyed.

This quantum model is harmonic in the sense that it is exact for quadratic
potentials $U$, but it incorporates anharmonicity through the effective
vibrational frequencies which were derived from the simulated covariance
matrix. Errors due to applying the quantum harmonic model should be small at
low temperatures, where motion generically becomes harmonic. Some prior
studies employ time-dependent velocity correlation functions, then Fourier
transform over time to obtain frequencies [20, 21, 22]. The systematics of
that approach differ markedly from ours, as in principle we do not require
time evolution at all; we only simulate trajectories for the sake of enlarging
our configurational ensemble.

The model Hamiltonian can be constructed in the actual harmonic limit of small
oscillations by evaluating the force constants within density functional
perturbation theory. This mode substantially underestimates the high
temperature entropy as it neglects thermal expansion. The quasiharmonic
approximation can be used to predict thermal expansion, resulting in improved
agreement, or better yet we can evaluate the force constants at the
experimental lattice parameters. As seen in Fig. 1 the quasiharmonic
approximation utilizing experimental lattice constants agree with experiment
about as well as our new method.

###  II.4 Ab-initio methods

Ab initio molecular dynamics (AIMD) simulations are performed for FCC Al in
supercells of size $4\times 4\times 4$ (256 atoms) and $6\times 6\times 6$
(864 atoms), and for BCC Na in a $6\times 6\times 6$ supercell (432 atoms). We
use the Vienna Ab initio Simulation Package (VASP [27]) using augmented plane
wave potentials [28] with the Perdew-Burke-Ernzerhof (PBE [29]) generalized
gradient exchange correlation functional. We use a single electronic $k$-point
and default plane wave energy cutoffs. When possible we use experimental
lattice constants at the appropriate temperatures. The molecular dynamics
simulations use Nosé thermostats with the default Nosé mass parameters. Our
time steps are 2fs, and our runs extend to 40ps for Al ($4\times 4\times 4$)
and 8ps for ($6\times 6\times 6$), and 7ps or greater for Na.

After allowing the simulated systems to approach equilibrium, the variances
and covariances are calculated from a continuing simulation by averaging ${\bf
u}_{i}{\bf u}_{j}$ over many samples. We also average over $\Omega$
reflection, rotation and translation symmetry operations ${\bf T}_{k}$ such
that $\boldmath{\sigma}_{i,j}=\frac{1}{\Omega}\sum_{k}{\bf T}_{k}{\bf
u}_{i}{\bf u}_{j}$ becomes symmetry invariant. In principle all the
information needed to evaluate the entropy is contained in just a single
representative structure of sufficient size, but the time averaging helps to
reduce statistical error.

We perform phonon calculations as implemented in phonopy [30] to obtain force
constants and vibrational frequencies, and then calculate vibrational entropy
as discussed in Section II.3. Rather than calculating the thermal expansion
ab-initio, as in the traditional quasiharmonic approximation [31], we simply
evaluate the force constants at the experimentally known temperature-dependent
lattice constant $a(T)$.

Electronic entropy is evaluated as

$\displaystyle S=$ $\displaystyle-k_{B}\int{\rm
d}ED(E)[f_{T,\mu}(E)\ln{f_{T,\mu}(E)}$ (15)
$\displaystyle+(1-f_{T,\mu}(E))\ln{(1-f_{T,\mu}(E))}]$

with $D(E)$ the electronic density of states calculated from a structure with
lattice constant $a(T)$, and $f_{T,\mu}$ the Fermi-Dirac occupation function.
The chemical potential $\mu$ is obtained as a function of $T$ using the
program Felect [32].

##  III Applications

###  III.1 Test cases: FCC Al and BCC Na

Our method successfully predicts vibrational entropy for Al and Na, as shown
in Fig. 1 parts (a) and (b). Fig. 1 parts (c) and (d) compare the residual
errors of various approximations by subtracting off the experimental
entropies. Curves labeled “Debye-Waller” and “single-site” neglect
correlations among the displacements of different atoms. In this case the
entropy reduces to

$S=\frac{3}{2}\ln{\left[2\pi e(\sigma_{x}/\Lambda)^{2}\right]}$ (16)

where $\sigma_{x}^{2}=\langle u_{x}^{2}\rangle$ is the mean square
displacement. This quantity is related to the Debye-Waller factor [33] that
diminishes the diffraction intensity of a peak or wavevector $\bm{q}$ by the
factor $\exp{(-q^{2}\langle u^{2}\rangle/3)}$. The displacements are sometimes
given in terms of $B=8\pi^{2}\langle u^{2}\rangle/3$. In Fig. 1 we compare the
experimental entropies of Al and Na with the prediction of Eq. (16) using
experimental values of the B-factor. Given the seeming disparity between
crystallographic and thermodynamic methods, the agreement is quite striking.

Note that the Debye-Waller and single-site entropies exceed the experimental
values. The displacement of a single atom applies forces that displace nearby
atoms, reducing the total amount of information needed to specify a given
configuration ${\bm{\mathcal{U}}}$. A similar effect is found in the entropy
of liquids, where the information content of pair correlation functions
reduces the entropy below the value for an ideal gas at the same overall
density [13, 14, 15]. The improvement upon including the full covariance
matrix is evident in the curve labeled “classical”, but as discussed
previously it suffers an unavoidable $\ln{T}$ negative divergence. This
divergence is aleviated at low temperatures through quantum model (Sec. II.3).
The quasiharmonc model is also quite accurate.

To better understand how the covariance matrix and entropy are influenced by
the range of correlations, and by our finite MD simulation cells, we study the
convergence of covariance matrix elements and corresponding entropy of Al,
including only matrix elements $\bm{\sigma}_{hkl}$ of pairs separated by
$R_{hkl}=|h{\bf a}+k{\bf b}+l{\bf c}|$. Fig. 3 (a) and (b) show that the
absolute value of $\det{(\bm{\sigma}_{hkl})}$ drops rapidly with increasing
the bond length, suggesting our simulation cell size is sufficient to capture
the dominant collective motions of the solid, although some indication of cell
size dependence can be seen in the excess correlation around $[hkl]=004$ at
$T=300$K. Similar decay of correlations was observed in other simulations [6,
34]. Comparing $T=900$K with $T=300$K, we see similar variation with
$R_{hkl}$, while the values at high temperature are nearly two orders of
magnitude larger. Comparing convergence of the $4\times 4\times 4$ (256 atom)
cell with the $6\times 6\times 6$ (864 atom) cell in Fig. 3 (c) and (d)
suggests the $4\times 4\times 4$ is fully adequate for entropy calculation at
high temperatures, but just barely sufficient at low temperature.

Figure 3: Left: Comparison of covariance matrix elements $\bm{\sigma}_{hkl}$
as bond length increases for Al at (a) $T=300$K and (b) $T=900$K. Right:
Convergence of entropy $S$ after including all covariance matrix elements
$\bm{\sigma}_{hkl}$ of pairs within $R<R_{hkl}$ at (c) T=300K and (d)
$T=900$K.

###  III.2 High entropy alloy: vibrational entropies of MoNbTaW

Although high entropy alloys (HEAs) acquire their name from the entropy of
chemical substitution, their vibrational entropy may exceed their
substitutional entropy by a considerable margin. Substitutional entropy is
relevant for stability mainly because the vibrational entropy of the mixture
lies close to the average vibrational entropy of the elements [35]. Here, we
investigate the applicability of our covariance method to calculate the
vibrational entropy of MoNbTaW [23]. Since chemical substitution is prevalent
in HEAs, we have to choose what specific arrangement of atoms to take. We will
take as representative structures the final configurations from hybrid MC/MD
simulations [24], which reflect the temperature variation of chemical order.

Figure 4: (a) Comparison plot of the average of experimental elemental
entropies [36, 37, 38, 39], $S_{\rm avg}$, with the entropy $S^{v}_{\rm
qha}+S^{e}$ calculated using the quasiharmonic approximation, the entropy
$S^{v}_{\rm cm}+S^{e}$ calculated from the covariance matrix, and the entropy
$S^{v}_{\rm pacm}+S^{e}$ calculated from the pair-averaged covariance matrix.
(b) Convergence of vibrational entropy $S^{v}_{\rm cm}$ as simulation time
increases compared with vibrational entropies $S^{v}_{\rm qha}$ and
$S^{v}_{\rm pacm}$.

We calculate the vibrational entropy $S^{v}_{\rm cm}$ of a specific chemical
configuration at each temperature using the covariance matrix $\Sigma_{\bf u}$
obtained from from MD simulations. Fig. 4 (a) plots entropies $S^{v}_{\rm
cm}+S^{e}$ of theses structures. We compare our prediction with the average
experimental entropies of pure elements, $S_{\rm avg}$, and with the
quasiharmonic vibrational entropies $S_{\rm qha}^{v}+S^{e}$ of a cF16
(Heusler) MoNbTaW structure at the same lattice parameters as our MD
simulations. These temperature-dependent lattice parameters were determined by
varying the volume until the simulated total pressures vanish on average. It
is seen from Fig. 4 (a) that both quasiharmonic and covariance matrix
entropies are close to, but slightly smaller than, the averaged entropy
$S_{\rm avg}$ of pure elements, consistent with prior calculations [40].

The vibrational entropy derived from the covariance matrix converges slowly
because these chemically disordered structures lack symmetry and we cannot
employ symmetry averaging as discussed in Section II.4. As a result, the
covariance matrix has poor statistics and is hard to converge as illustrated
in Fig. 4 (b). Unfortunately we lack an extrapolation formula for entropy vs.
simulation time. At long times these entropies converge towards entropies
calculated from the quasiharmonic approximation $S_{\rm qha}$.

In an effort to alleviate the poor statistics, we introduce a “pair averaged
covariance matrix”, $\bar{\Sigma}_{\bf u}$, that maintains the chemical
identities at each site while averaging of their chemical environments. The
$(i,j)$ element of the full covariance matrix $\Sigma_{\bf u}$ is the $3\times
3$ matrix $\bm{\sigma}_{i,j}^{\alpha\beta}$, where the superscripts remind us
that the chemical species at site $i$ is $c(i)=\alpha$ and the chemical
species at site $j$ is $c(j)=\beta$. Let the $\mathcal{P}_{i,j}^{\alpha\beta}$
be the set of all pairs $(i^{\prime},j^{\prime})$ such that ${\bf
R}_{i^{\prime},j^{\prime}}={\bf R}_{i,j}$ and $c(i^{\prime})=\alpha$ and
$c(j^{\prime})=\beta$. We define the $(i,j)$ element of $\bar{\Sigma}_{\bf u}$
as

$\bar{\Sigma}_{\bf
u}(i,j)=\bar{\bm{\sigma}}_{i,j}^{\alpha\beta}=\frac{1}{N_{i,j}^{\alpha\beta}}\sum_{i^{\prime},j^{\prime}}\bm{\sigma}_{i^{\prime},j^{\prime}}^{\alpha\beta}$
(17)

where the sum runs over the set $\mathcal{P}_{i,j}^{\alpha\beta}$ containing
$N_{i,j}^{\alpha\beta}$ elements. The entropy computed from $\bar{\Sigma}_{\bf
u}$ is expected to provide a close upper bound on $S_{\rm cm}^{v}$.

###  III.3 BCC to HCP phase transition in titanium

Certain elements and compounds are so strongly anharmonic that the entropy
simply cannot be calculated within the harmonic or quasiharmonic
approximation. Elements in columns 3 and 4 of the Periodic Table undergo
diffusionless (Martensitic) phase transformations from BCC ($\beta$-phase)
stable at high temperature to HCP ($\alpha$-phase) stable at low temperature.
Harmonic analysis predicts their BCC states to be mechanically unstable at low
temperature because they exhibit imaginary vibrational frequency modes.
Eigenvectors of these modes describe the transformation pathway [41, 42]. The
instability prevents application of conventional harmonic or quasiharmonic
calculations of the entropy. Our calculation method circumvents this
difficulty because it does not require the calculation of vibrational
frequencies.

These structural phase transitions are of practical importance, motivating
considerable efforts to predict transition temperatures and understand their
mechanisms [43, 44, 45, 46, 47, 48]. Proposed methods include phase space
partitioning [43, 44, 45], effective force constant averaging [48], and an
“augmented lattice” model [44]. Predicted transition temperatures range from
1095K to 1114K, in general agreement with in agreement with the experimental
transition temperature $T_{c}=1166K$ [49]. We apply our covariance matrix
method to calculate vibrational entropy and predict the transition temperature
$T_{c}=1060$K.

Figure 5:  (a) Lattice constant $a$ of BCC Ti. [50, 51, 52, 53, 54] (b) and
(c) Lattice constants $a$ and $c$ of HCP Ti.[50, 55, 56, 53, 57] Symbols in
black are selected as fitting database. (d) Comparison of density $\rho$
between BCC and HCP Ti.

We perform AIMD simulations for both BCC and HCP Ti at lattice constants that
are fitted to experimental measurements with quadratic functions as shown in
Fig. 5. Considering the scattering of experimental measurements of lattice
constants, we choose to fit Ref. [50, 53] for lattice parameters of BCC Ti and
Ref.[50, 55, 56] for HCP Ti. To minimize size effect, we prepare simulation
cells with the same atomic number—an orthorhombic 256-atom 4x4x4 supercell
based on a 4-atom unit cell $(\bm{a}=a,0,0;\bm{b}=0,a,-a;\bm{c}=0,a,a)$ for
BCC Ti, and an orthorhombic 256-atom 4x4x4 supercell based on a 4-atom unit
cell $(\bm{a}=a,0,0;\bm{b}=0,\sqrt{3}a,0;\bm{c}=0,0,c)$ for HCP Ti.

Figure 6:  Comparison of calculated entropy and experimental entropy of
$\alpha$-phase Ti and $\beta$-phase Ti. (b) Comparison of calculated and
experimental enthalpy. (c) Comparison of experimental and calculated free
energies. Yellow and green backgrounds shade regions of stability of the
$\alpha$ and $\beta$ phases, respectively, as determined by experiment.
Experimental data is plotted with solid lines in regions of stability, and
dashed lines in regions of instability. Experimental data comes from NIST-
JANAF Thermochemical Tables.[49, 58, 59]. Calculated entropies
$S=S^{v}+S^{e}$.

A comparison of calculated total entropy $S^{v}_{cm}+S^{e}$ and experimental
entropy is illustrated in Fig. 6 (a). Electronic entropies $S^{e}$ are
calculated from Eq. 15 with electronic density states obtained at the given
volume for each temperature. As shown in Fig. 7, BCC Ti has a substantially
higher electronic entropy than HCP Ti due to the pseudogap at the Fermi energy
of the HCP density of electronic states. Formation of the pseudogap drives the
Burger’s distortion from BCC to HCP [42]. Entropy of HCP Ti from our work
compares well to the experimental entropy except one value at $T=1400$K which
falls in the region where HCP is thermodynamically unstable. The entropy of
BCC Ti, however, is overestimated by an amount of $0.5k_{\rm B}$ to $1.0k_{\rm
B}$ at all temperatures.

Enthalpies are obtained by averaging energies over our MD simulations. To
place enthalpies on the experimental scale, we shifted all of our calculated
enthalpy values so that our enthalpy of $\alpha$ matched the experimental
value at $T=800$K. For both phases our simulation matches well with
measurement at temperatures below the $\alpha\to\beta$ temperature while at
higher temperatures it falls below the experimental enthalpy. Finally, we
compute the Gibbs free energy $G=H-TS$ based on our calculated entropy and
enthalpy of HCP Ti at T=$800$K, $1000$K, $1200$K and BCC Ti at T=$1200$K,
$1400$K, $1600$K and predict $\alpha$-$\beta$ phase transition temperature
$T_{c}^{\rm pred}=1060K$ (see Fig. 6).

Figure 7:  (a) Calculated electronic entropies $S^{e}$ of BCC and HCP Ti. (b)
Electronic densities of states of BCC and HCP Ti evaluated at their volumes at
$T=1200$K. Fermi smearing of width $\sigma=k_{\rm B}T=0.103$eV has been
applied.

To understand the overestimate of BCC entropy, which leads to a low estimate
of $T_{c}$, we compare calculated phonon spectra and vibrational density of
states derived from our force constant matrix (see Fig. 8) with results from
Ref. [54]. Note that our effective vibrational frequencies fall systematically
below the experiment, explaining the overestimate of entropy. We tested to see
if this could be due to errors in lattice constant, but the impact of volume
changes was not sufficient to explain our disagreement. Presumably the fault
lies in some aspect of our simulation method. Below we investigate possible
explanations in finite size effects, or anharmonicity, but these also turn out
to be too small to explain the discrepancy.

Figure 8:  (a) Calculated vibrational dispersion spectral of BCC Ti at T=1200K
in a cubic cell with 250 atoms. (b) Calculated vibrational density of states.
Red and blue dots in (a), and red curve in (b) come from Petry [54].

####  III.3.1 Finite site effect

To evaluate the impact of simulated cell size on the entropy of BCC Ti, we
perform entropy calculation for three sizes: 54-atoms, 128-atoms, and
250-atoms. Fig. 9 shows a linear relation between entropy $S^{v}_{cm}$ and
inverse size $1/N$. With larger cells, entropy increases, and so does the
disagreement with experiment. This finite size effect for BCC Ti resembles the
finite size effect in high-pressure high-temperature BCC Fe [60], so we
believe the effect is real.

Figure 9:  Calculated vibrational entropy $S^{v}_{cm}$ versus the inverse of
the simulated number of atoms, $1/N$, for BCC Ti at $T=1200$K. Experiment [54]
is at $T=1208$K.

####  III.3.2 Anharmonicity

We investigate the effect of anharmonic corrections to the single site
probability density

$p^{h}({\bf
u})=\frac{e^{-\frac{x^{2}+y^{2}+z^{2}}{2\sigma^{2}}}}{\sqrt{(2\pi\sigma^{2})^{3}}}=\frac{e^{-\frac{R^{2}}{2\sigma^{2}}}}{\sqrt{(2\pi\sigma^{2})^{3}}}.$
(18)

At the lowest order of anharmonicity, the probability density $p^{a}({\bf u})$
includes the isotropic term

$I({\bf u})=(x^{2}+y^{2}+z^{2})^{2},$ (19)

and the anisotropic term

$A({\bf
u})={{x}^{4}}+{{y}^{4}}+{{z}^{4}}-3\left({{x}^{2}}{{y}^{2}}+{{x}^{2}}{{z}^{2}}+{{y}^{2}}{{z}^{2}}\right),$
(20)

which are invariant under cubic symmetry operations. The anharmonic
probability density is hence approximated by

$p^{a}({\bf
u})=\frac{1}{Z}\exp\left(-\frac{R^{2}}{2\sigma^{2}}\right)\exp\left(-\frac{a}{4\sigma^{4}}I({\bf
u})-\frac{b}{4\sigma^{4}}A({\bf u})\right),$ (21)

where $Z$ is the normalization factor

$Z=\int_{V}\,d{\bf u}~{}p^{a}({\bf u}),$ (22)

and the integration volume $V$ is the Wigner-Seitz cell of an atom.

In practice we cut off the integration at the cube
$V=[-8\sigma,+8\sigma]^{3}$, as justified by the rapid vanishing of
$p^{a}({\bf u})$. We calculate averages ${\langle R^{2}\rangle}$, ${\langle
R^{4}\rangle}$, and ${\langle A\rangle}$ during our simulation, then we fit
values of $\sigma$, $a$, and $b$ by solving the simultaneous nonlinear
equations

$\displaystyle{\langle R^{2}\rangle}$ $\displaystyle=$
$\displaystyle\int_{V}\,d{\bf u}~{}R^{2}~{}p^{a}({\bf u})$ (23)
$\displaystyle{\langle R^{4}\rangle}$ $\displaystyle=$
$\displaystyle\int_{V}\,d{\bf u}~{}R^{4}~{}p^{a}({\bf u})$ (24)
$\displaystyle{\langle A\rangle}$ $\displaystyle=$
$\displaystyle\int_{V}\,d{\bf u}~{}A({\bf u})~{}p^{a}({\bf u})$ (25)

where the probability $p^{a}({\bf u})$ is given by Eq. (21. Finally, the
positional part of the anisotropic entropy, $S^{a}$, is calculated from

$S^{a}=-k_{\rm B}\int_{V}\,d{\bf u}~{}p^{a}(\bm{u})\ln p^{a}(\bm{u}).$ (26)

|  $\left<R^{2}\right>$ [Å] |  $\left<R^{4}\right>$ [Å4] |  $\left<I\right>$ [Å4] |  $\sigma^{2}$ [Å2] | $a$ | $b$ |  $S^{a}-S^{h}$ [$k_{\rm B}$]
---|---|---|---|---|---|---|---
FCC Al | 0.05933 | 0.005945 | 0.00003873 | 0.01918 | -0.005718 | -0.001890 | -0.00017975
BCC Ti | 0.19061 | 0.06040 | -0.002597 | 0.06433 | 0.003210 | 0.01596 | -0.001227
Table 1:  Statistical average of $\left<R\right>$, $\left<R^{2}\right>$ from
MD simulations. Correction to harmonic entropy and parameters $\sigma^{2}$,
$a$, $b$, and $S^{a}-S^{h}$ of FCC Al at $T=500$K and BCC Ti at $T=1200$K.

Table. 1 compares the influence of anharmonicity in FCC Al with BCC Ti, and
presents numerical values of the averages in Eq. (23) and the solutions for
$\sigma^{2},a,b$ and entropy $S$. Anharmonicity tends to reduce the entropy
for both FCC Al and BCC Ti, yet the reduction is insufficient to explain our
entropy excess in BCC Ti. Differences in the signs of the $a$ and $b$
parameters between FCC Al and BCC Ti imply opposite deviations of our harmonic
model from the simulated distribution. In FCC Al, the simulated distribution
is more narrow, with a higher probability at origin than in our harmonic
model; in BCC Ti, the simulated distribution is broader and lower at the
origin. Our anharmonic model captures these deviations, as shown in the
marginal distributions $p(x)$ in Fig. 10.

Figure 10:  Marginal probability distributions $p(x)$ of (a) FCC Al at T=500K
and (b) BCC Ti at T=1200K. Crosses are histograms of the simulated data, red
lines are fits to the harmonic model $p^{h}$ (Eq. (18)), and blues lines are
fits to the anharmonic model $p^{a}$ (Eq. (21)).

To examine the anisotropies, we plot the marginal distributions $p^{a}(x,y)$
in Fig. 11. In FCC Al, atomic displacements are reduced in the near-neighbor
directions [110] and correspondingly enhanced in the [100] directions. In BCC
the displacements are reduced in the nearest-neighbor directions [111] (not
shown). Four such bonds project onto the [100] directions, while only two
project onto [110], explaining the observed pattern. Overall, FCC Al is more
isotropic than BCC Ti and hence has a smaller angular entropy correction.

Figure 11: Difference of harmonic and anharmonic marginal probability
distributions, $p^{a}(x,y)-p^{h}(x,y)$ for FCC Al (Left) and BCC Ti (Right).

##  IV Conclusion

We apply the information-theoretic entropy formula Eq. (4) to evaluate the
vibrational entropies of solids from the variance and covariance of atomic
displacements. This approach generalizes prior work on the information-based
entropy of liquids [13, 14, 15]. In the case of liquids, the single atom
entropy (ideal gas term) overestimates the entropy and must be corrected by
removing the mutual information of the pair correlation functions. In the case
of solids, the variance of individual atomic displacements can be measured
through diffraction experiments that yield the Debye-Waller B-factor. Thus we
find a crystallographic approach to estimate the thermodynamic entropy.
However, as in the case of liquids, the one-body approximation overestimates
the entropy by the information content of correlation functions, and we can
improve the entropy estimate by including the covariance of atom pairs. This
might be possible to achieve through diffraction experiments that measure the
second-order thermal diffuse scattering [33]. It is easy to achieve through
AIMD simulations of the atomic displacement covariance matrix, as we
demonstrate in this paper for elemental Al and Na.

The method applies generally to solids, but the particular implementation
given here relies on the accuracy of a Gaussian approximation to the
distribution function. Hence it is most likely to work when the atomic
displacements are small, and it is likely to fail in molecular solids where
coherent bond rotations are present. Although we mainly demonstrated the
method for elemental solids, it also holds in principle for complex
crystalline and noncrystalline solids. We give an example of such an
application for the MoNbTaW high entropy alloy.

The quasiharmonic method may be equally accurate and more efficient than our
AIMD method when anharmonicity mainly enters through thermal expansion, but a
simulation-based approach in principle includes additional anharmonic
contributions. Doing so may require correlations beyond those captured by the
Gaussian approximation (see Sec. III.3.2). Our simulation-based approach seems
most useful when the simulation has already been completed for other purposes.
Then the entropy comes essentially for free on top of whatever other
information was sought.

In some cases the quasiharmonic method cannot be applied due to the presence
of imaginary frequency vibrational modes. The high temperature BCC phases of
columns 3 and 4 of the periodic table exhibit such modes; they achieve
mechanical stability only through their entropies. For elemental Ti, our AIMD
method is capable of estimating the vibrational entropy, although the modes
seem slightly softer than observed in experiment. We also point out the
unexpected strong contribution to stability from the electronic entropy.

Our simulation approach based on the probability distribution is more flexible
than the velocity-velocity correlation method because it does not rely on an
underlying harmonic model, at least in the high temperature limit. Further,
because it does not depend upon the dynamics, it can be used in conjunction
with Monte Carlo simulation in addition to molecular dynamics. It requires
only a single representative configuration provided the cell is sufficiently
large, rather than demanding a long continuous trajectory.

###### Acknowledgements.

This work was supported by the Department of Energy under grant DE-SC0014506.
We thank Don Nicholson for useful discussions and for alerting us to Reference
[6], and we thank Michael C. Gao for prior collaboration on the entropy of
liquid metals. We benefited from computer time at the National Energy Research
Scientific Computing Center (NERSC) using award ERCAP0015745.

## References

  * Morris and Ho [1995] J. R. Morris and K. M. Ho, Calculating accurate free energies of solids directly from simulations, Phys. Rev. Lett. 74, 940 (1995).
  * Meirovitch  _et al._ [2009] H. Meirovitch, S. Cheluvaraja, and R. P. White, Methods for calculating the entropy and free energy and their application to problems involving protein flexibility and ligand binding, Curr. Prot. Pept. Sci. 10, 229 (2009).
  * Widom [2016] M. Widom, Entropy and diffuse scattering: comparison of NbTiVZr and CrMoNbV, Met. Mat. Trans. A (2016).
  * Sutton and Levchenko [2020] C. Sutton and S. V. Levchenko, First-principles atomistic thermodynamics and configurational entropy, Front. Chem. 8, 757 (2020).
  * Nir  _et al._ [2020] A. Nir, E. Sela, R. Beck, and Y. Bar-Sinai, Machine-learning iterative calculation of entropy for physical systems, Proc. Nat. Acad. Sci. 117, 30234 (2020).
  * Nicholson  _et al._ [2021] D. M. Nicholson, C. Y. Gao, M. T. McDonnell, C. C. Sluss, and D. J. Keffer, Entropy pair functional theory: Direct entropy evaluation spanning phase transitions, Entropy 23, 234 (2021).
  * Frenkel and Smit [2002] D. Frenkel and B. Smit,  _Understanding Molecular Simulation: From Algorithms to Applications_ , 2nd ed. (Academic Press, 2002).
  * Kastner and Thiel [2005] J. Kastner and W. Thiel, Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: ”umbrella integration”, J. Chem. Phys. 123, 144104 (2005).
  * Grabowski  _et al._ [2019] B. Grabowski, Y. Ikeda, P. Srinivasan, F. Körmann, C. Freysoldt, A. I. Duff, A. Shapeev, and J. Neugebauer, Ab initio vibrational free energies including anharmonicity for multicomponent alloys, npj Comp. Mat. 5, 80 (2019).
  * Fultz [2010] B. Fultz, Vibrational thermodynamics of materials, Prog. Mat. Sci. 55, 247 (2010).
  * Widom [2018] M. Widom, Modeling the structure and thermodynamics of high entropy alloys, J. Mat. Resch. 33, 2881 (2018).
  * Ma  _et al._ [2015] D. Ma, B. Grabowski, F. Körmann, J. Neugebauer, and D. Raabe, Ab initio thermodynamics of the cocrfemnni high entropy alloy: Importance of entropy contributions beyond the configurational one, Acta Mat. (2015).
  * Gao and Widom [2018] M. C. Gao and M. Widom, Information entropy of liquid metals, J. Phys. Chem. B 122, 3550 (2018).
  * Widom and Gao [2019] M. Widom and M. C. Gao, First principles calculation of the entropy of liquid aluminum, Entropy 21, 131 (2019).
  * Huang  _et al._ [2021] Y. Huang, M. C. Gao, and M. Widom, Ab-initio free energies of liquid metal alloys: application to the phase diagrams of Li-Na and K-Na, Phys. Rev. Mater. 6, 013802 (2021).
  * Allison [1996a] T. Allison, JANAF Thermochemical Tables, NIST Standard Reference Database 13 (1996a), https://janaf.nist.gov/tables/Al-001.html.
  * Allison [1996b] T. Allison, JANAF Thermochemical Tables, NIST Standard Reference Database 13 (1996b), https://janaf.nist.gov/tables/Na-001.html.
  * Nakashima [2019] P. N. Nakashima, The crystallography of aluminum and its alloys, Encyclopedia of Aluminum and Its Alloys 10.1201/9781351045636-140000245 (2019).
  * Foust [1972] O. J. Foust,  _Sodium-Nak engineering handbook_ (Gordon and Breach, New York, 1972).
  * Rahman [1964] A. Rahman, Correlations in the motion of atoms in liquid argon, Phys. Rev. (1964).
  * Dickey and Paskin [1969] J. M. Dickey and A. Paskin, Computer simulation of the lattice dynamics of solids, Phys. Rev. 188, 1407 (1969).
  * Lin  _et al._ [2003] S.-T. Lin, M. Blanco, and W. A. Goddard, The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids, J. Chem. Phys. (2003).
  * Senkov  _et al._ [2010] O. N. Senkov, G. B. Wilks, D. B. Miracle, C. P. Chuang, and P. K. Liaw, Refractory high-entropy alloys, Intermetallics 18, 1758 (2010).
  * Widom  _et al._ [2014] M. Widom, W. Huhn, S. Maiti, and W. Steurer, Hybrid Monte Carlo/molecular dynamics simulation of a refractory metal high entropy alloy, Met. Mat. Trans. A 45, 196 (2014).
  * Gibbs [1902] J. W. Gibbs,  _Elementary Principles in Statistical Mechanics_ (Yale, 1902).
  * Shannon [1948] C. E. Shannon, A Mathematical Theory of Communication, Bell System Technical Journal 27, 379 (1948).
  * Kresse and Furthmuller [1996] G. Kresse and J. Furthmuller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, First principles phonon calculations in materials science, Scr. Mater. 108, 1 (2015).
  * Togo  _et al._ [2010] A. Togo, L. Chaput, I. Tanaka, and G. Hug, First-principles phonon calculations of thermal expansion in Ti3SiC2, Ti3AlC2, and Ti3GeC2, Phys. Rev. B 81, 174301 (2010).
  * [32] Electronic free energy (Felect), http://alloy.phys.cmu.edu/Felect.html, accessed: 2021-08-30.
  * Warren [1969] B. E. Warren,  _X-ray diffraction_ (Addison-Wesley, Reading, MA, 1969).
  * Jeong  _et al._ [2003] I.-K. Jeong, R. H. Heffner, M. J. Graf, and S. J. L. Billinge, Lattice dynamics and correlated atomic motion from the atomic pair distribution function, Phys. Rev. B 67, 104301 (2003).
  * Gao  _et al._ [2017] M. C. Gao, Z. C., P. Gao, F. Zhang, L. Z. Ouyang, M. Widom, and J. A. Hawk, Thermodynamics of concentrated solid solution alloys, Curr. Opin. Solid State and Mat. Sci. 21, 238 (2017).
  * Allison [1996c] T. Allison, JANAF Thermochemical Tables, NIST Standard Reference Database 13 (1996c), https://janaf.nist.gov/tables/Mo-001.html.
  * Allison [1996d] T. Allison, JANAF Thermochemical Tables, NIST Standard Reference Database 13 (1996d), https://janaf.nist.gov/tables/Nb-001.html.
  * Allison [1996e] T. Allison, JANAF Thermochemical Tables, NIST Standard Reference Database 13 (1996e), https://janaf.nist.gov/tables/Ta-001.html.
  * Allison [1996f] T. Allison, JANAF Thermochemical Tables, NIST Standard Reference Database 13 (1996f), https://janaf.nist.gov/tables/W-001.html.
  * Körmann  _et al._ [2017] F. Körmann, Y. Ikeda, B. Grabowski, and M. H. F. Sluiter, Phonon broadening in high entropy alloys, npj Comp. Mat. 3, 36 (2017).
  * Burgers [1934] W. Burgers, On the process of transition of the cubic-body-centered modification into the hexagonal-close-packed modification of zirconium, Physica 1, 561 (1934).
  * Feng and Widom [2018] B. Feng and M. Widom, Band structure theory of the bcc to hcp Burgers distortion, Phys. Rev. B 98, 174108 (2018).
  * van de Walle  _et al._ [2015] A. van de Walle, Q. Hong, S. Kadkhodaei, and R. Sun, The free energy of mechanically unstable phases, Nat. Comm. 6, 1 (2015).
  * Kadkhodaei  _et al._ [2017] S. Kadkhodaei, Q.-J. Hong, and A. van de Walle, Free energy calculation of mechanically unstable but dynamically stabilized bcc titanium, Phys. Rev. B 95, 064101 (2017).
  * Kadkhodaei and van de Walle [2020] S. Kadkhodaei and A. van de Walle, Software tools for thermodynamic calculation of mechanically unstable phases from first-principles data, Comp. Phys. Comm. 246, 106712 (2020).
  * Mei  _et al._ [2009] Z.-G. Mei, S.-L. Shang, Y. Wang, and Z.-K. Liu, Density-functional study of the thermodynamic properties and the pressure–temperature phase diagram of ti, Phys. Rev. B 80, 104116 (2009).
  * Sangiovanni  _et al._ [2019] D. Sangiovanni, J. Klarbring, D. Smirnova, N. Skripnyak, D. Gambino, M. Mrovec, S. Simak, and I. Abrikosov, Superioniclike diffusion in an elemental crystal: Bcc titanium, Phys. Rev. Lett. 123 (2019).
  * Hellman  _et al._ [2011] O. Hellman, I. Abrikosov, and S. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84, 180301 (2011).
  * Allison [1996g] T. Allison, JANAF Thermochemical Tables, NIST Standard Reference Database 13 (1996g), https://janaf.nist.gov/tables/Ti-001.html.
  * Aurelio  _et al._ [2002] G. Aurelio, A. Fernández Guillermet, G. J. Cuello, and J. Campo, Metastable phases in the Ti-V system: Part i. neutron diffraction study and assessment of structural properties, Met. Mat. Trans. A 33, 1307 (2002).
  * Eppelsheimer and Penman [1950] D. S. Eppelsheimer and R. R. Penman, Accurate Determination of the Lattice of Beta-Titanium at 900° C, Nature 166, 960 (1950).
  * Senkov  _et al._ [2001] O. Senkov, B. Chakoumakos, J. Jonas, and F. Froes, Effect of temperature and hydrogen concentration on the lattice parameter of beta titanium, Mat. Resch. Bulletin 36, 1431 (2001).
  * Spreadborough and Christian [1959] J. Spreadborough and J. W. Christian, The Measurement of the Lattice Expansions and Debye Temperatures of Titanium and Silver by X-ray Methods, Proc. Phys. Soc. 74, 609 (1959).
  * Petry  _et al._ [1991] W. Petry, A. Heiming, J. Trampenau, M. Alba, C. Herzig, H. R. Schober, and G. Vogl, Phonon dispersion of the bcc phase of group-IV metals. i, bcc titanium, Phys. Rev. B 43, 10933 (1991).
  * Pawar and Deshpande [1968] R. R. Pawar and V. T. Deshpande, The anisotropy of the thermal expansion of $\alpha$-titanium, Acta Cryst. A 24, 316 (1968).
  * Roberts [1962] W. Roberts, Preferred orientation and anisotropy in titanium, J. Less-common Metals 4, 345 (1962).
  * Willens [1962] R. H. Willens, Vacuum x‐ray diffractometer for high temperature studies of metals sensitive to contamination by oxygen and nitrogen, Rev. Sci. Inst. 33, 1069 (1962).
  * Allison [1996h] T. Allison, JANAF Thermochemical Tables, NIST Standard Reference Database 13 (1996h), https://janaf.nist.gov/tables/Ti-002.html.
  * Allison [1996i] T. Allison, JANAF Thermochemical Tables, NIST Standard Reference Database 13 (1996i), https://janaf.nist.gov/tables/Ti-003.html.
  * Belonoshko  _et al._ [2017] A. B. Belonoshko, T. Lukinov, J. Fu, J. Zhao, S. Davis, and S. I. Simak, Stabilization of body-centred cubic iron under inner-core conditions, Nat. Geoscience 10, 312 (2017).