# Comparing four hard-sphere approximations
for the low-temperature WCA melting line

Eman Attia  Jeppe C. Dyre  Ulf R. Pedersen<EMAIL_ADDRESS>Glass and Time, IMFUFA,
Department of Science and Environment, Roskilde University, P. O. Box 260,
DK-4000 Roskilde, Denmark

###### Abstract

By combining interface-pinning simulations with numerical integration of the
Clausius-Clapeyron equation we determine accurately the melting-line
coexistence pressure and fluid/crystal densities of the Weeks-Chandler-
Andersen (WCA) system covering four decades of temperature. The data are used
for comparing the melting-line predictions of the Boltzmann, Andersen-Weeks-
Chandler, Barker-Henderson, and Stillinger hard-sphere approximations. The
Andersen-Weeks-Chandler and the Barker-Henderson theories give the most
accurate predictions, and they both work excellently in the zero-temperature
limit for which analytical expressions are derived here.

##  I Introduction

While systems of purely repulsive particles are rarely found in nature, they
provide convenient models for both fluids and solids Hansen and McDonald
(2013). Examples are the inverse-power law (IPL) systems based on a
homogeneous pair potential that varies with distance $r$ as $(r/\sigma)^{-n}$
in which $\sigma$ is a length Hoover  _et al._ (1971, 1972); Heyes and Branka
(2007); Branka and Heyes (2011) and the exponential repulsive (EXP) pair
potential that varies with distance as $\exp(-r/\sigma)$ Bacher  _et al._
(2019, 2018); Pedersen  _et al._ (2019). The oldest and most important purely
repulsive system is that of hard spheres (HS) Alder and Wainwright (1957);
Wood and Jacobson (1957); Alder and Wainwright (1959, 1960), which despite its
simplicity provides a good zeroth-order model of realistic systems with both
repulsive and attractive interactions Hirschfelder  _et al._ (1954); Bernal
(1964); Widom (1967); Barker and Henderson (1976); Weeks  _et al._ (1971);
Chandler  _et al._ (1983). A purely repulsive system has a single fluid phase
and no gas-liquid phase transition. In contrast, the symmetry-breaking liquid-
solid transition is present in all purely repulsive systems and, because of
the absence of a gas phase, here the liquid-solid phase boundary extends to
zero temperature.

This paper studies the noted Weeks, Chandler, and Andersen (WCA) purely
repulsive system Weeks  _et al._ (1971); Andrews (1976); Chandler  _et al._
(1983); Speedy  _et al._ (1989); Ben-Amotz and Herschbach (1990); Hess  _et
al._ (1998); Ben-Amotz and Stell (2004a, b); Heyes and Okumura (2006a);
Berthier and Tarjus (2009); Pedersen  _et al._ (2010); Berthier and Tarjus
(2011); Khrapak  _et al._ (2011); Bøhling  _et al._ (2012); Dell and Schweizer
(2015); Valdès  _et al._ (2018); Chattoraj and Ciamarra (2020); Zhou  _et al._
(2020); Attia  _et al._ (2021); Khali  _et al._ (2021); Banerjee and Wales
(2021); Toxvaerd (2021); Singh and Singh (2021); Heyes  _et al._ (2021); Zhou
_et al._ (2022); de Kuijper  _et al._ (1990); Ahmed and Sadus (2009);
Mirzaeinia  _et al._ (2017), which is arrived at by cutting and shifting the
Lennard-Jones (LJ) interaction at its minimum Weeks  _et al._ (1971). In
contrast to the IPL and EXP systems, the WCA pair potential has a finite range
beyond which pair forces are zero, like those of the HS system. At the cutoff,
the WCA pair potential and pair forces are smooth, and at low temperatures one
expects HS approximations to apply because only insignificant “overlaps” are
possible. Thus studies of the low-temperature melting line of the WCA system
provides an excellent testing ground for comparing different HS
approximations. This motivates the present study. In Sec. II we introduce the
WCA system and the four HS approximations considered, as well as give a few
simulation details. Section III details how we determined the WCA melting line
by interface pinning and Clausius-Clapeyron integration. The predictions of
the different HS approximations in regard to pressure and fluid/solid
densities at melting are compared in Sec. IV. Finally, Sec. V provides a brief
outlook.

##  II The WCA system and hard-sphere approximations

###  II.1 The WCA system

We consider mono-disperse systems. Let ${\bf R}=({\bf r}_{1},{\bf
r}_{2},\ldots,{\bf r}_{N})$ be the collective coordinate vector of $N$
particles with mass $m$ confined to the volume $V$ (with periodic boundaries)
and define the number density by $\rho\equiv N/V$. The potential part
$\mathcal{U}({\bf R})$ of the Hamiltonian, $\mathcal{H}({\bf
R})=\mathcal{U}({\bf R})+\mathcal{K}({\bf\dot{R}})$, is assumed to be a sum of
pair-potential contributions,

$\mathcal{U}({\bf R})=\sum_{n>m}^{N}v(|{\bf r}_{m}-{\bf r}_{n}|)\,.$ (1)

We recall that the LJ pair potential is defined Lennard-Jones (1924a, b) by

$v(r)\equiv 4\varepsilon\left[(r/\sigma)^{-12}-(r/\sigma)^{-6}\right]\,$ (2)

in which $\varepsilon$ has units of energy and $\sigma$ units of length. The
WCA pair potential (Fig. 1) is defined by cutting and shifting the LJ
potential at its minimum, which leads to Weeks  _et al._ (1971)

$v(r)=4\varepsilon\left[(r/\sigma)^{-12}-(r/\sigma)^{-6}]+1/4\right]\text{ for
}r\leq r_{c}$ (3)

and zero otherwise where

$r_{c}=\sqrt[6]{2}\sigma\simeq 1.1225\sigma$ (4)

The WCA pair potential is purely repulsive since the pair force $-dv/dr\geq 0$
for all $r$’s, and it is smooth since both $v(r)$ and its first derivative are
continuous (the second derivative is discontinuous at $r_{c}$, though).

Figure 1:  (a) The WCA pair potential Eq. (3). (b) The same pair potential on
a logarithmic energy scale, showing a steep slope at low pair energies.

All quantities obtained by simulations are below reported in units derived
from $m$, $\sigma$, $\varepsilon$ and the Boltzmann constant $k_{B}$.

Simulations of the WCA system are conducted using the RUMD software package
version 3.5 Bailey  _et al._ (2017). An initial configuration is first
constructed by replicating $8\times 8\times 20$ face centered cubic (FCC) unit
cells, resulting in a system of $N=5120$ particles. This initial configuration
is then scaled uniformly to the desired density $\rho$. If a liquid
configuration is needed, the crystal is melted in a high-temperature
simulation. The Newtonian equations of motion are discretized using the leap-
frog algorithm Frenkel and Smit (2002) with the temperature-dependent time
step

$dt=0.001\,\frac{\sigma}{\sqrt{k_{B}T/m}}\,.$ (5)

Simulations in the $NVT$ ensemble Bailey  _et al._ (2017); Frenkel and Smit
(2002); Grønbech-Jensen (2019); Grønbech-Jensen and Farago (2014); Grønbech-
Jensen  _et al._ (2014) are realized using a Langevin thermostat with
relaxation time $t_{T}=m/\alpha$ where $\alpha$ is a friction coefficient
given by

$t_{T}=0.2\,\frac{\sigma}{\sqrt{k_{B}T/m}}\,.$ (6)

For $Np_{z}T$ Langevin simulations Bailey  _et al._ (2017); Grønbech-Jensen
and Farago (2014); Grønbech-Jensen  _et al._ (2014) we used the same
thermostat relaxation time and the barostat relaxation time

$t_{p}=100\,\frac{\sigma}{\sqrt{k_{B}T/m}}\,.$ (7)

We have found that introducing this $1/\sqrt{T}$ scaling to the relaxation
times Heyes and Okumura (2006b) provides a simple way to ensure stability and
efficiency of computations spanning several orders of magnitude in temperature
(see Ref. Ahmed and Sadus, 2009 for a different approach). Note that in this
way the average number of steps needed to travel the distance $\sigma$ for a
thermal particle is the same for all temperatures.

###  II.2 Hard-Sphere approximations to the WCA system

