# Effects of Structural Inhomogeneity on Equilibration Processes in Langevin
Dynamics

Omid Mozafar<EMAIL_ADDRESS>Applied Mathematics Department, The University
of Western Ontario, London, Ontario, Canada N6A 5B7  Colin Denniston
<EMAIL_ADDRESS>Physics and Astronomy Department, The University of Western
Ontario, London, Ontario, Canada, N6A 3K7

###### Abstract

In recent decades, computer experiments have led to an accurate and
fundamental understanding of atomic and molecular mechanisms in fluids, such
as different kinds of relaxation processes toward steady physical states. In
this paper, we investigate how exactly the configuration of initial states in
a molecular-dynamics simulation can affect the rates of decay toward
equilibrium for the widely-known Langevin canonical ensemble. For this
purpose, we derive an original expression relating the system relaxation time
$\tau_{sys}$ and the radial distribution function $g(r)$ in the near-zero and
high-density limit. We found that for an initial state which is slightly
marginally inhomogeneous in the number density of atoms, the system relaxation
time $\tau_{sys}$ is much longer than that for the homogeneous case and an
increasing function of the Langevin coupling constant, $\gamma$. We also found
during structural equilibration, $g(r)$ at large distances approaches 1 from
above for the inhomogeneous case and from below for the macroscopically
homogeneous one. Keywords: System Relaxation Time, Molecular-Dynamics
Simulation, Langevin Thermostat, Wallace’s Entropy Expansion, Inhomogeneous
Structure, Radial Distribution Function.

##  I introduction

One of the fundamental physical processes in the world is the relaxation
process of many-body systems from any perturbation. Generally speaking, once a
system relaxes, it becomes temporally invariant and hence one would define
reliable and measurable quantities indicating the systems’ properties. The
word relaxation was originally utilized by the Scottish physicist, James C.
Maxwell, in 1867 to describe molecular processes [1]. Since then, many
researchers have put significant effort into generalizing and explaining the
basic concepts behind relaxation processes to be applicable to a wider range
of phenomena [2, 3, 4, 5, 6, 7, 8, 9]. In particular, in 1946, the Russian
physicist, Yacov Frenkel, shed some light on the relaxation processes in
liquids [10], leading to locating the Frenkel Line in supercritical fluids
[11, 12, 13, 14, 15].

Over the past few decades, the invention of novel methods of conducting
computer-based experiments with the advancement of algorithms and/or computers
has opened a new chapter in the study of different kinds of relaxation
processes [16, 17, 18]. However, computer experiments sometimes produce
nonphysical results due to imperfect algorithms or models, especially those
related to thermostats [19, 20, 21, 22]. In this work, by studying the
temporal evolution of the radial distribution function in molecular-dynamics
(MD) simulations [23, 24, 25, 26] we show possible impacts of initial
structural inhomogeneity on the relaxation processes in Langevin dynamics.

Measurement of the radial distribution function (RDF) in MD simulations (with
a constant number of atoms) has a long history [27, 28, 29, 30, 31]. Early
efforts suffered greatly from the inherently small system sizes accessible to
molecular dynamics. This spawned numerous works on accounting for all kinds of
finite-size effects when measuring the RDF [32, 33, 34, 35, 36]. Recent
advances in computer technology now allows the study of larger and larger
systems, i.e., $N>>1$. An unappreciated effect of this is the much longer
relaxation times needed for long-wavelength structural degrees of freedom. We
study these kinds of effects in this paper.

This paper is organized as follows: Sec. II reviews some of the physical
concepts; Sec. III derives an original expression connecting the RDF and
$\tau_{sys}$; Sec. IV provides details on the computer simulations and system
features; Sec. V presents and discusses the results obtained from our MD
simulations; Sec. VI summarizes the paper along with some conclusions and
suggestions.

##  II Background

The radial distribution function is the normalized local density distribution
within the system when one looks radially outwards from any particle. It can
also be thought of as a measure of the probability of finding one particle of
any shape and orientation located at a distance of $r$ from the center of mass
(c.m.) of a specified reference particle. The RDF is of great significance in
condensed matter physics as it can directly be related to the static structure
factor, $S(k)$, and hence, determined experimentally from radiation scattering
experiments, such as those using x rays and neutrons [37]. For an infinite
isotropic and homogeneous system, the relation between the RDF and $S(k)$ is
given by [38, 39, 40]

$\displaystyle
S(k)=1+4\pi\rho\int_{0}^{\infty}r^{2}[g(r)-1]\frac{\sin{kr}}{kr}dr,$ (1)

where $\rho$ is the system averaged number density of atoms. Physically,
$S(k)$ describes the system density response at wavelength $2\pi/k$ to a weak
enough, external perturbation [41]. Note that $S(k)$ is always nonnegative in
equilibrium.

The RDF is also used to link microscopic structural details to macroscopic
properties under the Kirkwood-Buff (KB) solution theory [42]. In the canonical
ensemble, for example, the potential of mean force ($w$) between any pair of
particles in the fluid is related to the RDF via [43]

$\displaystyle w_{N}(r)=-k_{B}T\ln{g_{N}(r)},$ (2)

where the subscript $N$ is added to highlight that the total number of
particles is constant, and $T$ is the equilibrium temperature. The $w$ is
conveniently written as a sum of two terms if the total potential energy is
approximated by a sum of identical, independent pair potential energies,
$u(r)$:

$\displaystyle w_{N}(r)=u(r)+\delta F_{N}(r),$ (3)

where $\delta F_{N}(r)$ includes the effects of the solvent and is the
canonical-ensemble average change of $F_{N}$, the Helmholtz free energy, of
the fluid introduced by bringing two atoms from infinity to a (finite)
distance $r$. If the process is done adiabatically, $\delta F_{N}(r)$ can be
considered as the work done on the solvent, that is the remaining $N-2$ atoms,
during the process. Generally speaking, $\delta F_{N}$ is nonzero at finite
number densities. Using finite-size corrections [44, 45, 46, 47, 48, 49], it
can be shown for $N>>1$ and $r$ greater than the correlation length $\xi$ at
any density, or all $r$ at $\rho\to 0^{+}$,

$\displaystyle\delta F_{N}(r)=\frac{k_{B}T\chi^{\infty}_{T}}{N},$ (4)

where $\chi_{T}^{\infty}=S(0)$ is the reduced isothermal compressibility of an
open system [in the thermodynamic limit (TL), that is $N,L\rightarrow\infty$,
while $\rho$ is constant], which shares the same equilibrium state [47].

There is a great number of research articles in the literature attempting to
determine short- or long-range behaviours of the RDF in order to understand a
wider range of phenomena, such as the wetting phenomenon [50, 51]. It has
rigorously been shown for a system in equilibrium with an interparticle
potential which either decays faster than a power law or is truncated at a
finite cutoff radius, the RDF for $r>r_{c}$, where $0<r_{c}<\xi$ is the cutoff
radius or the (effective) range of the pair potential energy, $u(r)$, is given
by (in the TL) [52, 53, 54, 55]

$\displaystyle g(r)=1+\left(\frac{R_{\xi}\cos{\beta_{\xi}r}}{2\pi
r}\right)e^{-r/\xi},$ (5)

where $R_{\xi}$ is the residue of the Fourier transform of $[g(r)-1]$
corresponding to the pole(s) with the smallest positive imaginary part, and
$\beta_{\xi}$ is some (real) constant. Equation (5) shows $r[g(r)-1]$ in
equilibrium decays asymptotically to zero exponentially, either monotonically
(with $\beta_{\xi}=0$) or sinusoidally (with $\beta_{\xi}\neq 0$) [56, 57,
58].