Perturbation theories have proven successful for describing many fluids near
freezing Zwanzig (1954); Rowlinson (1964); Frisch  _et al._ (1966); Barker and
Henderson (1967); Widom (1967); Weeks  _et al._ (1971); Andersen  _et al._
(1971); Verlet and Weis (1972); Barker and Henderson (1976); Andrews (1976);
Chandler  _et al._ (1983); Lado (1984); Kahl (1989); Speedy  _et al._ (1989);
Ben-Amotz and Herschbach (1990); Chang and Sandler (1994); Hess  _et al._
(1998); Hansen and McDonald (2013); Ben-Amotz and Stell (2004a, b);
Trokhymchuk  _et al._ (2005); Heyes and Okumura (2006a); Henderson (2009);
Dyre (2016); Solana (2019); Akhouri and Solana (2022); van Westen  _et al._
(2022). The basic assumption is that the pair interaction can be written as

$v(r)=v_{0}(r)+v_{1}(r)$ (8)

in which $v_{0}(r)$ is the pair potential of a (well-known) reference system
and $v_{1}(r)$ is a small perturbation potential. Often, the HS system is used
as the reference. Several suggestions have been made for choosing the
appropriate HS diameter, $d$. Below we list four well-known HS criteria; in
Sec. IV these are evaluated with respect to their ability to locate the solid-
liquid coexistence line.

In the zero-temperature limit ($T\rightarrow 0$) the WCA pair potential
approaches that of a HS Alder and Wainwright (1957); Wood and Jacobson (1957);
Alder and Wainwright (1959, 1960) system with diameter $d=r_{c}$, i.e., the
system described by

$v_{d}(r)=\infty\text{ for }r<d$ (9)

and zero otherwise. While this may not be intuitively obvious since the WCA
pair potential goes smoothly to zero at the cutoff, it becomes clear when the
WCA potential is shown in a log-plot (Fig. 1(b)). The simplest way of
assigning an effective HS diameter to a WCA particle is to use the truncation
distance

$d=r_{c}\,.$ (10)

This criterion is exact for $T\rightarrow 0$. At finite temperatures, however,
the effective HS diameter will be smaller, and here one needs to make some
physical assumption in order to improve Eq. (10) to arrive at better
approximations. We list below four such approximations.

####  II.2.1 Boltzmann’s hard-sphere criterion

In his 1890 _Lectures on Gas Theory_ Boltzmann (1864) Boltzmann suggested that
the effective HS diameter $d$ should be identified with the distance of
closest approach when the velocities of two head-on colliding particles are
given by their average kinetic energy at far distances. This criterion can be
written as

$v(d)=k_{B}T\,,$ (11)

which for the WCA system results in

$d=\frac{r_{c}}{\sqrt[6]{1+\sqrt{k_{B}T/\varepsilon}}}\,.$ (12)

Boltzmann’s idea, which provides the simplest HS approximation, has been used
to estimate the effective HS diameter of the WCA system by a number of authors
Andrews (1976); Speedy  _et al._ (1989); Ben-Amotz and Herschbach (1990); Hess
_et al._ (1998); Ben-Amotz and Stell (2004b); Heyes and Okumura (2006a).

####  II.2.2 The Andersen-Weeks-Chandler hard-sphere criterion

A more sophisticated HS criterion was suggested in 1971 by Andersen, Weeks,
and Chandler (AWC) Andersen  _et al._ (1971). Their motivation was to match as
well as possible the Helmholtz free energy of the pair potential in question
to the associated HS system. The AWC criterion may be summarized as follows.
If

$e(r)=\exp(-v(r)/k_{B}T)$ (13)

is the pair-potential Boltzmann probability factor, the AWC effective HS
diameter $d$ is identified from

$\int_{0}^{\infty}r^{2}y_{d}(r)\Delta e(r)dr=0$ (14)

in which $\Delta e(r)=e(r)-e_{d}(r)$ is the blip function and $y_{d}(r)$ the
cavity function of the HS fluid. In the Percus-Yevic approximation the cavity
function is given analytically Wertheim (1963); Kahl (1989); Chang and Sandler
(1994); Hansen and McDonald (2013); Trokhymchuk  _et al._ (2005); Smith  _et
al._ (2008); Henderson (2009), which is convenient for applications of Eq.
(14). The appearance of the blip function in Eq. (14) effectively limits the
AWC integral to values near $d$. Thus it is sufficient to consider the zeroth
and first shell of $y_{d}(r)$ to evaluate the AWC integral of Eq. (14) with a
high accuracy. We used the following implementation of the cavity function in
the determination of the HS diameter $d$ via Eq. (14) Chang and Sandler
(1994). If $s\equiv r/d$,

$\displaystyle y_{d}(s)=\begin{cases}c_{0}-c_{1}s+c_{3}s^{3}&\textrm{ for
}s<1\\\ H_{1}(s)/s&\textrm{ for }1<s<2\end{cases}$ (15)

where

$\displaystyle H_{1}(s)$ $\displaystyle=$ $\displaystyle a_{1}\exp A(s-1)r$
$\displaystyle+$ $\displaystyle a_{2}\exp B(s-1)\cos C(s-1)$ $\displaystyle+$
$\displaystyle a_{3}\exp B(s-1)\sin C(s-1)\,.$

The parameters depend on the packing fraction $\eta$ (see Eqs. (6) and
(15)-(17) in Ref. Chang and Sandler, 1994). For the coexistence packing
fraction $\eta=0.4909$, corresponding to the density
$\rho_{l}=0.9375\sigma^{-3}$, we have $c_{0}=58.4514$, $c_{1}=67.9928$,
$c_{3}=14.3461$, $A=1.58498$, $B=-3.68494$, $C=3.85160$ $a_{1}=0.56770$,
$a_{2}=4.23705$ and $a_{3}=-1.41141$. We evaluated the AWC integral
numerically using the Python module SciPy’s et al. (2020) implementation of
QUADPACK Piessens  _et al._ (1983).

####  II.2.3 The Barker- Henderson hard-sphere criterion

The Barker and Henderson (BH) theory Barker and Henderson (1967), which
predates the AWC theory, can be viewed as a simplification of the AWC theory
Hansen and McDonald (2013). Specifically, it is assumed that $r$-squared times
the cavity-function is a constant, $r^{2}y_{d}=$const., implying that Eq. (14)
can be written

$0=\int_{0}^{\infty}[1-e(r)]-[1-e_{d}(r)]dr\,.$ (17)

Since the integral of $1-e_{d}(r)$ is $d$, one arrives at the following HS
criterion

$d=\int_{0}^{\infty}[1-e(r)]dr\,.$ (18)

The $r^{2}y_{d}=$const. assumption is reasonable since the blip function
limits the integral to values near $d$ where $y_{d}$ does not change much when
the temperature is sufficiently low. As $T$ is lowered, the blip function
narrows; thus the AWC diameter reduces to the BH criterion when $T\rightarrow
0$. Note that the BH criteria depends on temperature but not on density (the
AWC criteria depends on both temperature and density). The BH integral of Eq.
(18) is easily evaluated numerically using, e.g., the Python module SciPy’s et
al. (2020) implementation of QUADPACK Piessens  _et al._ (1983).

####  II.2.4 Stillinger’s hard-sphere criterion

At low temperatures, the integrand of the BH criterion Eq. (18) changes
rapidly from nearly unity for $r<d$ to nearly zero for $r>d$. This motivates
the HS criterion proposed by Stillinger in 1976 Stillinger (1976); Likos
(2001); Heyes  _et al._ (2021). He pragmatically identified the HS diameter
from where the pair-potential Boltzmann factor equals one half, i.e.,

$e(d)=\frac{1}{2}\,.$ (19)

Stillinger introduced this in connection with his study of the Gaussian-core
model Stillinger (1976). The same idea can also be applied to the WCA
potential, however, leading Heyes  _et al._ (2021) to

$d=\frac{r_{c}}{\sqrt[6]{1+\sqrt{k_{B}T\ln(2)/\varepsilon}}}\,.$ (20)

The functional form of this HS criterion is identical to that of Boltzmann if
$T$ is replaced by $T\ln(2)$: The factor 2 is here from Eq. (19), and with
$e(d)=1/\exp(1)$ one arrives at Boltzmann’s criterion Ben-Amotz and Stell
(2004b).