In the canonical ensemble, the average entropy $\langle S_{N}\rangle$ of an
atomic fluid at a temperature $T$ with $N$ indistinguishable atoms, described
by the canonical space and momentum coordinates
$(\textbf{r}_{1},\textbf{p}_{1}),\dots,(\textbf{r}_{N},\textbf{p}_{N})$, can
be calculated from the generalized Gibbs-Shannon formula [59, 60]:

$\displaystyle\begin{split}\langle S_{N}\rangle=-\frac{k_{B}}{N!}\int
d\textbf{r}_{1}d\textbf{p}_{1}\dots
d\textbf{r}_{N}d\textbf{p}_{N}\mathcal{P}_{N}\ln{[h^{3N}\mathcal{P}_{N}]},\end{split}$
(6)

where $\mathcal{P}_{N}$ is the (unnormalized) phase-space probability density,
$k_{B}$ is the Boltzmann’s constant, $h$ is the Planck’s constant, and $N!$
accounts for indistinguishablity of particles. If $H_{N}$, the total
Hamiltonian, separates into additive terms for the potential and kinetic
energy, the phase-space probability density can be factorized as

$\displaystyle\mathcal{P}_{N}=g_{N}(\textbf{r}_{1},\dots,\textbf{r}_{N})\prod_{i=1}^{N}f_{1}(\textbf{p}_{i}),$
(7)

where $f_{1}(\textbf{p}_{i})$ is the $1$-particle momentum probability
density:

$\displaystyle f_{1}(\textbf{p}_{i})=\rho(2\pi
mk_{B}T)^{-3/2}e^{-|\textbf{p}_{i}|^{2}/2mk_{B}T}$ (8)

Equation (7) serves as a definition for the $N$-particle positional
distribution function, $g_{N}(\textbf{r}_{1},\dots,\textbf{r}_{N})$;
Physically, it is a measure of the joint probability of finding the particle
$1$ at position $\textbf{r}_{1},\dots$, and particle $N$ at position
$\textbf{r}_{N}$. Utilizing the generalized Kirkwood superposition
approximation, one can factorize $g_{N}(\textbf{r}_{1},\dots,\textbf{r}_{N})$
as [61, 62]

$\displaystyle\begin{split}g_{N}(\textbf{r}_{1},\dots,\textbf{r}_{N})=g_{N}(\textbf{r}_{1},\textbf{r}_{2})\times\dots\times
g_{N}(\textbf{r}_{N-1},\textbf{r}_{N})\\\ \times\delta
g_{N}(\textbf{r}_{1},\textbf{r}_{2},\textbf{r}_{3})\times\dots\times\delta
g_{N}(\textbf{r}_{1},\dots,\textbf{r}_{N}),\end{split}$ (9)

where, for example,

$\displaystyle\delta
g_{N}(\textbf{r}_{1},\textbf{r}_{2}.\textbf{r}_{3})\equiv\frac{g_{N}(\textbf{r}_{1},\textbf{r}_{2},\textbf{r}_{3})}{g_{N}(\textbf{r}_{1},\textbf{r}_{2})g_{N}(\textbf{r}_{1},\textbf{r}_{2})g_{N}(\textbf{r}_{2},\textbf{r}_{3})}.$

By substituting Eq.(7) into Eq.(6) and using Eqs.(8) and (9), an entropy
expansion is obtained, which is well approximated by [63, 64, 65, 66, 67, 68]

$\displaystyle\langle S_{N}\rangle=\langle S_{N}^{id}\rangle+S_{N}^{(2)}$ (10)

for high (near the freezing point) and near-zero densities, where the first
term on the r.h.s is the ideal-gas mean entropy, and the second term is the
$2$-body (excess) entropy.

For a homogeneous and isotropic fluid, one can express the 2-body (excess)
entropy in terms of the RDF [69] and $w$ [by using Eq.(2)]:

$\displaystyle S_{N}^{(2)}=\frac{N\rho}{2T}\int
d\mathbf{r}g_{N}(r)w_{N}(r)-\frac{Nk_{B}}{2}$ (11)

Substituting $w_{N}(r)$ in Eq.(11) with the expression in Eq.(3), the
following expression is obtained:

$\displaystyle S_{N}^{(2)}=\frac{\langle U_{N}\rangle}{T}+\delta
S_{N}^{(2)}-\frac{Nk_{B}}{2},$ (12)

where $\langle U_{N}\rangle$ is the system’s total potential energy,
$Nk_{B}/2$ is the entropy contribution of the canonical ensemble, and

$\displaystyle\delta S_{N}^{(2)}=\frac{N\rho}{2T}\int
d\mathbf{r}g_{N}(r)\delta F_{N}(r),$ (13)

which is the solvent entropy contribution and referred to as the residual
2-body (excess) entropy. For low densities, one can show for $N>>1$,

$\displaystyle\lim_{\>\rho\>\to\>0^{+}}\delta
S_{N}^{(2)}=\frac{Nk_{B}\chi_{T}^{\infty}}{2}$ (14)

by direct substitution of Eq.(4) into Eq.(13). Utilizing Eqs.(12) and (14),
the near-zero density limit of Eq.(10) can be written for $N>>1$ as

$\displaystyle\lim_{\>\rho\>\to\>0^{+}}\langle S_{N}\rangle=\langle
S_{N}^{id}\rangle+\frac{\langle
U_{N}\rangle}{T}+\frac{Nk_{B}}{2}(\chi_{T}^{\infty}-1).$ (15)

In the ideal-gas limit, $\langle U_{N}\rangle=0$ and $\chi_{T}^{\infty}=1$.

##  III  Theory

A closed system is in thermodynamic equilibrium, provided that its Helmholtz
free energy is globally minimized [70, 71]. The Helmholtz free energy for a
system in contact with a heat bath is mathematically defined as [72, 73]

$\displaystyle F_{N}:=\langle E_{N}\rangle-T\langle S_{N}\rangle,$ (16)

where $\langle E_{N}\rangle=\langle E^{id}_{N}\rangle+\langle U_{N}\rangle$ is
the sum of the kinetic and potential energy, and -$T\langle S_{N}\rangle$ is
the mean heat exchanged between the system and its environment. Replacing
$\langle E_{N}\rangle$ and $T\langle S_{N}\rangle$ with their equivalent
expressions in Eq.(16) and using Eq.(12), the following expression for a
closed $NVT$ system at high (near-freezing) or near-zero densities is
obtained:

$\displaystyle F_{N}=F_{N}^{id}(T)+T\left[\frac{Nk_{B}}{2}-\delta
S_{N}^{(2)}\right],$ (17)

where $F_{N}^{id}(T)$ is the ideal-gas Helmholtz free energy.

The system relaxation time, denoted $\tau_{sys}$, is defined as the inverse
rate of approach to (thermodynamic) equilibrium. Based on Eq.(17), the
Helmholtz free energy equilibrates once the temperature $T$ and residual two-
body excess entropy reach equilibrium. Thus, for such systems

$\displaystyle\tau_{sys}=max(\tau_{T},\tau_{\delta S}),$ (18)

where $\tau_{T}$ and $\tau_{\delta S}$ are the relaxation times for the
temperature $T$ and residual two-body excess entropy, respectively. In
statistical mechanics, the (kinetic) temperature $T$ can be defined via the
equipartition theorem [74] as the mean kinetic energy per atom ($T:=2\langle
E^{id}_{N}\rangle/3Nk_{B}$). It is, therefore, a local quantity which should,
in general, equilibrate faster than the residual two-body entropy, whose
equilibration depends on the equilibration of all the particles’ positions
within the system, (at least for local and particle-wise thermostating). Thus,
we typically always have $\tau_{\delta S}>\tau_{T}$ and

$\displaystyle\tau_{sys}=\tau_{\delta S}.$ (19)

In nonequilibrium thermodynamics [75], the definition of entropy, i.e.,
Eq.(6), is also used for systems near equilibrium. However, this assumption
may not be valid for systems in which the local equilibration hypothesis is
never satisfied, such as in glasses [76, 77]. For nonglassy (atomic) systems
in contact with a local thermostat, such as the Langevin or dissipative
particle dynamics (DPD) thermostat, one would expect the temperature $T$ is
equilibrated before local equilibration is reached (for systems with local
thermostats, in local equilibrium, temperature is in global equilibrium,
whereas structure is in local equilibrium). Hence, for such systems, Eq.(13)
is generalized for $t>\tau_{T}$, ensuring the systems have most likely reached
local equilibrium, as

$\displaystyle\delta S_{N}^{(2)}(t)=\frac{N\rho}{2T}\int
d\mathbf{r}g_{N}(r,t)\delta F_{N}(r,t),$ (20)

with

$\displaystyle\begin{split}g_{N}(r,t)&=e^{-w_{N}(r,t)/k_{B}T}\\\
&=e^{-u(r)/k_{B}T}e^{-\delta F_{N}(r,t)/k_{B}T},\end{split}$ (21)

where the second line comes from generalizing Eq.(3). Equation (20) implies
that $\tau_{\delta S}$ is equal to $\tau_{\delta F}$, the relaxation time for
$\delta F_{N}(r,t)$. Thus, for nonglassy atomic systems at near-freezing or
zero densities, $\tau_{sys}=\tau_{\delta F}$ [refer to Eq.(19)].

In generic physical systems, it may not be always feasible to model the
relaxation processes of dynamical quantities with an exponentially decaying
function, especially when those systems undergo at least one phase transition.
Nonetheless, one would expect that most of the effects of phase transitions on
relaxation processes disappear once local equilibrium is reached. If $\delta
F_{N}(r,t)$ in Eq.(21) is substituted with (the exponential ergodicity
hypothesis [78])

$\displaystyle\delta F_{N}(r,t)=\delta F_{N}(r)+\delta
F^{\prime}_{N}(r)e^{-t/\tau_{\delta F}}$ (22)

(for $t>\tau_{T}$), where the second term on the RHS indicates the extra
(positive or negative) adiabatic work needed in nonequilibrium conditions to
bring the two particles from infinity to a distance $r$, the instantaneous RDF
for a nonglassy system (at any density) is written as follows:

$\lim_{\>t\>>\>\tau_{T}}g_{N}(r,t)=g_{N}(r)e^{\alpha_{N}(r)e^{-t/\tau_{\delta
F}}}\\\ =g_{N}(r)\\{1+\alpha_{N}(r)e^{-t/\tau_{\delta
F}}+O(e^{-2t/\tau_{\delta F}})\\},$ (23)

where

$\displaystyle\alpha_{N}(r):=-\frac{\delta F_{N}^{\prime}(r)}{k_{B}T}$

is a signed dimensionless physical quantity including nonequilibrium
(structural) information. In the above Taylor series, the higher-order terms
are neglected as they decay much faster than the first two terms in the
brackets. This paper is only concerned with situations in which $N>>1$.
Hereafter, we drop the subscript $N$ from all the quantities to stress finite-
size effects are negligible for $N>>1$.

In computer simulations, if the initial structure is more expanded compared to
the final structure, one would expect $\delta F^{\prime}$ to be, on average,
positive (equivalent to $\alpha<0$); This is because the initial number
density distribution is more uniform than it should be. Consequently, less
spots are, on average, available to bring two atoms from infinity to a
distance $r$ without disrupting the structural stability. Hence, one needs
some extra positive work to make space for them. On the other hand, if the
initial structure is more compacted than the final structure, one would expect
$\delta F^{\prime}$ to be, on average, negative (equivalent to $\alpha>0$);
This is because there will be more marginal space initially available (in the
simulation box) than in the final state, which makes moving two atoms from
infinity to a finite distance without changing the position of other atoms
more likely, leading to averagely less total (adiabatic) work or $\delta
F^{\prime}<0$ [refer to Eq.(22); $\delta F(r,t)<\delta F(r)\Rightarrow\delta
F^{\prime}<0$].

Equation (23) can be further simplified to (at any density)

$\lim_{\>t\>>\>\tau_{T}}\lim_{\>r\>>\>\xi}g(r,t)=1+\alpha e^{-t/\tau_{\delta
F}}$ (24)

since for $r>\xi$, the RDF $g(r)\approx 1$ and $\alpha(r)$ is expected to be
(effectively) $r$ independent, provided the structural effects of the system-
thermostat coupling are small enough in order not to significantly perturb the
system’s natural dynamics. At near-zero densities, $\xi\approx r_{c}$ and
$\tau_{\delta F}\approx\tau_{sys}$. Substituting these values for near-zero
densities in Eq.(24) gives

$\displaystyle\lim_{\>t\>>\>\tau_{T}}\lim_{\>r\>>\>r_{c}}g(r,t)=1+\alpha
e^{-t/\tau_{sys}}.$ (25)

The sign of $\alpha$ determines how the RDF (during equilibration) decays to
one in time. Equation (25) is of particular significance since it can be used
to study (possible) effects of initial structural inhomogeneity on $\alpha$
and the time for the system relaxation process $\tau_{sys}$ in the Langevin
model. In principle, Eq.(25) can also be used at high densities (near the
freezing point) for $r>\xi\>(\neq r_{c})$. In practice, $\tau_{sys}$ decreases
with increasing number density $\rho$ at a constant temperature [79]. Hence,
it is more difficult to observe the (long-range) evolution of the RDF in time
at high densities compared to near-zero ones. We mainly focus on dilute
systems in this paper.

##  IV  Methodology and Simulation Details

This article deals primarily with how the initial structure can affect the
system relaxation time in the Langevin model. For this purpose, we performed
many molecular-dynamics (MD) simulations utilizing the LAMMPS software package
[80] and computed the time-averaged RDF, $g(r)$, after the system’s
temperature $T$ is relaxed for three different well-known and well-used
initial structures. The pair potential energy $u(r)$ is a truncated and
shifted 6-12 Lennard-Jones (LJ) potential:

$\displaystyle u(r)=\begin{cases}V_{LJ}(r)-V_{LJ}(r_{c}),&\text{for $0<r\leq
r_{c}$}\\\ 0,&\text{for $r>r_{c}$}\end{cases},$ (26)

where $r_{c}=2.5\sigma$ is the cut-off radius, and

$\displaystyle
V_{LJ}(r)=4\varepsilon\>\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],$

with $\sigma$ and $\varepsilon$ being the effective atomic radius and the
dispersion energy (depth of the potential well), respectively. The 6-12-LJ
equation approximates well the interactions in the noble massive gases (e.g.,
Ar and Kr), whose interactions are dominated by van der Waals forces. Note
that the constant term $V_{LJ}(r_{c})$ in the above equation is added to avoid
any discontinuity at $r=r_{c}$ in the potential, leading to an impulsive
contribution to the (thermodynamic) pressure [81].

In a computer experiment, it is often convenient to express physical
quantities in units other than the SI units. For example, by expressing these
quantities in LJ reduced units allows one to benefit from the corresponding-
states principle [82]. For the LJ reduced units,
$k_{B}=m=\sigma=\varepsilon=1$. Also, using LJ reduced units reduces the
floating-point/round-off error (caused by the limited precision of
computational processors) as parameters are of order one [81]. We note that
all the physical quantities quoted below are in LJ units.