Figure 2:  Determination of the coexistence pressure at the temperature $T_{0}=2\varepsilon/k_{B}$ (red diamond) by means of the interface-pinning method Pedersen (2013); Pedersen  _et al._ (2013); Thapar and Escobedo (2014); Pedersen  _et al._ (2015); Cheng  _et al._ (2019); Newman  _et al._ (2019); Steenbergen  _et al._ (2017); Sharma and Escobedo (2018); Zou  _et al._ (2020); Zhu  _et al._ (2021). See the text for details. The inset shows an interface-pinned configuration where the colors indicate the rotational bond order parameter $\bar{q}_{4}$ defined in Ref. Lechner and Dellago, 2008. With this coloring crystalline particles are reddish, while fluid particles are greenish. Table 1:  Selected state points on the coexistence line determined with the interface pinning (IP) method and by numerical integration of the Clausius-Clapeyron (CC) relation (the Supplementary Material gives all computed data). The numbers in parenthesis give the statistical uncertainty (95% confidence interval). $T\>[\varepsilon/k]$ | $p\>[\varepsilon/\sigma^{3}]$ | $\rho_{l}\>[1/\sigma^{3}]$ | $\rho_{s}\>[1/\sigma^{3}]$ | Method
---|---|---|---|---
$20$ | $633.309$ | $1.78328$ | $1.85850$ | CC
$2$ | $31.8086(66)$ | $1.08441(5)$ | $1.15192(6)$ | IP
$2$ | $31.7532$ | $1.08413$ | $1.15163$ | CC
$0.2$ | $2.05169(33)$ | $0.80004(3)$ | $0.87356(4)$ | IP
$0.2$ | $2.05118$ | $0.79992$ | $0.87358$ | CC
$0.02$ | $0.174944(47)$ | $0.70638(5)$ | $0.77889(6)$ | IP
$0.002$ | $0.016687(3)$ | $0.67717(3)$ | $0.74791(3)$ | IP
$0.002$ | $0.016680$ | $0.67705$ | $0.74792$ | CC
Figure 3:  Coexistence pressure as a function of the temperature. (a) The
solid black line shows the reduced coexistence pressure $p/k_{B}T$ as a
function of the temperature (this study). The black dashed line gives the
$T\rightarrow 0$ HS limit, $p_{\bullet}/k_{B}T$, and the colored dots
represent literature coexistence pressures de Kuijper  _et al._ (1990); Ahmed
and Sadus (2009); Mirzaeinia  _et al._ (2017). The red diamonds were computed
with the interface-pinning method (this study). The blue dashed line shows
that at low temperatures the pressure scales as $T^{3/2}$, as expected from HS
theories (see the text). (b) The absolute value of the coexistence pressure in
excess of its $T\rightarrow 0$ limit. Figure 4:  Fluid density at freezing and
solid density at melting as functions of the temperature. (a) The solid black
line shows the density of the fluid at coexistence (this study). The dashed
line marks the $T\rightarrow 0$ limit and the colored dots are literature data
de Kuijper  _et al._ (1990); Ahmed and Sadus (2009); Mirzaeinia  _et al._
(2017). The red diamonds are densities computed with the interface-pinning
method. (b) The solid black line shows the density of the solid at coexistence
(this study), the dashed line is the $T\rightarrow 0$ limit, and the colored
dots represent literature data de Kuijper  _et al._ (1990); Ahmed and Sadus
(2009); Mirzaeinia  _et al._ (2017).

##  III Numerical determination of the phase transition line

The interface pinning method Pedersen (2013); Pedersen  _et al._ (2013);
Thapar and Escobedo (2014); Pedersen  _et al._ (2015); Cheng  _et al._ (2019);
Newman  _et al._ (2019); Steenbergen  _et al._ (2017); Sharma and Escobedo
(2018); Zou  _et al._ (2020); Zhu  _et al._ (2021) is used to compute the
solid-liquid chemical potential difference $\Delta\mu$ for isothermal state-
points at temperatures 0.002$\varepsilon/k_{B}$, 0.02$\varepsilon/k_{B}$,
0.2$\varepsilon/k_{B}$, and 2$\varepsilon/k_{B}$. For a given temperature, we
first set up a FCC crystal elongated in the $z$-direction with the given
density and compute the equilibrium pressure in an $NVT$ simulation. We then
construct a configuration of half-crystal and half-fluid by a high-temperature
simulation, where particle positions are only updated for half of the
particles (resulting in melting for these particles). This produces a
configuration similar to the one shown in the inset of Fig. 2. We then perform
an $Np_{z}T$ simulation by adding a harmonic bias-field to the potential part
of the Hamiltonian,

$\mathcal{U}_{\textrm{IP}}({\bf R})=\mathcal{U}({\bf
R})+\frac{\kappa}{2}\left(Q({\bf R})-a\right)^{2},$ (21)

which forces the system toward configurations with a fluid-crystal interface.
Here, $\kappa$ and $a$ are parameters of the bias-field, and $Q({\bf R})$ is
an order-parameter that measure crystallinity though long-range order (see Eq.
(15) in Ref. Pedersen, 2013). The chemical potential difference between the
two phases, $\Delta\mu$, is computed from the average force, $\kappa(\langle
Q({\bf R})\rangle-a)$, which the bias field results in on the system (see Eq.
(9) in Ref. Pedersen, 2013). This is then repeated for several FCC densities
(and thus pressures) near coexistence. As an example, Fig. 2 shows the
pressures versus the computed chemical potentials at $2\varepsilon/k_{B}$,
considering eleven pressures slightly above $31.7\varepsilon/\sigma^{3}$. The
coexistence state point at $\Delta\mu=0$ is determined by linear regression,
compare the solid line on Fig. 2. From this we find the coexistence pressure
$p=31.8086(66)\varepsilon/\sigma^{3}$ where the number in parenthesis gives
the statistical error on the last two digits using a 95% confidence interval.
Table 1 reports the thermodynamic coexistence data obtained by the interface-
pinning (IP) method and numerical integration of the Clausius-Clapeyron (CC)
relation as detailed below.

While the interface-pinning method is accurate and provides specific error
estimates, it can be computationally expensive because long simulations are
needed to properly represent interface fluctuations, which are usually
significantly slower than fluctuations of the bulk solid and fluid Pedersen
(2013). As an alternative, we determine most points on the coexistence line by
numerical integration of the Clausius-Clapeyron relation (below $s$ and $v$
are the entropy and volume per particle)

$\frac{dp}{dT}=\frac{\Delta s}{\Delta v}\,.$ (22)

This is an example of the Gibbs-Duhem integration methods discussed by Kofke
Kofke (1993a, b), which do not involve slow fluctuations of an interface. The
volume difference $\Delta v=v_{l}-v_{s}$ and the entropy difference $\Delta
s=s_{l}-s_{s}=(\Delta u+p\Delta v-\Delta\mu)/T$ can both be evaluated from
standard $NpT$ simulations of the two bulk phases at coexistence (where of
course $\Delta\mu=0$).

We use a trapezoidal predictor-corrector method to compute coexistence
pressures at the temperatures $T_{i}=0.02\times 10^{(i/24)}$ where $i$ is an
integer, compare the solid black line on Fig. 3. Substituting $t=T$ and $y=p$
we can write the first-order differential equation to be solved in the
standard form

$y^{\prime}=f(t,y)$ (23)

where $f$ is the slope evaluated as $\Delta s/\Delta v$ (Eq.(22)). Suppose one
knows the point $(t_{i},y_{i})$ on the coexistence line, either from the
interface-pinning method or from a previous step of the Clausius-Clapeyron
integration, and wish to compute the next point $(t_{i+1},y_{i+1})$. Let the
step length in $t$ be $h=t_{i+1}-t_{i}$. The prediction of the simple Euler
algorithm is

$y^{(0)}_{i+1}=y_{i}+hf(t_{i},y_{i})\,.$ (24)

A better estimate is provided using Heun’s method:

$y^{(1)}_{i+1}=y_{i}+\frac{h}{2}[f(t_{i},y_{i})+f(t_{i}+h,y^{0}_{i+1})]$ (25)

The next estimate in an iterative predictor-corrector approach is

$y^{(2)}_{i+1}=y_{i}+\frac{h}{2}[f(t_{i},y_{i})+f(t_{i}+h,y^{1}_{i+1})]$ (26)

or, in general,

$y^{k+1}_{i+1}=y_{i}+\frac{h}{2}[f(t_{i},y_{i})+f(t_{i}+h,y^{k}_{i+1})]\,.$
(27)

In the limit of large $k$’s this converges to the trapezoidal rule of
integration where forward and backward integrations yield the same result.

Which criterion to use in order to determine when the predictor-corrector
iterations have converged? To answer this we note that since slopes are
evaluated from finite $NpT$ simulations, we expect a significant statistical
error on the $f$’s used above. If $\bar{f}(t,y)$ is the theoretical slope,
$f(t,y)=\bar{f}(t,y)+e_{f}$ where $e_{f}$ is drawn from a normal distribution
with standard deviation $\sigma_{f}$. This error is estimated by dividing
$NpT$ simulations into statistically independent blocks Flyvbjerg and Petersen
(1989). The error on $y^{k+1}_{i+1}$ is $e_{y}=he_{f}$ and
$\sigma_{y}=|h|\sigma_{f}$. It is sensible to terminate the predictor-
corrector iteration when

$|y^{k+1}_{i+1}-y^{k}_{i+1}|<\sigma_{y}\,$ (28)

since this indicates that changes of $y_{i+1}$’s are mainly due to the
statistical uncertainty on the slopes.

In summary, numerical integration of the Clausius-Clapeyron relation comes
with errors from ignoring higher-order terms and from the statistical
uncertainty of the slopes. To quantify the overall error of the integration we
can compare to the accurate estimates from interface pinning at selected state
points. As an example, for $T_{48}=2\varepsilon/k_{B}$ from the Clausius-
Clapeyron integration we estimate the coexistence pressure to be
$31.7532\varepsilon/\sigma^{3}$, which should be compared to
$31.8086(66)\varepsilon/\sigma^{3}$ for the interface-pinning method, see
Table 1. The error of the computed phase-transition line is not visible in
most figures of this paper, with notable exceptions at low temperatures (error
bars are shown in the below figures whenever errors are significant).

Figures 3 and 4 show coexistence pressures and densities, respectively, from
this study and from the literature de Kuijper  _et al._ (1990); Ahmed and
Sadus (2009); Mirzaeinia  _et al._ (2017). We note that the low-temperature
estimates of Ref. Ahmed and Sadus, 2009 are not accurate, while the high-
temperature estimates of Refs. de Kuijper  _et al._ , 1990; Ahmed and Sadus,
2009; Mirzaeinia  _et al._ , 2017 are consistent with our results. As a
consistency check, we note that the computed coexistence line reaches the HS
limit Pieprzyk  _et al._ (2019) when $T\rightarrow 0$ (the dashed lines on
Figs. 3 and 4 show the HS limits).

Figure 5: Melting-line pressure compared to HS predictions. (a) The solid
black line shows the reduced coexistence pressure, $p/k_{B}T$. The dashed
lines show predictions of the HS theories (see text for details). The red
diamonds show coexistence pressures computed with the interface-pinning
method. (b) $\alpha_{p}(T)=2(p/p_{\bullet}-1)/\sqrt{k_{B}T/\varepsilon}$ (Eq.
(43)) along the computed phase transition line (black solid) and the
theoretical predictions also shown in the upper panel (dashed lines). The blue
dashed line ($\alpha_{0}=0.89(1)$) is the $T\rightarrow 0$ limit determined
from coexistence densities, see Fig. 6(b). AWC and BH give accurate
predictions in the low-temperature limit. The red diamonds are the results of
the interface-pinning method where blue error bars indicate the statistical
error. We note a systematic inaccuracy of the Clausius-Clapeyron integration
(solid black) at the lowest temperatures. Figure 6:  Density-temperature phase
diagram. (a) The solid black lines are the coexistence densities, compare Fig.
4. The vertical black dashed lines mark the $T\rightarrow 0$ HS limits, i.e.,
the quantities $\rho_{\bullet}^{(l)}$ and $\rho_{\bullet}^{(s)}$. The
turquoise, green, yellow, and red dashed curves are predictions of the HS
theories (see the text for details). The two blue dashed lines are the low-
temperature fits
$\rho_{l}=\rho_{\bullet}^{(l)}[1+0.445\sqrt{k_{B}T/\varepsilon}]$ and
$\rho_{s}=\rho_{\bullet}^{(s)}[1+0.445\sqrt{k_{B}T/\varepsilon}]$. (b) The
black $+$’s show
$\alpha_{\rho}(T)=2({\rho}/{\rho_{\bullet}}-1)/\sqrt{k_{B}T/\varepsilon}$
where the densities $\rho$ and $\rho_{\bullet}$ refer to the fluid. The green
$\times$’s is $\alpha_{\rho}(T)$ using the solid densities. The red diamonds
are densities computed with the interface-pinning method. The blue error bars
indicate the 95% confidence interval. We find that the zero-temperature limit
gives $\alpha_{0}=\lim_{T\rightarrow 0}\alpha(T)=0.89(1)$. The turquoise,
green, yellow, and red dashed curves are predictions of the HS theories. The
AWC and BH give the correct low-temperature limit within the statistical
accuracy.

##  IV Comparing the predictions of the different hard-sphere theories

Having accurately located the WCA phase transition, we can now use this to
test the four HS theories by comparing their predictions to the low-
temperature WCA melting-line data.

###  IV.1 Coexistence pressure and densities

Starting with the coexistence pressure, we first need coexistence information
on the HS system. Fernandez et al. Fernandez  _et al._ (2012) estimate that
the HS coexistence pressure is given by $p_{d}=11.5727(10)k_{B}T/d^{3}$. This
value is consistent with

$p_{d}=11.5712(10)\ k_{B}T/d^{3}$ (29)

computed more recently by Pieprzik et al. Pieprzyk  _et al._ (2019); we use
the latter value in this paper. In the zero-temperature limit, the HS diameter
of the WCA interaction is $d=r_{c}$, which gives the coexistence pressure

$p_{\bullet}=8.1821(7)\ k_{B}T/\sigma^{3}\,.\\\ $ (30)

The bullet subscript “$\bullet$” refers throughout the paper to the HS limit
of the WCA model that is approached when $T\rightarrow 0$, i.e., setting
$d=r_{c}$.

The solid black line in Fig. 5(a) shows the coexistence pressure divided by
the thermal energy, $p/k_{B}T$, and the black dashed line shows the $d=r_{c}$
prediction. The predicted pressure is too low since the effective HS diameter
is smaller than $r_{c}$ at finite temperature where particles are allowed to
overlap. In Fig. 5(a) we also consider other criteria for $d$’s (by insertions
into Eq. (29)). At $T=0.02\varepsilon/k_{B}$ the $d=r_{c}$ criterion
underestimates the coexistence pressure by 7%, while both the AWC and BH
criterion give only a 1% error. Thus, the HS theories give a significant
improvement of the predicted coexistence pressure. It is hard to decide from
Fig. 5 which theory is best since this depends on the temperature. We return
below to the low-temperature limit that provides a definite answer. First, we
turn to the HS theories’ predictions of the melting- and freezing densities.

Turning next to the freezing density, we first note that the HS fluid freezing
density has been computed recently by Moir, Lue and Bannerman to the value
Moir  _et al._ (2021)

$\rho_{d}^{(l)}=0.93890(7)/d^{3}$ (31)

and the melting density of the solid to

$\rho_{d}^{(s)}=1.03715(9)/d^{3}$ (32)

In the zero-temperature limit of the WCA system ($d=r_{c}$) we get

$\rho_{\bullet}^{(l)}=0.66390(5)/\sigma^{3}$ (33)

and

$\rho_{\bullet}^{(s)}=0.73337(6)/\sigma^{3}$ (34)

When inserting the $d$’s of the above HS criterions we get the temperature-
dependent density predictions shown in Fig. 6(a) as colored dashed lines.

###  IV.2 Analytical treatment of the low-temperature limit

Inspired by the functional form of Stillingers’s and Boltzmann’s HS criteria
(Eq. (20) and Eq. (12)) we write the low-temperature limit of the effective HS
diameter as