A pure monatomic fluid at number density $\rho_{0}=3.95\times
10^{-2}=\rho_{cr}^{Ar}/8<<1$ ($\rho_{cr}^{Ar}$ is the critical density of
Argon) was simulated in the Langevin canonical ensemble with zero total linear
momentum and with periodic boundaries to remove surface effects. In Refs. [83,
84], it is shown that in MD experiments with zero total linear momentum and
periodic boundaries, the infinitesimal Galilean boost is conserved, which
ensures each atom remains an inertia reference frame. In the Langevin model,
the equation of motion for the $i$-th particle relative to the system center
of mass reference frame ($R_{c.m.}=0$) is given by [85]

$\ddot{\mathbf{r}}_{i}(t)=\mathbf{F}_{i}+\mathbf{g}_{i}(t)\quad(i=1,2,\dots,N),$
(27)

where $\mathbf{F}_{i}$ is the net conservative force on the $i$-th particle,
and $\textbf{g}_{i}(t)$ is a fictitious force, which modifies the dynamics of
the system to account for the presence of the Langevin thermostat.
$\textbf{g}_{i}(t)$ is responsible for the variations leading to the system
thermalization and is given by [85]

$\displaystyle\textbf{g}_{i}(t)dt=-\gamma\dot{\textbf{r}}_{i}(t)dt+\sqrt{2\gamma
T}d\textbf{W}_{i}(t),$ (28)

where $\gamma>0$ is the Langevin coupling constant, and $d\textbf{W}_{i}$ is a
vector of independent Wiener processes, satisfying

$\begin{gathered}\langle
dW_{ij}(t)dW_{kl}(t^{\prime})\rangle=\delta_{ik}\delta_{jl}\delta(t-t^{\prime})dt\\\
\langle dW_{ij}(t)\rangle=0\\\ (j,l=x,y,z)\text{ and
}(i,k=1,\dots,N),\end{gathered}$ (29)

where $\delta(t-t^{\prime})$ is the Dirac delta distribution, $\delta_{ik/jl}$
is the Kronecker delta function, and $\langle\dots\rangle$ indicates an
average over an ensemble of simultaneous, independent, and similar
experiments.

In our computer experiments, the equations of motion (27) are integrated
numerically utilizing the GJF-2GJ algorithm, which is quite accurate,
especially for the kinetic sampling of the phase space [86]. It has also been
proven that the GJF-2GJ algorithm is capable of providing exact thermodynamic
responses for constant and harmonic potentials for any timestep size $\delta
t$ within the verlet stability criteria. In the GJF-2GJ algorithm, atom’s
velocities are computed at middle points of time steps. The $i$-th atom’s
half-step velocity relative to the system’s c.m. ($R_{c.m.}=0$) is then given
by [87]

$\displaystyle\textbf{v}_{i}(t+\frac{\delta
t}{2})=\frac{\textbf{r}_{i}(t+\delta t)-\textbf{r}_{i}(t)}{\sqrt{b}\delta
t}\quad(i=1,2,\dots,N),$ (30)

where $\textbf{r}_{i}$ is the on-site atomic position, $\delta t$ is the
timestep size, and

$\displaystyle b=(1+\frac{\gamma\delta t}{2})^{-1}$ (31)

to preserve semi-symplecticity. In computer experiments, the resulting
equilibrium distributions are different from the true theoretical
distributions as a result of discretization or truncation errors [88, 89]. To
avoid losing physical information due to such systematic errors, the timestep
should be much smaller than the inverse of the fastest vibrational frequency
in the system. In this paper, we chose $\delta t=0.0007$. Note that in
numerical calculations with computers, the roundoff error $\propto 1/\delta t$
[90].

We first created $N_{1}\sim 10^{5}$ LJ particles at temperature $T=0$ in the
form of a face-centered cubic crystal, consisting of $43^{3}$ cells with
lattice constant $a=(4/\rho_{0})^{1/3}=4.66$ and locating at the center of a
simulation box with a volume of $V=(45a)^{3}$. This structure is slightly
marginally inhomogeneous in the number density of atoms [see Fig. 1(a)] due to
the empty gap around the system’s edges. Then, we created $N_{2}\sim 10^{5}$
LJ particles under the same conditions, except this structure now fills the
entire simulation box. Such a structure is macroscopically homogeneous in the
number density of particles [see Fig. 1(b)]. Both of these initial structures
are commonly used in MD simulations. The second for the reason that
macroscopically homogeneous densities are expected to be desirable initial
conditions, and the first as a common approximation to the second without
having to worry about lattice commensurability effects at the periodic
boundaries. For Fig. 1(b), LAMMPS is careful to put only one particle at the
boundaries to avoid any unwanted atom overlap (using the box style of the
create_atoms command in LAMMPS).

(a) marginally inhomogeneous structure (FCC crystal)

(b) uniform structure (FCC crystal)

(c) uniform structure (amorphous solid)

Figure 1: 2D illustration of the initial structures.

Unless otherwise stated, for each case, the phase-space trajectories were
first run for $t_{e}=\Delta t=175$ to equilibrate the system’s temperature at
a supercritical temperature $T=3.936$ for a few multiples of a basic Langevin
thermal coupling constant of $\gamma=\gamma_{0}=5/7$. As will be confirmed in
the next section, $t_{e}=175\>(>>\tau_{T})$ is more than enough time to
equilibrate the system’s temperature for all the $\gamma$ values used. The
trajectories were then run for a duration of $t_{pro}=\Delta t=525$ so as to
calculate the (time-averaged) RDF, $g(r)$, up to a maximum distance of
$r_{m}=10r_{c}=25$. Finally, to ensure that the results are not limited to
only those systems with crystalline initial structures, we also performed a
number of extra MD simulations with amorphous initial structures (randomly-
distributed particles), which occupied the entire simulation box [see Fig.
1(c)]. This is the only distinction between the second and third batches of
simulations. It should also be noted that this kind of initial structure is
typically (energetically) unstable because randomly-generated particles are
often highly overlapped. Therefore, we had to perform a local pre-energy
minimization using the Conjugate-Gradient (CG) algorithm for $\Delta t\simeq
2.1$ for this initial state.

###  IV.1 Measurement of the time-averaged RDF

In an MD experiment, the time-averaged RDF is measured as follows: First, a
random atom is selected. Next, the algorithm computes $\delta n(r,t)$. Then,
the algorithm increases the distance from $r$ to $r+\delta r$ and does exactly
the same measurements [91]. This procedure continues until $r$ reaches its
maximum, that is $r_{m}=10\times r_{c}=25$. The algorithm repeats the whole
process for another randomly-chosen atom, and so on. In the end, the algorithm
determines the instantaneous RDF, defined as

$\displaystyle g(r,t):=\frac{1}{4\pi N\rho_{0}r^{2}}\sum_{i=1}^{N}\frac{\delta
n_{i}(r,t)}{\delta r},$ (32)

where the summation is over particles and $\delta n_{i}(r,t)$ is the number of
particles within the spherical shell of thickness $\delta r=25/256$ centered
at the $i$-th atom between time $t$ and $t+\delta t$. The time-averaged RDF is
computed by averaging $g(r,t)$ over time:

$\displaystyle g(r)=\frac{\delta t}{t_{pro}}\sum_{j=1}^{t_{pro}/\delta
t}{g(r,t_{e}+j\delta t)}$ (33)

with $t_{e}=175$ and $t_{pro}=525$ being the equilibration and production
time, respectively. In the previous section, we derived an (approximate)
expression for $g(r,t)$ for $r>r_{c}$ and $t>\tau_{T}$. By substituting
Eq.(24) into Eq.(33), one finds that if the system is not well equilibrated,
the time-averaged RDF undergoes a shift. For high (near the freezing point)
and near-zero densities,

$\displaystyle\lim_{\>t_{e}\>>\>\tau_{T}}\lim_{\>r\>>\>\xi}[g(r)-1]=Ae^{-t_{e}/\tau_{sys}},$
(34)

where $\xi\approx r_{c}$ at low densities, and

$A=-\frac{\delta
F^{\prime}}{k_{B}T}\frac{\tau_{sys}}{t_{pro}}\left(1-e^{-t_{pro}/\tau_{sys}}\right)\propto-\delta
F^{\prime}.$

In Langevin dynamics, Eq.(34) may be generalized as follows to incorporate
effects of the system-thermostat coupling (weak coupling):

$\displaystyle\lim_{\>t_{e}\>>\>\tau_{T}(\gamma)}\lim_{\>r\>>\>\xi}[g(r)-1]=A(\gamma)e^{-t_{e}/\tau_{sys}(\gamma)}$
(35)

In this paper, we study Eq.(35) for some (common) values of $\gamma\;(\sim 1)$
to see how the RHS may respond to the three well-used initial structures in a
molecular-dynamics simulation.

##  V  Results and Discussion

###  V.1 Temperature relaxation time

In Fig. 2(a), we plotted the evolution of instantaneous temperature $T$ for
the homogeneous crystalline structure shown in Fig. 1(b) for some multiples of
a Langevin coupling constant, $\gamma_{0}=5/7=1/2000\delta t$, to demonstrate
the system is in thermal equilibrium for all $t\geq 175$. This figure shows
the stronger the temperature coupling between the system and the Langevin
thermostat, the faster the system reaches thermal equilibrium. As expected,
for an exponentially-ergodic system, the time evolution of quantities, such as
temperature $T$, can be approximated and well fitted by an exponential curve
[22]. The relaxation time for quantity $X$, denoted by $\tau_{X}$, is then the
(characteristic) time constant of the fitted exponential function. We have
computed the temperature relaxation time for some $\gamma$ values and plotted
them versus $(\gamma)^{-1}$ in the inset. As is expected, the temperature
relaxation time $\tau_{T}$ varies linearly with the inverse of the thermostat
coupling constant $\gamma$ as $\tau_{T}\approx(\gamma)^{-1}/2$. In Langevin
dynamics, by using the Ito’s Lemma Chain Rule [92], we obtain the following
expression for the system’s kinetic temperature $T$ during equilibration [85]:

$\displaystyle\langle T(t)\rangle=(T_{i}-T_{f})e^{-t/\tau_{T}}+T_{f},$ (36)

where $\tau_{T}=1/2\gamma$ is the temperature relaxation time, and $T_{f}$
($T_{i}$) is the final (initial) equilibrium temperature. For this paper,
$T_{i}=0$ and $T_{f}=3.936=3\times T_{cr}^{Ar}$ ($T_{cr}^{Ar}$ is the critical
temperature of argon).

(a) uniform initial structure

(b) marginally inhomogeneous initial structure

Figure 2: Temporal evolution of the instantaneous (kinetic) temperature $T$
during the equilibration phase for both kinds of initial structures and some
multiples of $\gamma_{0}=5/7$. As can be observed, the higher the value of the
Langevin coupling constant $\gamma$, the faster the system reaches thermal
equilibrium. Both insets depict that the (kinetic) temperature relaxation
time, $\tau_{T}$, varies linearly with $1/\gamma$. To be precise,
$\tau_{T}\approx 1/2\gamma$. (Note that all the values are in the LJ units.)

In Fig. 2(b), similar to the homogeneous structure, we plotted the temporal
evolution of temperature during the system’s temperature equilibration for the
same $\gamma$ values, but now for the inhomogeneous structure shown in Fig.
1(a). Then, we have calculated the temperature relaxation time for each
$\gamma$ and plotted them versus $1/\gamma$ in the inset. As expected, the
temperature relaxation time varies with the Langevin coupling constant
$\gamma$ as $\tau_{T}\approx 1/2\gamma$, which is in excellent agreement with
Eq.(36). The simulations with the amorphous initial structure [see Fig. 1(c)]
also show the same outcomes (not displayed here). From the above figures, it
is clear that $t_{e}=175>>\tau_{T}$ and the temperature relaxation time is
independent of the initial structure in Langevin dynamics.

(a) uniform initial structure

(b) marginally inhomogeneous initial structure

Figure 3: $r[g(r)-1]$ for both kinds of initial structures, $t_{e}=175$, and
some multiples of $\gamma_{0}=5/7$. These figures reveal that, although for
the uniform initial structure, $r[g(r)-1]$ is independent of $\gamma$ at
thermal equilibrium, it is, indeed, highly $\gamma$ dependent for the
marginally inhomogeneous initial structure. (Note that the insets are a
zoomed-in version of the figures and all the values are in the LJ units.)

###  V.2 Examining Eq.(35)

For dilute systems and $r>r_{c}$ we expect $[g(r)-1]$ to be infinitesimally
small if the system is fully equilibrated. To make any deviation from the
equilibrium more easily visible, in this section, we plot $r[g(r)-1]$, instead
of $[g(r)-1]$ (i.e., we scale it by $r$). For Fig. 3(a), we plotted
$r[g(r)-1]$ for the case of the uniform structure of Fig. 1(b), the
equilibration time $t_{e}=175>>\tau_{T}$ and some multiples of
$\gamma=\gamma_{0}\approx 0.7$. This figure demonstrates that the RDF is
independent of $\gamma$ for $t_{e}=175$. Thus, $t_{e}>>\tau_{sys}$ for
$\gamma=\gamma_{0},\>2\gamma_{0},\>4\gamma_{0}$ [based on Eq.(35)]. At a very
high value of $\gamma=70\gamma_{0}\>(=50)$, we found that the RDF is slightly
out of equilibrium (not shown) which means the system relaxation time
$\tau_{sys}$ has increased for $\gamma=70\gamma_{0}$ such that $t_{e}$ is not
much longer than $\tau_{sys}$ anymore. Afterwards, we repeated the simulations
for the case of the randomly-dispersed structure of Fig. 1(c). As is expected,
we found exactly the same behaviour for $r[g(r)-1]$ (not shown here). This
would suggest that the system relaxation times of initial structures which
are, on average, uniformly distributed throughout the simulation box are of
the same order and only weakly dependent on the Langevin coupling constant
(increasing behaviour).

(a) nonuniform initial structure

(b) nonuniform initial structure

Figure 4: (a) Plot of $r[g(r)-1]$ versus $r$ for $\gamma=2\gamma_{0}=10/7$ and
some multiples of $t_{e}=175$. As seen, the long-range slope decreases
continuously with increasing the equilibration time. (Note that the inset is a
zoomed-in version of the figure and all the values are in the LJ units.) (b)
The slope of $\lim_{\>r\>>\>r_{c}}r[g(r)-1]$ vs. the equilibration time
$t_{e}$ for $\gamma=\gamma_{0}=5/7$ and $2\gamma_{0}=10/7$. The dotted and
dashed lines are the exponential fits to the data from our MD simulations. The
time constant of each exponential fit is the system relaxation time [based on
Eq.(35)]. Figure 5: How and in which direction the RDF evolves over time for
two different initial structures at instant $t=2.8$. The RDF graph for the
inhomogeneous case, unlike the homogeneous one, approaches one from above,
which is the consequence of initial structural inhomogeneity. (Note that all
the values are in the LJ units.)