$d_{\alpha}=r_{c}\left(1-\frac{\alpha_{0}}{6}\sqrt{k_{B}T/\varepsilon}\right)\textrm{
for }T\rightarrow 0\,,$ (35)

which implies

$d_{\alpha}^{-3}=r_{c}^{-3}\left(1+\frac{\alpha_{0}}{2}\sqrt{k_{B}T/\varepsilon}\right)\textrm{
for }T\rightarrow 0\,.$ (36)

For the Boltzmann criterion one has $\alpha_{0}=1$ while Stillinger’s
criterion gives $\alpha_{0}=\sqrt{\ln(2)}\simeq 0.83$.

Since $d$ is the same for the AWC and BH criteria in the $T\rightarrow 0$
limit (see Sec. II.2.3), the $\alpha_{0}$’s are also identical. To evaluate
$\alpha_{0}$ we first note that the BH integral defining the HS diameter (Eq.
(18)) can be written

$d=r_{c}-\int_{0}^{r_{c}}\exp(-v(r)/k_{B}T)dr\,.$ (37)

Since the WCA pair potential is purely repulsive, it reaches its minimum at
zero when $r=r_{c}$. Thus at low temperatures the above integral is centered
near $r_{c}$, i.e., near $x=0$ where $x=r_{c}-r$. Keeping the first non-
vanishing term in a Taylor expansion we get

$v(x)=\frac{1}{2}k_{2}x^{2}\textrm{ for }T\rightarrow 0$ (38)

with (see Ref. Attia  _et al._ , 2021)

$k_{2}\equiv\left.\frac{d^{2}v}{dr^{2}}\right|_{r_{c}}=36\sqrt[3]{4}\varepsilon/\sigma^{2}\,.$
(39)

Finding $d$ from Eq. (37) involves solving a Gaussian integral in $x$.
Expanding the upper limit of the integral to infinity, which is exact as
$T\rightarrow 0$, we find

$d=r_{c}-\sqrt{\frac{\pi k_{B}T}{2k_{2}}}\,.$ (40)

By equating $d=d_{\alpha}$ (Eqs. (40) and (35)) we get

$\alpha_{0}=\frac{6}{r_{c}}\sqrt{\frac{\pi\varepsilon}{2k_{2}}}$ (41)

or $\alpha_{0}=\sqrt{\pi}/2\cong 0.886227$. The theoretical $\alpha_{0}$
values are summarized in Table 2.

Table 2:  $\alpha_{0}$ values. From simulations |  | $0.89(1)$
---|---|---
Boltzmann |  | $1$
AWC and BH | $\frac{1}{2}\sqrt{\pi}=$ | $0.886227\ldots$
Stillinger |  $\sqrt{\ln(2)}$ = | $0.832555\ldots$

To estimate $\alpha_{0}$ from the simulations we insert $d_{\alpha}^{-3}$ of
Eq. (36) into Eq. (29) for the coexistence pressure, leading to

$p=p_{\bullet}[1+\frac{\alpha_{0}}{2}\sqrt{k_{B}T/\varepsilon}]\textrm{ for
}T\rightarrow 0\,.$ (42)

Thus, a way to determine $\alpha_{0}$ is to define the function (Fig. 5(b))

$\alpha_{p}(T)=\frac{2}{\sqrt{k_{B}T/\varepsilon}}\left[\frac{p(T)}{p_{\bullet}}-1\right]$
(43)

for which we note that $\alpha_{0}=\alpha(T)$ for $T\rightarrow 0$. Similarly,
we get for the densities $\rho=\rho_{l}$ or $\rho=\rho_{s}$

$\rho=\rho_{\bullet}(1+\frac{\alpha_{0}}{2}\sqrt{k_{B}T/\varepsilon})\textrm{
for }T\rightarrow 0$ (44)

and define

$\alpha_{\rho}(T)=\frac{2}{\sqrt{k_{B}T/\varepsilon}}\left[\frac{\rho(T)}{\rho_{\bullet}}-1\right]\,.$
(45)

Figure 6(a) shows the temperature dependence of the fluid and solid densities
at coexistence (solid lines). These densities yield the $\alpha_{\rho}(T)$’s
shown with black $+$’s and green $x$’s, respectively, on Fig. 6(b). From the
low-temperature points we estimate $\alpha_{0}=0.89(1)$. The colored dashed
lines show the predictions of the HS theories (the $T\rightarrow 0$ limits
agree with the values of Table 2). We conclude that the AWC and BH theories
gives excellent agreement as $T\rightarrow 0$. Figure 5(b) shows
$\alpha_{p}(T)$ computed using the coexistence pressure. In agreement with the
results for the $\alpha_{\rho}(T)$’s we find that $\alpha_{0}=0.89(1)$ (blue
dashed line).

The success of the AWC and BH theories suggests writing the coexistence
pressure and densities as follows (inserting $\alpha_{0}=\sqrt{\pi}/2$ into
Eqs. (42) and (44))

$p=p_{\bullet}\left(1+\sqrt{\frac{\pi k_{B}T}{16\varepsilon}}\right)$ (46)

and

$\rho=\rho_{\bullet}\left(1+\sqrt{\frac{\pi k_{B}T}{16\varepsilon}}\right),$
(47)

respectively, see the blue dashed lines of Figs. 3-6. Interestingly, this low-
temperature approximation gives better predictions than the _neat_ HS theories
– even at high temperatures (with the exception of Boltzmann’s criterion near
$T\simeq 0.5\varepsilon/k_{B}$). We do not have an explanation for this.

Equations (46) and (47) summarize an important result of this paper, providing
an analytical HS approximation for the low-temperature freezing of the WCA
fluid. This can be generalized to any other purely repulsive pair-potential
that is truncated smoothly at $r=r_{c}$ by the following steps:

  1. 1.

Compute $k_{2}$ using Eq. (39) and

  2. 2.

derive $\alpha_{0}$ within the BH theory by inserting $k_{2}$ into Eq. (41).

  3. 3.

Low-temperature predictions for coexistence pressure and densities are then
provided by inserting $\alpha_{0}$ into Eqs. (42) and (44), respectively.

Figure 7:  The solid black line shows the reduced coexistence pressure
$p/k_{B}T$ as a function of the temperature. The red and green dashed curves
are isomorphs of the fluid, i.a., lines along which the excess entropy,
$S_{ex}$, is constant. The isomorphs touch the phase-transition line at
$T_{\star}=0.02\varepsilon/k_{B}$ and $T_{\star}=2\varepsilon/k_{B}$,
respectively. The turquoise dashed line is the prediction of the AWC theory.
Figure 8:  (a) The radial distribution $g(r)$ for the fluid at coexistence.
(b) The radial distribution as a function of the reduced distance
$\tilde{r}=r\sqrt[3]{\rho}$ for the fluid at coexistence. (c) The radial
distribution as a function of the reduced distance $\tilde{r}=r\sqrt[3]{\rho}$
for a fluid isomorph that touch the coexistence line at
$T_{\star}=2\varepsilon/k_{B}$.  Figure 9:  The solid lines show the mean-
square displacement $\langle|{\bf r}(t)-{\bf r}(0)|^{2}\rangle$ for selected
state-points along the coexistence line (see Table 1). The dashed lines are
long-time fits to $\langle|{\bf r}(t)-{\bf r}(0)|^{2}\rangle=6Dt$ where $D$ is
the diffusion coefficient. The dots in the inset show the reduced diffusion
coefficient $\tilde{D}=D\rho^{1/3}T^{1/2}$. The red $x$’s and green $+$’s is
the reduced diffusion coefficient for state points along the isomorphs with
$T_{\star}=2\varepsilon/k:B$ and $T_{\star}=0.02\varepsilon/k_{B}$,
respectively.

##  V Outlook