In Fig. 3(b), we plotted $r[g(r)-1]$ for the initial structure of Fig. 1(a)
and the equilibration time $t_{e}=175\;(>>\tau_{T}$; the system is assumed to
have already reached local thermal equilibrium), and some multiples of
$\gamma_{0}=5/7$. As is expected from Eq.(35), $\lim_{\>r\>>\>r_{c}}r[g(r)-1]$
is $\gamma$ dependent in local equilibrium. In particular,
$\lim_{\>r\>>\>r_{c}}r[g(r)-1]$ varies linearly with a $\gamma$ dependent
slope for weak thermal couplings and nonlinearly for strong thermal couplings.
Such deviation from zero is because (structural) entropy is not maximal and
$\tau_{T}<<t_{e}<\tau_{sys}$. By comparing these two figures, we can conclude
(hom: homogeneous, inhom: inhomogeneous)

$\tau_{sys}^{inhom}>>\tau_{sys}^{hom}.$ (37)

Therefore, an initial small marginal inhomogeneity in the Langevin model
significantly increases the time required to equilibrate the system’s
structure. Further, increasing the Langevin constant (which reduces the
thermal equilibration time) makes the equilibration significantly worse. The
inhomogeneity in the marginal inhomogeneous initial state is primarily of a
long-wavelength nature and the Langevin thermostat suppresses or slows down
the evolution of the periodic and long-wavelength density fluctuations needed
to equilibrate these modes.

In Fig. 4(a), we plotted $r[g(r)-1]$ for the inhomogeneous structure for
$\gamma=2\gamma_{0}=10/7$ and several multiples of $t_{e}=175$. This figure
clearly reveals that the slope decreases with increasing the duration of the
equilibration phase, i.e., $t_{e}$. Hence, the non-zero slope is basically
because the amount of the equilibration time is not enough to equilibrate the
system’s structure. Next, in Fig. 4(b), we plotted the evolution of the slope
with the equilibration time for $\gamma=\gamma_{0},\;2\gamma_{0}$. We found
that the slope varies with $t_{e}$ in an exponential manner in local
equilibrium, as exactly predicted in Eq.(35). In Fig. 4(b),

$\text{Slope}=\frac{\int_{r_{c}=2.5}^{r_{m}=25}[g(r)-1]dr}{r_{m}-r_{c}}=A(\gamma)e^{-t_{e}/\tau_{sys}(\gamma)},$
(38)

where $A(\gamma)$ was found to be positive. The numerical results obtained
from Fig. 4(b) for $\tau_{sys}(\gamma)$, i.e.,
$\tau_{sys}(\gamma_{0})=(0.013)^{-1}$ and
$\tau_{sys}(2\gamma_{0})=(0.006)^{-1}\approx 2\times\tau_{sys}(\gamma_{0})$,
suggest an increasing behaviour for the system relaxation time with increasing
the Langevin constant $\gamma$. This is in contrast to the (kinetic)
temperature relaxation time $\tau_{T}$ which we saw decrease with increasing
$\gamma$ [$\tau_{sys}(\gamma)\propto\gamma$ vs. $\tau_{T}(\gamma)\propto
1/\gamma$]. This increasing behaviour was, indeed, predictable due to the
structural disruptive effect caused by the Langevin thermostat.

In Fig. 5, we plotted the RDF for $\gamma=\gamma_{0}$ at instant
$t=2.8=4\tau_{T}(\gamma_{0})$ for both homogeneous [see Fig. 1(b)] and
inhomogeneous [see Fig. 1(a)] initial structures. As is seen, in the
marginally inhomogeneous case, unlike the uniform one, $g(r)$ approaches one
in time from above over the system’s relaxation. To be specific, based on
Eq.(35), for the inhomogeneous structure (i.e., the compacted structure),
$A>0$, and for the homogeneous case (the more expanded structure), $A<0$. This
sign change was, indeed, predictable since $A\propto-\delta F^{\prime}$, the
extra (positive/negative) adiabatic work to move two atoms during the
structural equilibration from infinity to a distance $r$. In the Theory
section, we showed that the sign of $\delta F^{\prime}$ depends on the initial
structure, and in the inhomogeneous case, unlike the uniform one, it is
negative ($\delta F_{inhom}^{\prime}<0$). This sign change can be considered
as the second effect of initial structural inhomogeneity on the system
structural relaxation. [$g(r,t)$ graphs evolve in the directions shown in Fig.
5 before the system reaches equilibrium.]

Finally, to confirm that the low-density results are also applicable to higher
densities, we have plotted $r[g(r)-1]$ in Fig. 6 at density
$\rho=16\times\rho_{0}$ for homogeneous [Fig. 1(b)] and inhomogeneous [Fig.
1(a)] initial structures, $\gamma=\gamma_{0},\>2\gamma_{0},\>4\gamma_{0}$,
$t_{e}=7\geq 10\tau_{T}$, and $t_{pro}=4.9\geq 7\tau_{T}$. As expected, the
homogeneous graphs, unlike the inhomogeneous ones, are all fully equilibrated
for the same values of $t_{e}$ and $t_{pro}$. Hence, one would conclude,
assuming exponential ergodicity, $\tau_{sys}^{inhom}>\tau_{sys}^{hom}$ (and
$\uparrow\tau_{sys}\sim\;\uparrow\gamma$) also holds at higher densities. An
interesting point about this figure is that $\lim_{r>\xi}A$ for
$\gamma=4\gamma_{0}=2.86$ (similar to the case of $\gamma=70\gamma_{0}=50$ at
$\rho=\rho_{0}$) is $r$ dependent for $r>\xi$; this is because this system is
now in the Brownian dynamics regime, that is strong thermal couplings. These
observations imply that at higher densities, the minimum $\gamma$ value
required for Brownian dynamics is smaller than that for low densities. Figure
6 demonstrates the validity of our theory and its predictions at high
densities.

Figure 6: $r[g(r)-1]$ at density $\rho=16\times\rho_{0}=0.632$ for both kinds
of initial structures, $t_{e}=7$, $t_{pro}=4.9$, and some multiples of
$\gamma_{0}=5/7$. This figure demonstrates the RDF in all the homogeneous
cases is fully equilibrated regardless of the value of $\gamma$. However, for
the marginally inhomogeneous cases, $r[g(r)-1]$ is $\gamma$ dependent,
indicating the RDF has not been well equilibrated (Note that the inset is a
zoomed-in version of the figure, and all the values are in the LJ units.)

##  VI  Conclusion

In this paper, we found that the system relaxation time is not always
independent of the initial state in molecular-dynamics simulations. In
particular, in Langevin dynamics, for a marginally inhomogeneous initial
structure compared to macroscopically homogeneous initial structures, the rate
at which the structure equilibrates is much lower; This is probably because
imposing marginal inhomogeneity in a system under periodic boundary conditions
introduces long-wavelength inhomogeneous modes whose fluctuations are
suppressed by the Langevin thermostat thereby slowing equilibration. While
these results might seem intuitive, we have placed them on a solid analytical
foundation in this paper. In addition, we found initial structural
inhomogeneity makes the RDF approach one in time from above at large
distances. These effects of initial structural inhomogeneity in the Langevin
$NVT$ ensemble should be taken into account while investigating the structural
evolution of systems at near-zero and high densities in (MD) simulations. As
DPD and Langevin thermostats are similar [93], we would expect almost the same
behaviours for the DPD thermostat.