We have shown that HS theories give excellent predictions of the WCA melting
line at low temperatures, in particular for the AWC and BH approximations. At
higher temperatures the HS theories are less accurate. This is not surprising
because the WCA model only resembles a HS system at low temperatures. How to
predict the WCA melting-line pressures and coexistence densities at high
temperatures? One possibility is to generalize the low-temperature HS
approximation by considering the lines of constant excess entropy
$S_{\textrm{ex}}$ (this is the entropy in excess of the ideal gas entropy at
the same density and temperature, a negative quantity that in some textbooks
Prigogine and Defay (1954) is referred to as the residual entropy). For the HS
system these lines are determined entirely by the density, i.e., they are
vertical in the density-temperature phase diagram. In Ref. Attia  _et al._ ,
2021 it has been shown that the WCA system’s structure and dynamics are near-
invariant along the lines of constant excess entropy, which are referred to as
isomorphs Gnan  _et al._ (2009); Schrøder and Dyre (2014). An isomorph can be
computed by numerical integration in the $\ln T$-$\ln\rho$ plane (e.g., using
the fourth-order Runge-Kutta method (RK4) Attia  _et al._ (2021)) for which
the required slope is $f=1/\gamma$ where Gnan  _et al._ (2009); Dyre (2014,
2018)

$\gamma\equiv\left(\frac{\partial\ln
T}{\partial\ln\rho}\right)_{S_{\textrm{ex}}}\,.$ (48)

The “density-scaling exponent” $\gamma$ may be computed from virial- and
potential-energy fluctuations in the $NVT$ ensemble via the general
statistical-mechanical identity $\gamma=\langle\Delta W\Delta
U\rangle/\langle(\Delta U)^{2}\rangle$ Gnan  _et al._ (2009). Figure 7 shows
the reduced pressure $p/k_{B}T$ of two fluid isomorphs that overlap with the
coexistence line at $T_{\star}=0.02\varepsilon/k_{B}$ and
$T_{\star}=2\varepsilon/k_{B}$, respectively (dashed green and red lines). For
comparison, the turquoise dashed line shows the prediction of the reduced
coexistence pressure of the AWC theory. For the entire temperature span the
isomorphs gives predictions with an overall accuracy comparable to that of the
best HS approximation (AWC).

Figures 8 and 9 show the structure and dynamics along the melting line and the
fluid isomorph in reduced units Gnan  _et al._ (2009). Interestingly, the
physics is more invariant along the coexistence lines than along the isomorph.
This is in contrast to previous findings for the LJ system, where the opposite
applies Pedersen  _et al._ (2016). We note, however, that isomorphs only
follow the coexistence lines to a first approximation. For the LJ system,
accurate predictions for the thermodynamics of freezing and melting can be
arrived at within the isomorph-theoretical perturbation framework proposed in
Ref. Pedersen  _et al._ , 2016 – we hope to apply the same method to the WCA
system in the near future.

## Acknowledgment

This work was supported by the VILLUM Foundation’s Matter grant (No. 16515).

## Author Declarations

### Conflict of Interest

The authors have no conflicts to disclose.

## Data availability

The data that support the findings of this study are openly available in
Zenodo at http://doi.org/10.5281/zenodo.6505218, reference number 6505218.

## References

  * Hansen and McDonald (2013) J.-P. Hansen and I. R. McDonald,  _Theory of Simple Liquids: With Applications to Soft Matter_ , 4th ed. (Academic, New York, 2013).
  * Hoover  _et al._ (1971) W. G. Hoover, S. G. Gray, and K. W. Johnson, “Thermodynamic properties of the fluid and solid phases for inverse power potentials,” J. Chem. Phys. 55, 1128–1136 (1971).
  * Hoover  _et al._ (1972) W. G. Hoover, D. A. Young, and R. Grover, “Statistical mechanics of phase diagrams. i. inverse power potentials and the close-packed to body-packed cubic transition,” J. Chem. Phys. 56, 2207–2210 (1972).
  * Heyes and Branka (2007) D. M. Heyes and A. C. Branka, “Physical properties of soft repulsive particle fluids,” Phys. Chem. Chem. Phys. 9, 5570–5575 (2007).
  * Branka and Heyes (2011) A. C. Branka and D. M. Heyes, “Pair correlation function of soft-sphere fluids,” J. Chem. Phys. 134, 064115 (2011).
  * Bacher  _et al._ (2019) A. K. Bacher, T. B. Schrøder, and J. C. Dyre, “The EXP pair-potential system. I. Fluid phase isotherms, isochores, and quasiuniversality,” J. Chem. Phys. 149, 114501 (2019).
  * Bacher  _et al._ (2018) A. K. Bacher, T. B. Schrøder, and J. C. Dyre, “The EXP pair-potential system. II. Fluid phase isomorphs,” J. Chem. Phys. 149, 114502 (2018).
  * Pedersen  _et al._ (2019) U. R. Pedersen, A. K. Bacher, T. B. Schrøder, and J. C. Dyre, “The EXP pair-potential system. III. Thermodynamic phase diagram,” J. Chem. Phys. 150, 174501 (2019).
  * Alder and Wainwright (1957) B. J. Alder and T. E. Wainwright, “Phase transition for a hard sphere system,” J. Chem. Phys. 27, 1208–1209 (1957).
  * Wood and Jacobson (1957) W. W. Wood and J. D. Jacobson, “Preliminary results from a recalculation of the Monte Carlo equation of state of hard spheres,” J. Chem. Phys. 27, 1207–1208 (1957).
  * Alder and Wainwright (1959) B. J. Alder and T. E. Wainwright, “Studies in molecular dynamics. I. general method,” J. Chem. Phys. 31, 459–466 (1959).
  * Alder and Wainwright (1960) B. J. Alder and T. E. Wainwright, “Studies in molecular dynamics. II. behavior of a small number of elastic spheres,” J. Chem. Phys. 33, 1439–1451 (1960).
  * Hirschfelder  _et al._ (1954) J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird,  _Molecular Theory of Gases and Liquids_ (J. Wiley & Sons (New York), 1954).
  * Bernal (1964) J. D. Bernal, “The Bakerian lecture, 1962. The structure of liquids,” Proc. R. Soc. London Ser. A 280, 299–322 (1964).
  * Widom (1967) B. Widom, “Intermolecular forces and the nature of the liquid state,” Science 157, 375–382 (1967).
  * Barker and Henderson (1976) J. A. Barker and D. Henderson, “What is ”liquid”? Understanding the states of matter,” Rev. Mod. Phys. 48, 587–671 (1976).
  * Weeks  _et al._ (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, “Role of repulsive forces in determining the equilibrium structure of simple liquids,” J. Chem. Phys. 54, 5237–5247 (1971).
  * Chandler  _et al._ (1983) D. Chandler, J. D. Weeks, and H. C. Andersen, “Van der Waals picture of liquids, solids, and phase transformations,” Science 220, 787–794 (1983).
  * Andrews (1976) F. C. Andrews, “A simple approach to the equilibrium statistical mechanics of two-dimensional fluids,” J. Chem. Phys. 64, 1941–1947 (1976).
  * Speedy  _et al._ (1989) R. J. Speedy, F. X. Prielmeier, T. Vardag, E. W. Lang, and H.-D. Lüdemann, “Diffusion in simple fluids,” Molecular Physics 66, 577–590 (1989).
  * Ben-Amotz and Herschbach (1990) D. Ben-Amotz and D. R. Herschbach, “Estimation of effective diameters for molecular fluids,” J. Phys. Chem. 94, 1038–1047 (1990).
  * Hess  _et al._ (1998) S. Hess, M. Kröger, and H. Voigt, “Thermomechanical properties of the WCA-Lennard-Jones model system in its fluid and solid states,” Physica A 250, 58–82 (1998).
  * Ben-Amotz and Stell (2004a) D. Ben-Amotz and G. Stell, “Hard sphere perturbation theory for fluids with soft-repulsive-core potentials,” J. Chem. Phys. 120, 4844–4851 (2004a).
  * Ben-Amotz and Stell (2004b) D. Ben-Amotz and G. Stell, “Reformulation of weeks-chandler-andersen perturbation theory directly in terms of a hard-sphere reference system,” J. Phys. Chem. B 108, 6877–6882 (2004b).
  * Heyes and Okumura (2006a) D. M. Heyes and H. Okumura, “Equation of state and structural properties of the Weeks-Chandler-Andersen fluid,” J. Chem. Phys. 124, 164507 (2006a).
  * Berthier and Tarjus (2009) L. Berthier and G. Tarjus, “Nonperturbative effect of attractive forces in viscous liquids,” Phys. Rev. Lett. 103 (2009), 10.1103/physrevlett.103.170601.
  * Pedersen  _et al._ (2010) U. R. Pedersen, T. B. Schrøder, and J. C. Dyre, “Repulsive reference potential reproducing the dynamics of a liquid with attractions,” Phys. Rev. Lett. 105 (2010), 10.1103/physrevlett.105.157801.
  * Berthier and Tarjus (2011) L. Berthier and G. Tarjus, “The role of attractive forces in viscous liquids,” J. Chem. Phys. 134, 214503 (2011).
  * Khrapak  _et al._ (2011) S. A. Khrapak, M. Chaudhuri, and G. E. Morfill, “Communication: Universality of the melting curves for a wide range of interaction potentials,” J. Chem. Phys. 134, 241101 (2011).
  * Bøhling  _et al._ (2012) L. Bøhling, A. A. Veldhorst, T. S Ingebrigtsen, N. P. Bailey, J. S. H.en, S. Toxvaerd, T. B. Schrøder, and J. C Dyre, “Do the repulsive and attractive pair forces play separate roles for the physics of liquids?” J. Condens. Matter Phys. 25, 032101 (2012).
  * Dell and Schweizer (2015) Z. E. Dell and K. S. Schweizer, “Microscopic theory for the role of attractive forces in the dynamics of supercooled liquids,” Phys. Rev. Lett. 115 (2015), 10.1103/physrevlett.115.205702.
  * Valdès  _et al._ (2018) L.-C. Valdès, J. Gerges, T. Mizuguchi, and F. Affouard, “Crystallization tendencies of modelled Lennard-Jones liquids with different attractions,” J. Chem. Phys. 148, 014501 (2018).
  * Chattoraj and Ciamarra (2020) J. Chattoraj and M. P. Ciamarra, “Role of attractive forces in the relaxation dynamics of supercooled liquids,” Phys. Rev. Lett. 124 (2020), 10.1103/physrevlett.124.028001.
  * Zhou  _et al._ (2020) Y. Zhou, B. Mei, and K. S. Schweizer, “Integral equation theory of thermodynamics, pair structure, and growing static length scale in metastable hard sphere and Weeks-Chandler-Andersen fluids,” Phys. Rev. E 101 (2020), 10.1103/physreve.101.042121.
  * Attia  _et al._ (2021) E. Attia, J. C. Dyre, and U. R. Pedersen, “Extreme case of density scalig: The Weeks-Chandler-Andersen system at low temperatures,” Phys. Rev. E 103 (2021), 10.1103/physreve.103.062140.
  * Khali  _et al._ (2021) S. S. Khali, D. Chakraborty, and D. Chaudhuri, “Two-step melting of the Weeks-Chandler-Anderson system in two dimensions,” Soft Mat. 17, 3473–3485 (2021).
  * Banerjee and Wales (2021) A. Banerjee and D. J. Wales, “Energy landscapes for a modified repulsive Weeks-Chandler-Andersen potential,” J. Condens. Matter Phys. 34, 034004 (2021).
  * Toxvaerd (2021) S. Toxvaerd, “Role of the attractive forces in a supercooled liquid,” Phys. Rev. E 103 (2021), 10.1103/physreve.103.022611.
  * Singh and Singh (2021) A. Singh and Y. Singh, “How attractive and repulsive interactions affect structure ordering and dynamics of glass-forming liquids,” Phys. Rev. E 103 (2021), 10.1103/physreve.103.052105.
  * Heyes  _et al._ (2021) D. M. Heyes, S. Pieprzyk, and A. C. Brańka, “Application of cell models to the melting and sublimation lines of the Lennard-Jones and related potential systems,” Phys. Rev. E 104 (2021), 10.1103/physreve.104.044119.
  * Zhou  _et al._ (2022) Y. Zhou, B. Mei, and K. S. Schweizer, “Activated relaxation in supercooled monodisperse atomic and polymeric WCA fluids: Simulation and ECNLE theory,” J. Chem. Phys. 156, 114901 (2022).
  * de Kuijper  _et al._ (1990) A. de Kuijper, J. A. Schouten, and J. P. J. Michels, “The melting line of the Weeks-Chandler-Anderson Lennard-Jones reference system,” J. Chem. Phys. 93, 3515–3519 (1990).
  * Ahmed and Sadus (2009) A. Ahmed and R. J. Sadus, “Phase diagram of the Weeks-Chandler-Andersen potential from very low to high temperatures and pressures,” Phys. Rev. E 80 (2009), 10.1103/physreve.80.061101.
  * Mirzaeinia  _et al._ (2017) A. Mirzaeinia, F. Feyzi, and S. Majid Hashemianzadeh, “Equation of state and Helmholtz free energy for the atomic system of the repulsive Lennard-Jones particles,” J. Chem. Phys. 147, 214503 (2017).
  * Lennard-Jones (1924a) J. E. Lennard-Jones, “On the determination of molecular fields. – I. from the variation of the viscosity of a gas with temperature,” Proc. R. Soc. Lond. A 106, 441–462 (1924a).
  * Lennard-Jones (1924b) J. E. Lennard-Jones, “On the determination of molecular fields. – II. from the equation of state of a gas,” Proc. R. Soc. Lond. A 106, 463–477 (1924b).
  * Bailey  _et al._ (2017) N. Bailey, T. Ingebrigtsen, J. S. Hansen, A. Veldhorst, L. Bøhling, C. Lemarchand, A. Olsen, A. Bacher, L. Costigliola, U. R. Pedersen, H. Larsen, J. C. Dyre, and T. Schrøder, “RUMD: A general purpose molecular dynamics package optimized to utilize GPU hardware down to a few thousand particles,” SciPost Phys. 3 (2017), 10.21468/scipostphys.3.6.038.
  * Frenkel and Smit (2002) D. Frenkel and B. Smit,  _Understanding molecular simulation : from algorithms to applications_ , 2nd ed. (Academic Press, San Diego, 2002).
  * Grønbech-Jensen (2019) N. Grønbech-Jensen, “Complete set of stochastic Verlet-type thermostats for correct Langevin simulations,” Mol. Phys. 118, e1662506 (2019).
  * Grønbech-Jensen and Farago (2014) N. Grønbech-Jensen and O. Farago, “Constant pressure and temperature discrete-time Langevin molecular dynamics,” J. Chem, Phys. 141, 194108 (2014).
  * Grønbech-Jensen  _et al._ (2014) N. Grønbech-Jensen, N. R. Hayre, and O. Farago, “Application of the G-JF discrete-time thermostat for fast and accurate molecular simulations,” Comput. Phys. Commun. 185, 524–527 (2014).
  * Heyes and Okumura (2006b) D. M. Heyes and H. Okumura, “Equation of state and structural properties of the Weeks-Chandler-Andersen fluid,” J. Chem. Phys. 124, 164507 (2006b).
  * Zwanzig (1954) R. W. Zwanzig, “High-temperature equation of state by a perturbation method. i. nonpolar gases,” J. Chem. Phys. 22, 1420–1426 (1954).
  * Rowlinson (1964) J. S. Rowlinson, “An equation of state of gases at high temperatures and densities,” Mol. Phys. 7, 349–361 (1964).
  * Frisch  _et al._ (1966) H. L. Frisch, J. L. Katz, E. Praestgaard, and J. L. Lebowitz, “High-temperature equation of state – argon,” J. Phys. Chem. 70, 2016–2020 (1966).
  * Barker and Henderson (1967) J. A. Barker and D. Henderson, “Perturbation theory and equation of state for fluids. II. a successful theory of liquids,” J. Chem. Phys. 47, 4714–4721 (1967).
  * Andersen  _et al._ (1971) H. C. Andersen, J. D. Weeks, and D. Chandler, “Relationship between the hard-sphere fluid and fluids with realistic repulsive forces,” Phys. Rev. A 4, 1597–1607 (1971).
  * Verlet and Weis (1972) L. Verlet and J.-J. Weis, “Perturbation theory for the thermodynamic properties of simple liquids,” Mol. Phys. 24, 1013–1024 (1972).
  * Lado (1984) F. Lado, “Choosing the reference system for liquid state perturbation theory,” Mol. Phys. 52, 871–876 (1984).
  * Kahl (1989) G. Kahl, “Analytic representation for the pair-correlation function of a hard-sphere Yukawa system,” Mol. Phys. 67, 879–889 (1989).
  * Chang and Sandler (1994) J. Chang and S. I. Sandler, “A real function representation for the structure of the hard-sphere fluid,” Mol. Phys. 81, 735–744 (1994).
  * Trokhymchuk  _et al._ (2005) A. Trokhymchuk, I. Nezbeda, J. Jirsák, and D. Henderson, “Hard-sphere radial distribution function again,” J. Chem. Phys. 123, 024501 (2005).
  * Henderson (2009) D. Henderson, “Analytic methods for the Percus-Yevick hard sphere correlation functions,” Condens. Matter Phys. 12, 127–135 (2009).
  * Dyre (2016) J. C. Dyre, “Simple liquids’ quasiuniversality and the hard-sphere paradigm,” J. Phys. Condens. Matter 28, 323001 (2016).
  * Solana (2019) J. R. Solana,  _Perturbation theories for the thermodynamic properties of fluids and solids_ (CRC Press, London, England, 2019).
  * Akhouri and Solana (2022) B. P. Akhouri and J. R. Solana, “On the choice of the effective diameter in the high-temperature expansion for the Lennard-Jones fluid,” Mol. Phys. (2022), 10.1080/00268976.2022.2028918.
  * van Westen  _et al._ (2022) T. van Westen, M. Hammer, B. Hafskjold, A. Aasen, J. Gross, and Ø. Wilhelmsen, “Perturbation theories for fluids with short-ranged attractive forces: A case study of the Lennard-Jones spline fluid,” J. Chem. Phys. 156, 104504 (2022).
  * Boltzmann (1864) L. Boltzmann,  _Lectures on Gas Theory_ (Dover Publications, 1864).
  * Wertheim (1963) M. S. Wertheim, “Exact solution of the Percus-Yevick integral equation for hard spheres,” Phys. Rev. Lett. 10, 321–323 (1963).
  * Smith  _et al._ (2008) W. R. Smith, D. J. Henderson, P. J. Leonard, J. A. Barker, and E. W. Grundke, “Fortran codes for the correlation functions of hard sphere fluids,” Mol. Phys. 106, 3–7 (2008).
  * et al. (2020) P. Virtanen et al., “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nat. Methods 17, 261–272 (2020).
  * Piessens  _et al._ (1983) R. Piessens, E. De Doncker-Kapenga, C. Uberhuber, and D. K. Kahaner,  _Quadpack_ , Springer series in computational mathematics (Springer, Berlin, Germany, 1983).
  * Stillinger (1976) F. H. Stillinger, “Phase transitions in the gaussian core system,” J. Chem. Phys. 65, 3968–3974 (1976).
  * Likos (2001) C. N. Likos, “Effective interactions in soft condensed matter physics,” Phys. Rep. 348, 267–439 (2001).
  * Pedersen (2013) U. R. Pedersen, “Direct calculation of the solid-liquid Gibbs free energy difference in a single equilibrium simulation,” J. Chem. Phys. 139, 104102 (2013).
  * Pedersen  _et al._ (2013) U. R. Pedersen, F. Hummel, G. Kresse, G. Kahl, and C. Dellago, “Computing Gibbs free energy differences by interface pinning,” Phys. Rev. B 88 (2013), 10.1103/physrevb.88.094101.
  * Thapar and Escobedo (2014) V. Thapar and F. A. Escobedo, “Extensions of the interfacial pinning method and application to hard core systems,” J. Chem. Phys. 141, 124117 (2014).
  * Pedersen  _et al._ (2015) U. R. Pedersen, F. Hummel, and C. Dellago, “Computing the crystal growth rate by the interface pinning method,” J. Chem. Phys. 142, 044104 (2015).
  * Cheng  _et al._ (2019) B. Cheng, E. A. Engel, J. B., C. Dellago, and M. Ceriotti, “Ab initio thermodynamics of liquid and solid water,” PNAS 116, 1110–1115 (2019).
  * Newman  _et al._ (2019) R. S. Newman, S. N., J. Dshemuchadse, and S. C. Glotzer, “Shape-controlled crystallisation pathways in dense fluids of ccp-forming hard polyhedra,” Mol. Phys. 117, 3819–3826 (2019).
  * Steenbergen  _et al._ (2017) K. G. Steenbergen, E. Pahl, and P. Schwerdtfeger, “Accurate, large-scale density functional melting of Hg: Relativistic effects decrease melting temperature by 160 K,” J. Phys. Chem. Lett. 8, 1407–1412 (2017).
  * Sharma and Escobedo (2018) A. K. Sharma and F. A. Escobedo, “Nucleus-size pinning for determination of nucleation free-energy barriers and nucleus geometry,” J. Chem. Phys. 148, 184104 (2018).
  * Zou  _et al._ (2020) Y. Zou, S. Xiang, and C. Dai, “Investigation on the efficiency and accuracy of methods for calculating melting temperature by molecular dynamics simulation,” Comput. Mater. Sci. 171, 109156 (2020).
  * Zhu  _et al._ (2021) L. Zhu, J. Janssen, S. Ishibashi, F. Körmann, B. Grabowski, and J. Neugebauer, “A fully automated approach to calculate the melting temperature of elemental crystals,” Comput. Mater. Sci. 187, 110065 (2021).
  * Lechner and Dellago (2008) W. Lechner and C. Dellago, “Accurate determination of crystal structures based on averaged local bond order parameters,” J. Chem. Phys. 129, 114707 (2008).
  * Kofke (1993a) D. A. Kofke, “Gibbs-Duhem integration: a new method for direct evaluation of phase coexistence by molecular simulation,” Mol. Phys. 78, 1331–1336 (1993a).
  * Kofke (1993b) D. A. Kofke, “Direct evaluation of phase coexistence by molecular simulation via integration along the saturation line,” J. Chem. Phys. 98, 4149–4162 (1993b).
  * Flyvbjerg and Petersen (1989) H. Flyvbjerg and H. G. Petersen, “Error estimates on averages of correlated data,” J. Chem. Phys. 91, 461–466 (1989).
  * Pieprzyk  _et al._ (2019) S. Pieprzyk, M. N. Bannerman, A. C. Brańka, M. Chudak, and D. M. Heyes, “Thermodynamic and dynamical properties of the hard sphere system revisited by molecular dynamics simulation,” Phys. Chem. Chem. Phys. 21, 6886–6899 (2019).
  * Fernandez  _et al._ (2012) L. A. Fernandez, V. Martin-Mayor, B. Seoane, and P. Verrocchio, “Equilibrium fluid-solid coexistence of hard spheres,” Phys. Rev. Lett. 108 (2012), 10.1103/physrevlett.108.165701.
  * Moir  _et al._ (2021) C. Moir, L. Lue, and M. N. Bannerman, “Tethered-particle model: The calculation of free energies for hard-sphere systems,” J. Chem. Phys. 155, 064504 (2021).
  * Prigogine and Defay (1954) I. Prigogine and R. Defay,  _Chemical thermodynamics_ (Prentice Hall Press, 1954).
  * Gnan  _et al._ (2009) N. Gnan, T. B. Schrøder, U. R. Pedersen, N. P. Bailey, and J. C. Dyre, “Pressure-energy correlations in liquids. IV. ’Isomorphs’ in liquid phase diagrams,” J. Chem. Phys. 131, 234504 (2009).
  * Schrøder and Dyre (2014) T. B. Schrøder and J. C. Dyre, “Simplicity of condensed matter at its core: Generic definition of a Roskilde-simple system,” J. Chem. Phys. 141, 204502 (2014).
  * Dyre (2014) J. C. Dyre, “Hidden scale envariance in condensed matter,” J. Phys. Chem. B 118, 10007–10024 (2014).
  * Dyre (2018) J. C. Dyre, “Perspective: Excess-entropy scaling,” J. Chem. Phys. 149, 210901 (2018).
  * Pedersen  _et al._ (2016) U. R. Pedersen, L. Costigliola, N. P. Bailey, T. B. Schrøder, and J. C. Dyre, “Thermodynamics of freezing and melting,” Nat. Commun. 7 (2016), 10.1038/ncomms12386.