This paper showed that the structural relaxation time, unlike the temperature
relaxation time $\tau_{T}$, is an increasing function of the Langevin
constant, $\gamma$, in Langevin dynamics. In (MD) simulations, to ensure the
equilibration time $t_{e}$ is enough for the system to equilibrate at low/high
densities, one should plot $r[g(r)-1]$ for $r>\xi$ [$\xi$ denotes the
(effective) correlation length]. Any state- or thermostat-dependent
(non)linear behaviour in $r[g(r)-1]$ indicates the system has not equilibrated
properly and hence some nonequilibrium values might have leaked into the
equilibrium time averaging of dynamical quantities, which leads to unreliable
averaged equilibrium values.

The data that supports the findings of this study are available within the
article in the figures.

##  VII Acknowledgements

This work was supported by the Natural Sciences and Engineering Council of
Canada (NSERC). We would like to thank the Shared Hierarchical Academic
Research Computing Network (SHARCNET) and Compute/Calcul Canada for the
computational resources.

## References

  * Maxwell [1867] J. C. Maxwell, Philosophical Transactions of the Royal Society of London 157, 49 (1867).
  * Robinson and Frosch [1963] G. W. Robinson and R. P. Frosch, The Journal of Chemical Physics 38, 1187 (1963).
  * Jortner [1970] J. Jortner, Pure and Applied Chemistry 24, 165 (1970).
  * Abragam [1983] A. Abragam,  _The Principles of Nuclear Magnetism_ , International Series of Monographs on Physics (Clarendon Press, 1983).
  * Karmakar  _et al._ [2009] S. Karmakar, C. Dasgupta, and S. Sastry, Proceedings of the National Academy of Sciences 106, 3675 (2009).
  * Rauch and Tremaine [1996] K. P. Rauch and S. Tremaine, New Astronomy 1, 149 (1996).
  * Gujrati [2018] P. D. Gujrati, Entropy 20, 3675 (2018).
  * Samanta and Richert [2015] S. Samanta and R. Richert, The Journal of Chemical Physics 142, 044504 (2015).
  * Torre  _et al._ [2004] R. Torre, P. Bartolini, and R. Righini, Nature 428, 296 (2004).
  * Frenkel [1946] J. Frenkel,  _Kinetic Theory of Liquids_ , 1st ed. (Oxford at the Clarendon Press, 1946).
  * Santoro and Gorelli [2008] M. Santoro and F. A. Gorelli, Physical Review B 77, 212103 (2008).
  * Brazhkin  _et al._ [2012] V. V. Brazhkin, A. G. Lyapin, V. N. Ryzhov, K. Trachenko, Y. D. Fomin, and E. N. Tsiok, Physics-Uspekhi 55, 1061 (2012).
  * Trachenko and Brazhkin [2015] K. Trachenko and V. V. Brazhkin, Reports on Progress in Physics 79, 016502 (2015).
  * Fomin  _et al._ [2016] Y. D. Fomin, V. N. Ryzhov, E. N. Tsiok, V. V. Brazhkin, and K. Trachenko, Journal of Physics. Condensed Matter: an Institute of Physics journal 28, 43LT01 (2016).
  * Ghosh and Krishnamurthy [2018] K. Ghosh and C. V. Krishnamurthy, Physical Review E 97, 012131 (2018).
  * Allen and Tildesley [2017] M. P. Allen and D. J. Tildesley,  _Computer Simulation of Liquids_ , 2nd ed. (Oxford University Press, 2017).
  * Leimkuhler and Matthews [2015] B. Leimkuhler and C. Matthews,  _Molecular Dynamics: With Deterministic and Stochastic Numerical Methods_ , Interdisciplinary Applied Mathematics (Springer International Publishing, 2015).
  * Gao  _et al._ [2016] X. Gao, J. Fang, and H. Wang, The Journal of Chemical Physics 144, 124113 (2016).
  * Braun  _et al._ [2018] E. Braun, S. M. Moosavi, and B. Smit, Journal of Chemical Theory and Computation 14, 5262 (2018).
  * Davidchack [2010] R. L. Davidchack, Journal of Computational Physics 229, 9323 (2010).
  * Rosta  _et al._ [2009] E. Rosta, N. Buchete, and G. Hummer, Journal of Chemical Theory and Computation 5, 1393 (2009).
  * Basconi and Shirts [2013] J. E. Basconi and M. R. Shirts, Journal of Chemical Theory and Computation 9, 2887 (2013).
  * Rapaport [2004] D. C. Rapaport,  _The Art of Molecular Dynamics Simulation_, 2nd ed. (Cambridge University Press, 2004).
  * Kuksin  _et al._ [2005] A. Y. Kuksin, I. V. Morozov, G. E. Norman, V. V. Stegailov, and I. A. Valuev, Molecular Simulation 31, 1005 (2005).
  * Dunn  _et al._ [2016] J. Dunn, E. Antillon, J. Maassen, M. Lundstrom, and A. Strachan, Journal of Applied Physics 120, 225112 (2016).
  * Komatsu  _et al._ [2004] T. Komatsu, N. Yoshii, S. Miura, and S. Okazaki, Fluid Phase Equilibria 226, 345 (2004).
  * Yarnell  _et al._ [1973] J. L. Yarnell, M. J. Katz, R. G. Wenzel, and S. H. Koenig, Physical Review A 7, 2130 (1973).
  * Renninger and Averbach [1973] A. L. Renninger and B. L. Averbach, Physical Review B 8, 1507 (1973).
  * Gruebel and Clayton [1967] R. W. Gruebel and G. T. Clayton, The Journal of Chemical Physics 46, 639 (1967).
  * Kacner and Bartell [1979] M. A. Kacner and L. S. Bartell, The Journal of Chemical Physics 71, 192 (1979).
  * Fichthorn and Weinberg [1991] K. A. Fichthorn and W. H. Weinberg, The Journal of Chemical Physics 95, 1090 (1991).
  * Kodama  _et al._ [2006] K. Kodama, S. Iikubo, T. Taguchi, and S. Shamoto, Acta Crystallographica Section A 62, 444 (2006).
  * Román  _et al._ [1999] F. L. Román, A. González, J. A. White, and S. Velasco, American Journal of Physics 67, 1149 (1999).
  * Dawass  _et al._ [2018] N. Dawass, P. Krüger, S. K. Schnell, D. Bedeaux, S. Kjelstrup, J. M. Simon, and T. J. H. Vlugt, Molecular Simulation 44, 599 (2018).
  * Villamaina and Trizac [2014] D. Villamaina and E. Trizac, European Journal of Physics 35, 035011 (2014).
  * Kim and Yamamoto [2000] K. Kim and R. Yamamoto, Physical Review E 61, R41 (2000).
  * Dinnebier and Billinge [2008] R. E. Dinnebier and S. J. L. Billinge, eds.,  _Powder Diffraction: Theory and Practice_ , 1st ed. (The Royal Society of Chemistry, 2008).
  * Barrat and Hansen [2003] J. L. Barrat and J. P. Hansen,  _Basic Concepts for Simple and Complex Liquids_ (Cambridge University Press, 2003).
  * Zhang [2016] K. Zhang, arXiv:1606.03610 [cond-mat.soft] (2016).
  * Sturm [1993] K. Sturm, Zeitschrift Naturforschung Teil A 48, 233 (1993).
  * Hansen and McDonald [2013] J. P. Hansen and I. R. McDonald,  _Theory of Simple Liquids: With Applications to Soft Matter_ , 4th ed. (Academic Press, 2013).
  * Newman [1994] K. E. Newman, Chemistry Society Reviews 23, 31 (1994).
  * Chandler [1987] D. Chandler,  _Introduction to Modern Statistical Mechanics_ (Oxford University Press, 1987).
  * Lebowitz and Percus [1961] J. L. Lebowitz and J. K. Percus, Physical Review 124, 1673 (1961).
  * Hill [1956] T. Hill,  _Statistical Mechanics: Principles and Selected Applications_ (McGraw-Hill, 1956).
  * Salacuse  _et al._ [1996] J. J. Salacuse, A. R. Denton, and P. A. Egelstaff, Physical Review E 53, 2382 (1996).
  * Heidari  _et al._ [2018] M. Heidari, K. Kremer, R. Potestio, and R. Cortes-Huerto, Molecular Physics 116, 3301 (2018).
  * Román  _et al._ [1997] F. L. Román, J. A. White, and S. Velasco, The Journal of Chemical Physics 107, 4635 (1997).
  * Salacuse [2008] J. Salacuse, Physica A: Statistical Mechanics and its Applications 387, 3073 (2008).
  * Henderson [1994] J. R. Henderson, Physical Review E 50, 4836 (1994).
  * Hughes  _et al._ [2014] A. P. Hughes, U. Thiele, and A. J. Archer, American Journal of Physics 82, 1119 (2014).
  * Dijkstra and Evans [2000] M. Dijkstra and R. Evans, Journal of Chemical Physics 112, 1449 (2000).
  * Stopper  _et al._ [2019] D. Stopper, H. Hansen-Goos, R. Roth, and R. Evans, Journal of Chemical Physics 151, 014501 (2019).
  * Vega  _et al._ [1995] C. Vega, L. F. Rull, and S. Lago, Physical Review E 51, 3146 (1995).
  * Savenko and Dijkstra [2005] S. V. Savenko and M. Dijkstra, Physical Review E 72, 021202 (2005).
  * Vorontsov and Kuts [2008] A. G. Vorontsov and D. A. Kuts, Journal of Physics: Conference Series 98, 012004 (2008).
  * Evans and Henderson [2009] R. Evans and J. R. Henderson, Journal of Physics: Condensed Matter 21, 474220 (2009).
  * Montero [2017] A. M. Montero, arXiv:1710.01118 [cond-mat.soft] (2017).
  * Plischke and Bergersen [2006] M. Plischke and B. Bergersen,  _Equilibrium Statistical Physics_ , 3rd ed. (World Scientific, 2006).
  * Landau and Lifshitz [1980] L. Landau and E. Lifshitz,  _Statistical Physics, part 1_, 3rd ed. (Butterworth-Heinemann, 1980).
  * Green [1952] H. Green,  _The Molecular Theory of Fluids_ (North-Holland, Amsterdam, 1952).
  * Kirkwood and Boggs [1942] J. G. Kirkwood and E. M. Boggs, The Journal of Chemical Physics 10, 394 (1942).
  * Wallace [1987] D. C. Wallace, The Journal of Chemical Physics 87, 2282 (1987).
  * Laird and Haymet [1992] B. B. Laird and A. D. J. Haymet, Physical Review A 45, 5680 (1992).
  * Widom and Gao [2019] M. Widom and M. Gao, Entropy 21, 131 (2019).
  * Baranyai and Evans [1990] A. Baranyai and D. J. Evans, Physical Review A 42, 849 (1990).
  * Borzsák and Baranyai [1992] I. Borzsák and A. Baranyai, Chemical Physics 165, 227 (1992).
  * Baranyai and Evans [1989] A. Baranyai and D. J. Evans, Physical Review A 40, 3817 (1989).
  * Gallo and Rovere [2015] P. Gallo and M. Rovere, Physical Review E 91, 012107 (2015).
  * Reif [2009] F. Reif,  _Fundamentals of Statistical and Thermal Physics_ (Waveland Press, 2009).
  * Greiner  _et al._ [1995] W. Greiner, L. Neise, H. Stocker, H. Stöcker, and D. Rischke,  _Thermodynamics and Statistical Mechanics_ , Classical Theoretical Physics (Springer, 1995).
  * Levine [2008] I. N. Levine,  _Physical Chemistry_ , 6th ed. (McGraw-Hill Education, 2008).
  * McCoy and Davis [1979] B. F. McCoy and H. T. Davis, Physical Review A 20, 1201 (1979).
  * Waterson and Beaufort [1851] J. J. Waterson and F. Beaufort, Abstracts of the Papers Communicated to the Royal Society of London 5, 604 (1851).
  * de Groot and Mazur [1984] S. R. de Groot and P. Mazur,  _Non-equilibrium Thermodynamics_, Dover Books on Physics (Dover Publications, 1984).
  * Vilar and Rubí [2001] J. M. G. Vilar and J. M. Rubí, Proceedings of the National Academy of Sciences 98, 11081 (2001).
  * Rubi and Perez-Vicente [1997] M. Rubi and C. Perez-Vicente, eds.,  _Complex Behaviour of Glassy Systems_, Lecture Notes in Physics (Springer-Verlag Berlin Heidelberg, 1997).
  * Davidchack  _et al._ [2009] R. L. Davidchack, R. Handel, and M. V. Tretyakov, Journal of Chemical Physics 130, 234101 (2009).
  * Bencivenga  _et al._ [2009] F. Bencivenga, A. Cimatoribus, A. Gessini, M. G. Izzo, and C. Masciovecchio, The Journal of Chemical Physics 131, 144502 (2009).
  * Plimpton [1995] S. Plimpton, Journal of Computational Physics 117, 1 (1995).
  * Frenkel and Smit [2002] D. Frenkel and B. Smit,  _Understanding Molecular Simulation: From Algorithms to Applications_ , 2nd ed. (Academic Press, 2002).
  * Tester and Modell [1996] J. Tester and M. Modell,  _Thermodynamics and Its Applications_ , 3rd ed. (Prentice Hall, 1996).
  * Ray and Zhang [1999] J. R. Ray and H. Zhang, Physical Review E 59, 4781 (1999).
  * White  _et al._ [2008] J. A. White, F. L. Román, A. González, and S. Velasco, Physica A: Statistical Mechanics and its Applications 387, 6705 (2008).
  * Bussi and Parrinello [2008] G. Bussi and M. Parrinello, Computer Physics Communications 179, 26 (2008).
  * Jensen and Grønbech-Jensen [2019] L. F. G. Jensen and N. Grønbech-Jensen, Molecular Physics 117, 2511 (2019).
  * Grønbech-Jensen [2020] N. Grønbech-Jensen, Molecular Physics 118, e1662506 (2020).
  * Mannella [2006] R. Mannella, SIAM Journal on Scientific Computing 27, 2121 (2006).
  * Batrouni  _et al._ [1985] G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage, B. Svetitsky, and K. G. Wilson, Physical Review D 32, 2736 (1985).
  * Kutz [2013] J. N. Kutz,  _Data-Driven Modeling and Scientific Computation: Methods for Complex Systems and Big Data_ (Oxford University Press, 2013).
  * Levine  _et al._ [2011] B. G. Levine, J. E. Stone, and A. Kohlmeyer, Journal of Computational Physics 230, 3556 (2011).
  * Gardiner [2009] C. Gardiner,  _Stochastic Methods: A Handbook for the Natural and Social Science_ , 4th ed., Springer Series in Synergetics (Springer, 2009).
  * Pastorino  _et al._ [2007] C. Pastorino, T. Kreer, M. Müller, and K. Binder, Phys. Rev. E 76, 026706 (2007).