New_connection Old_connection

# On the Application of Non-Gaussian Noise in Stochastic Langevin Simulations

Niels Grønbech-Jensen<EMAIL_ADDRESS>Department of Mechanical &
Aerospace Engineering
Department of Mathematics
University of California, Davis, CA 95616, U.S.A.

###### Abstract

In light of recent advances in time step independent stochastic integrators
for Langevin equations, we revisit the considerations for using non-Gaussian
distributions for the thermal noise term in discrete-time thermostats. We find
that the desirable time step invariance of the modern methods is rooted in the
Gaussian noise, and that deviations from this distribution will distort the
Boltzmann statistics arising from the fluctuation-dissipation balance of the
integrators. We use the GJ methods as the focus of our investigation since
these methods are the ones that contain the most accurate thermodynamic
measures of existing methods. Within this set of methods we find that any
distribution of applied noise which satisfies the two first moments given by
the fluctuation-dissipation theorem, will result in correct, time step
independent results that are generated by the first two moments of the system
coordinates. However, if non-Gaussian noise is applied, undesired deviations
in higher moments of the system coordinates will appear to the detriment of
several important thermodynamic measures that depend especially on the fourth
moments. The deviations, induced by non-Gaussian noise, become significant
with the one time step attenuation, thereby inhibiting the benefits of the
methods. Thus, we conclude that the application of Gaussian noise is necessary
for reliable thermodynamic results when using modern stochastic thermostats
with large time steps.

##  I Introduction

Over the past several decades, discrete-time Langevin and Brownian simulations
in computational statistical mechanics have been a centerpiece in many
specific simulation applications to molecular dynamics, e.g., materials
science, soft matter, and biomolecular modeling AllenTildesley ; Frenkel ;
Rapaport ; Hoover_book ; Leach . Inherent to simulations is the fundamental
complication of time discretization, which on one hand is necessary for
advancing time, but on the other hand inevitably corrupts the desired,
continuous-time system features of interest. An integral part of a simulation
task is therefore to explore and optimize the balance between the two
conflicting core objectives; namely simulation efficiency by increasing the
time step, and simulation accuracy by decreasing it.

Since discrete-time algorithms for time evolution produce only approximate
solutions to a continuous-time solution, and since the errors must vanish in
the limit of infinitesimal time step, most considerations for a useful
simulation environment have been considered in the limit of small time steps,
where the simulated system variables change only little in each step. In that
limit, the discrete-time system behaves much like the continuous-time system,
and algorithm development has therefore mostly been conducted by formulating
objectives and measures in continuous-time whereafter time has been
discretized and executed with relatively small steps in order to retain
accuracy.

For stochastic differential equations, notably Langevin equations, an added
complication is that it is not a specific trajectory which is of interest, but
instead the statistics of trajectories that determines if the desired
continuous-time behavior has been adequately represented. The traditional
approach to these problems has therefore been to acquire accurate statistics
from accurate trajectories obtained from small time steps, where continuous-
time behavior is well approximated. However, with the recent advent of
discrete-time stochastic algorithms that can produce correct statistics from
inaccurate, large time step trajectories GJF1 ; 2GJ ; GJ , the details of,
e.g., stochastic noise correlation and distribution now become essential for
obtaining the advantages of the statistical accuracy for all meaningful time
steps. The significance of these details arises from the discretization of
time, as outlined in the following.

We are here concerned with the Langevin equation Langevin ; Langevin_Eq

$\displaystyle m\dot{v}+\alpha\dot{r}$ $\displaystyle=$ $\displaystyle
f+\beta\;,$ (1)

where $m$ is the mass of an object with spatial (configurational) coordinate
$r$ and velocity (kinetic) coordinate $v=\dot{r}$. The object is subjected to
a force $f=-\nabla E_{p}(r)$, with $E_{p}(r)$ being a potential energy surface
for $r$, and linear friction that is represented by the non-negative constant
$\alpha$. The fluctuation-dissipation relationship specifies that the thermal
fluctuations $\beta$ can be represented by a distribution for which the
following first two moments are given Parisi :

$\displaystyle\langle\beta(t)\rangle$ $\displaystyle=$ $\displaystyle 0$ (2a)
$\displaystyle\langle\beta(t)\beta(t^{\prime})\rangle$ $\displaystyle=$
$\displaystyle 2\alpha\,k_{B}T\,\delta(t-t^{\prime})\;,$ (2b)

where $\delta(t)$ is Dirac’s delta function signifying that the fluctuations
are temporally uncorrelated, $k_{B}$ is Boltzmann’s constant, and $T$ is the
thermodynamic temperature of the heat-bath, to which the system is coupled
through the friction constant $\alpha$. A feature in statistical mechanics is
that the fluctuation-dissipation theorem for the noise distribution $\beta(t)$
only specifies the two moments given in Eq. (2) Parisi ; Langevin_Eq ;
comment_symmetry . Thus, there is a great deal of freedom in choosing the
distribution of the noise in the Langevin equation such that correct
thermodynamics is represented. The reason that the physical (statistical)
results do not depend on this freedom is that continuous-time demands a new,
uncorrelated noise value $\beta(t)$ be contributed and accumulated at every
instant, resulting in a Gaussian outcome at any non-zero time-scale due to the
central limit theorem (see, e.g., Ref. central_limit ). As a result, the
sampled (one-dimensional) Maxwell-Boltzmann velocity distribution
$\rho_{k,e}(v)$ and the configurational Boltzmann distribution $\rho_{c,e}(r)$
always become Langevin_Eq

$\displaystyle\rho_{k,e}(v)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{m}{2\pi
k_{B}T}}\exp\left(-\frac{\frac{1}{2}mv^{2}}{k_{B}T}\right)$ (3a)
$\displaystyle\rho_{c,e}(r)$ $\displaystyle=$ $\displaystyle{\cal
C}_{c}\exp\left(-\frac{E_{p}(r)}{k_{B}T}\right)\,,$ (3b) with the constant
${\cal C}_{c}$ determined by the normalization of probability
$\displaystyle\underset{\rm all\;space}{\int}\\!\\!\\!\rho_{c,e}(r)\,dr$
$\displaystyle=$ $\displaystyle 1\,,$ (3c)

regardless of which temporally uncorrelated distribution $\beta(t)$ is at play
in Eq. (1), and for as long as Eq. (2) is satisfied. Similarly, it is given
that diffusive transport for a flat surface is represented by a diffusion
constant AllenTildesley

$\displaystyle D_{E}$ $\displaystyle=$
$\displaystyle\frac{1}{d}\lim_{t\rightarrow\infty}\frac{\langle[r(t)-r(0)]^{2}\rangle}{2t}$
(4) $\displaystyle=$ $\displaystyle\frac{k_{B}T}{\alpha}\,,$ (5)

where $d$ is the dimension of the configurational space.

Given that the central limit theorem ensures a Gaussian noise distribution for
any appreciable integration of the instantaneous variable $\beta(t)$, it is
natural to choose $\beta(t)$ to be the Gaussian that conforms to Eq. (2).
However, in discrete-time simulations it is not uncommon to directly employ a
uniform distribution that satisfies Eq. (2). The rationale is rooted in the
computational efficiency of uniformly distributed pseudo-random generators
random , compared with the additional computational load required for
transforming those variables into Gaussian distributed variables before use. A
thorough investigation of using non-Gaussian valuables in discrete time was
presented in the 1988 publication Ref. Greiner_1 , where it was concluded
that, given the inherent inaccuracies of the stochastic discrete-time
integrators of the time, the scaling in statistical moments upon the applied
time step did not change adversely when using non-Gaussian noise
distributions, specifically the uniform. Building upon this work, Ref.
Dunweg_1  extended the investigation in 1991 to a broader class of noise
distributions, which would yield indistinguishable statistical results
compared to the use of Gaussian noise for vanishing time steps, where the
underlying systemic time step errors of the then contemporary methods would
vanish as well. The studies highlighted the benefits of computational
efficiencies in using simple, non-Gaussian noise given that the statistical
differences were deemed to be insignificant for small time steps, which in any
case were necessary for the accuracy of the stochastic integrators, even using
Gaussian noise. As pointed out in a recent article Finkelstein_1  that
compares the quality of contemporary stochastic integrators, this practice has
been widespread when using traditional algorithms, such as the Brünger-Brooks-
Karplus BBK  (also analyzed in Ref. Pastor_88 ) or the Schneider-Stoll SS
method, for as long as small time steps are applied (for a more expansive list
of traditional and recent Langevin integrators, please see Ref. GJ ). For
example, these two methods are used with non-Gaussian noise in, e.g., the
wide-distribution molecular dynamics simulation suites, LAMMPS Plimpton ;
Plimpton2 ; LAMMPS_uniform  and ESPResSo Espresso_4.2.0 .

It is the aim of this work to revisit the discussion of applying non-Gaussian
noise in light of the most modern stochastic integrators, which have been
designed to produce accurate, time step independent statistics and transport
through the use of Gaussian noise. We demonstrate that even if the low order
moments of the configurational and kinetic distributions are seemingly
unaffected by the details of the noise distribution, there may still be
significant differences in the sampled phase space that affect the statistical
sampling of the system. We base the analysis exclusively on the recently
formulated GJ set of methods GJ , since these methods have the most optimized
statistical properties that are invariant to the time step for linear systems
when using Gaussian noise. This also allows us to investigate the large time
step implications of using non-Gaussian noise without being subject to other
inherent time step errors of the integrators. The presentation is organized as
follows. Section II reviews the basic features and structure of the GJ
methods, Sec. III analyzes analytically the first four statistical moments of
the configurational and kinetic coordinates for linear systems, Sec. IV
presents companion simulations to the analysis conducted in Sec. III, and Sec.
V summarizes the discussion with specific conclusions.

##  II Review of the GJ methods

We adopt the complete set of GJ methods GJ  for solving Eq. (1) in the form

$\displaystyle u^{n+\frac{1}{2}}$ $\displaystyle=$
$\displaystyle\sqrt{c_{1}}\,v^{n}+\frac{\sqrt{c_{3}}\,\Delta{t}}{2m}\,f^{n}+\frac{\sqrt{c_{3}}}{2m}\,\beta^{n+1}$
(6a) $\displaystyle r^{n+1}$ $\displaystyle=$ $\displaystyle
r^{n}+\sqrt{c_{3}}\,\Delta{t}\,u^{n+\frac{1}{2}}$ (6b) $\displaystyle v^{n+1}$
$\displaystyle=$
$\displaystyle\frac{c_{2}}{\sqrt{c_{1}}}\,u^{n+\frac{1}{2}}+\sqrt{\frac{c_{3}}{c_{1}}}\,\frac{\Delta{t}}{2m}\,f^{n+1}+\sqrt{\frac{c_{3}}{c_{1}}}\,\frac{1}{2m}\,\beta^{n+1}\,,$

where a superscript $n$ on the variables $r$, $v$, $u$, and $f$ pertains to
the discrete time $t_{n}=t_{0}+n\,\Delta{t}$ at which the numerical
approximation is given, with $\Delta{t}$ being the discrete time step. Thus,
$f^{n}=f(t_{n},r^{n})$. The two represented velocities are, respectively, the
half-step velocity $u^{n+\frac{1}{2}}$ at time $t_{n+\frac{1}{2}}$, and the
on-site velocity $v^{n}$ at $t_{n}$, given by

$\displaystyle u^{n+\frac{1}{2}}$ $\displaystyle=$
$\displaystyle\frac{r^{n+1}-r^{n}}{\sqrt{c_{3}}\,\Delta{t}}$ (7a)
$\displaystyle v^{n}$ $\displaystyle=$
$\displaystyle\frac{r^{n+1}-(1-c_{2})r^{n}-c_{2}r^{n-1}}{2\sqrt{c_{1}c_{3}}\,\Delta{t}}+\frac{1}{4m}(\beta^{n}-\beta^{n+1})\,.$

The statistical definitions of what it means to be half-step and on-site are
given in Ref. GJ . The functional parameter
$c_{2}=c_{2}(\frac{\alpha\Delta{t}}{m})$ is the one-time-step attenuation,
which is a function of $\alpha\Delta{t}/m$. The two other functional
parameters, $c_{1}$ and $c_{3}$, are given by

$\displaystyle 2c_{1}$ $\displaystyle=$ $\displaystyle 1+c_{2}$ (8a)
$\displaystyle\frac{\alpha\Delta{t}}{m}\,c_{3}$ $\displaystyle=$
$\displaystyle 1-c_{2}\,.$ (8b)

It is the choice of the function $c_{2}$ that distinguishes the GJ methods
from each other (see Ref. GJ  for several choices of $c_{2}$). The basic
requirement for $c_{2}$ is that $c_{2}\rightarrow 1-\frac{\alpha\Delta{t}}{m}$
for $\frac{\alpha\Delta{t}}{m}\rightarrow 0$, ensuring that also
$c_{1}\rightarrow 1$ and $c_{3}\rightarrow 1$ for
$\frac{\alpha\Delta{t}}{m}\rightarrow 0$ in accordance with the frictionless
Verlet method Verlet ; Swope ; Beeman ; Buneman ; Hockney . Finally, the
discrete-time stochastic noise variable is defined by

$\displaystyle\beta^{n+1}$ $\displaystyle=$
$\displaystyle\int_{t_{n}}^{t_{n+1}}\beta(t)\,dt\,,$ (9)

such that the discrete-time fluctuation-dissipation theorem follows from Eqs.
(9) and (2) to become

$\displaystyle\langle\beta^{n}\rangle$ $\displaystyle=$ $\displaystyle 0$
(10a) $\displaystyle\langle\beta^{n}\beta^{\ell}\rangle$ $\displaystyle=$
$\displaystyle 2\alpha\Delta{t}\,k_{B}T\,\delta_{n,\ell}\,,$ (10b)

where $\delta_{n,\ell}$ is Kronecker’s delta function.

Eliminating the velocities from Eq. (6) yields the purely configurational
stochastic Verlet form

$\displaystyle r^{n+1}$ $\displaystyle=$ $\displaystyle
2c_{1}\,r^{n}-c_{2}\,r^{n-1}+\frac{c_{3}\Delta{t}^{2}}{m}\,f^{n}+\frac{c_{3}\Delta{t}}{2m}(\beta^{n}+\beta^{n+1})\,,$

which is the only stochastic Verlet form that 1) for a harmonic potential
$E_{p}(r)=-\kappa r$ ($\kappa>0$) and Gaussian fluctuations $\beta^{n}$, will
yield the correct Boltzmann distribution $\rho_{c,e}(r^{n})$ given by Eq.
(3b), such that

$\displaystyle\rho_{c,e}(r^{n})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\kappa}{2\pi
k_{B}T}}\exp\left(-\frac{\frac{1}{2}\kappa(r^{n})^{2}}{k_{B}T}\right)$ (12a)
$\displaystyle\langle r^{n}r^{n}\rangle$ $\displaystyle=$
$\displaystyle\frac{k_{B}T}{\kappa}\,,$ (12b)

for any time step within the stability range
($\Omega_{0}\Delta{t}<2\sqrt{c_{1}/c_{3}}$); 2) for a linear potential $f={\rm
const}$, yield the correct configurational drift velocity

$\displaystyle
v_{d}\;=\;\left\langle\frac{r^{n+1}-r^{n}}{\Delta{t}}\right\rangle$
$\displaystyle=$ $\displaystyle\frac{f}{\alpha}\,;$ (13)

and 3) for $f=0$, yield the correct Einstein diffusion given by Eq. (4), such
that

$\displaystyle D_{E}$ $\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{\left\langle(r^{n}-r^{0})^{2}\right\rangle}{2d\,n\Delta{t}}\;=\;\frac{k_{B}T}{\alpha}\,.$
(14)

Similarly, it was found GJ  that for a single Gaussian noise value $\beta^{n}$
per time step, the harmonic potential $E_{p}(r)=-\kappa r$ ($\kappa\geq 0$)
results in the unique half-step velocity $u^{n+\frac{1}{2}}$ in Eq. (7a) with
correct velocity distribution function $\rho_{k,e}(u^{n+\frac{1}{2}})$ from
Eq. (3a) such that

$\displaystyle\langle u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle$
$\displaystyle=$ $\displaystyle\frac{k_{B}T}{m}\,.$ (15)

In contrast, it was also found in Ref. GJ  that is not possible to define an
onsite velocity $v^{n}$ with the same time step independent property. Given
the statistical invariance of these equations upon the time step when using
Gaussian noise, we will now investigate the statistical properties of $r^{n}$
and $u^{n+\frac{1}{2}}$ when using other distributions than Gaussian for
$\beta^{n}$.

##  III Moment Analysis for Linear Systems

We will here investigate the relevant moments of the configurational
coordinate $r^{n}$ and half-step velocity $u^{n+\frac{1}{2}}$ as a function of
the applied noise distribution $\rho(\beta)$ of $\beta^{n}$, which satisfies
Eq. (10). For this analysis we will mostly assume a linear system with a
harmonic potential

$\displaystyle E_{p}(r)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\kappa\,r^{2}$ (16) $\displaystyle\Rightarrow\;\;f$
$\displaystyle=$ $\displaystyle-\kappa r\,,$ (17)

such that the linearized configurational stochastic GJ methods in Eq.
(LABEL:eq:GJ_c) can be written

$\displaystyle r^{n+1}$ $\displaystyle=$ $\displaystyle
2c_{1}X\,r^{n}-c_{2}\,r^{n-1}+\frac{c_{3}\Delta{t}}{2m}(\beta^{n}+\beta^{n+1})\,.$
(18)

with

$\displaystyle X$ $\displaystyle=$ $\displaystyle
1-\frac{c_{3}}{c_{1}}\frac{\Omega_{0}^{2}\Delta{t}^{2}}{2}\,.$ (19)

The natural frequency $\Omega_{0}$ of the oscillator is given by
$\Omega_{0}^{2}=\kappa/m$.

###  III.1 First and Second Moments

For the first moments of the variables, it is obvious that Eq. (10a) ensures
both $\langle r^{n}\rangle=0$ (from Eq. (18)) and subsequently $\langle
u^{n+\frac{1}{2}}\rangle=0$ (from Eq. (7a)) regardless of which distribution
is chosen for $\beta^{n}$. It similarly follows that Eq. (10a) also ensures
that the drift velocity in Eq. (13) for $f={\rm const}$ remains constant and
correct regardless of both time step and distribution $\beta^{n}$.

For the second moments, $\langle r^{n}r^{n}\rangle$ and $\langle
u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle$, we first find $\langle
r^{n}r^{n}\rangle$ as the solution to the system of equations given by the
statistical averages of Eq. (18) multiplied by, respectively $r^{n-1}$,
$r^{n}$, and $r^{n+1}$. This linear system of equations can be conveniently
written

$\displaystyle\left(\begin{array}[]{ccc}1&-2c_{1}X&c_{2}\\\
0&2c_{1}&-2c_{1}X\\\
c_{2}&-2c_{1}X&1\end{array}\right)\left(\begin{array}[]{c}\langle
r^{n-1}r^{n+1}\rangle\\\ \langle r^{n}r^{n+1}\rangle\\\ \langle
r^{n}r^{n}\rangle\end{array}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}0\\\ 1\\\
2c_{1}X+2\end{array}\right)\,\left(\frac{c_{3}\Delta{t}}{2m}\right)^{2}\,\langle\beta^{n}\beta^{n}\rangle\,.$
(29)

We immediately observe that this equation is unaffected by the details of the
distribution of $\beta^{n}$ for as long as Eq. (10) is satisfied. Thus, also
the second moment, $\langle r^{n}r^{n}\rangle$, is unaffected by the chosen
distribution of $\beta^{n}$, regardless of the time step $\Delta{t}$. The
corresponding second moment of the half-step velocity in Eq. (7a) is then
given by

$\displaystyle\langle u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle$
$\displaystyle=$ $\displaystyle\frac{2}{c_{3}\,\Delta{t}^{2}}(\langle
r^{n}r^{n}\rangle-\langle r^{n}r^{n+1}\rangle)\,,$ (30)

where both correlations on the right hand side are solutions found from Eq.
(29), which is unchanged for as long as Eq. (10b) is satisfied. Thus, the
second moment of the half-step velocity is also unaffected by the chosen
distribution of $\beta^{n}$, regardless of the time step $\Delta{t}$. The
averages of potential (Eq. (16)) and kinetic energies,

$\displaystyle\left\langle E_{p}(r^{n})\right\rangle$ $\displaystyle=$
$\displaystyle\frac{\kappa}{2}\langle r^{n}r^{n}\rangle\;=\;\frac{1}{2}k_{B}T$
(31) $\displaystyle\left\langle E_{k}(u^{n+\frac{1}{2}})\right\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2}\frac{m}{c_{3}\,\Delta{t}^{2}}(\langle
r^{n}r^{n}\rangle-\langle r^{n}r^{n+1}\rangle)$ (32) $\displaystyle=$
$\displaystyle\frac{1}{2}k_{B}T\,,$

are therefore unaffected for any applied stochastic variable $\beta^{n}$
satisfying Eq. (10).

Since the second moments of the coordinates $r^{n}$ and $u^{n+\frac{1}{2}}$ of
the GJ methods are invariant to the specific choice of noise distribution that
satisfies Eq. (10), we conclude that the GJ methods display time-step
invariance in both potential and kinetic energy, and therefore temperature,
regardless of the applied $\beta^{n}$ distribution. The invariance also
extends to the above-mentioned diffusion Eqs. (4) and (14), since this
quantity is rooted in only second moments of the noise variable.

Numerical validation of the time step independence of the second moments is
shown in Figs. 2ab and 6ab for a couple of non-Gaussian noise variables, and
simulations are discussed in Sec. IV. While the invariance of the second
moments against the choice of $\beta^{n}$ distribution is appealing, as it
gives the impression that temperature is correctly achieved, it may not be
sufficient to ensure proper thermodynamics sampling as we will see in the
following.

###  III.2 Third Moments

The third moments of $r^{n}$ and $u^{n+\frac{1}{2}}$ are generally not needed
for the calculations of thermodynamic quantities, especially not if the system
potentials are symmetric. However, they may be relevant to study in case
asymmetric fluctuations are applied to the system for some reason. Following
the same procedure that was used for calculating the second moments, we form
the statistical averages of the product of Eq. (18) with the six unique
combinations of $r^{\ell_{1}}r^{\ell_{2}}$, where $\ell_{i}=n-1,n,n+1$. The
result is the linear system

$\displaystyle{\cal A}_{3}\left(\begin{array}[]{c}\langle
r^{n-1}r^{n-1}r^{n+1}\rangle\\\ \langle r^{n-1}r^{n}r^{n+1}\rangle\\\ \langle
r^{n-1}r^{n+1}r^{n+1}\rangle\\\ \langle r^{n}r^{n}r^{n+1}\rangle\\\ \langle
r^{n}r^{n+1}r^{n+1}\rangle\\\ \langle
r^{n}r^{n}r^{n}\rangle\end{array}\right)$ $\displaystyle=$
$\displaystyle\left(\frac{k_{B}T}{\kappa}\right)^{\frac{3}{2}}\,{\cal
R}\;=\;\left(\begin{array}[]{c}\langle r^{n-1}r^{n-1}\beta^{n}\rangle\\\
\langle r^{n-1}r^{n}\beta^{n}\rangle\\\ \langle
r^{n-1}r^{n+1}(\beta^{n}+\beta^{n+1})\rangle\\\ \langle
r^{n}r^{n}\beta^{n}\rangle\\\ \langle
r^{n}r^{n+1}(\beta^{n}+\beta^{n+1})\rangle\\\ \langle
r^{n+1}r^{n+1}(\beta^{n}+\beta^{n+1})\rangle\end{array}\right)\frac{c_{3}\Delta{t}}{2m}\,,$
(45)

where

$\displaystyle{\cal A}_{3}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}1&0&0&-2c_{1}X&0&c_{2}\\\
0&1&0&c_{2}&-2c_{1}X&0\\\ c_{2}&-2c_{1}X&1&0&0&0\\\ 0&0&0&1&c_{2}&-2c_{1}X\\\
0&c_{2}&0&-2c_{1}X&1&0\\\ 0&0&c_{2}&0&-2c_{1}X&1\end{array}\right)\,.$ (52)

The right hand side ${\cal R}=\\{{\cal R}_{i}\\}$ can after some calculations
be expressed as

$\displaystyle{\cal R}_{1}$ $\displaystyle=$ $\displaystyle 0$ (53a)
$\displaystyle{\cal R}_{2}$ $\displaystyle=$ $\displaystyle 0$ (53b)
$\displaystyle{\cal R}_{3}$ $\displaystyle=$ $\displaystyle 0$ (53c)
$\displaystyle{\cal R}_{4}$ $\displaystyle=$
$\displaystyle(1-c_{2})^{3}\left(\frac{\kappa\Delta{t}}{\alpha}\right)^{\frac{3}{2}}\,\Gamma_{3}$
(53d) $\displaystyle{\cal R}_{5}$ $\displaystyle=$
$\displaystyle(2c_{1}X+1){\cal R}_{4}$ (53e) $\displaystyle{\cal R}_{6}$
$\displaystyle=$ $\displaystyle(2c_{1}X(2c_{1}X+2)+2){\cal R}_{4}\,,$ (53f)

where the fluctuations $\beta^{n}$ satisfy Eq. (10), and the skewness is
denoted by $\Gamma_{3}$,

$\displaystyle\Gamma_{3}$ $\displaystyle=$
$\displaystyle\frac{\left\langle\beta^{n}\beta^{n}\beta^{n}\right\rangle}{\left\langle\beta^{n}\beta^{n}\right\rangle^{\frac{3}{2}}}\,.$
(54)

As is the case when using any symmetric distribution for $\beta^{n}$, all odd
moments for the configurational coordinate $r^{n}$ will be zero for a Gaussian
distribution of $\beta^{n}$ applied to a symmetric potential, $E_{p}(r)$.

The third moment of the velocity $u^{n+\frac{1}{2}}$ is given by

$\displaystyle\left\langle
u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\right\rangle\;=\;\frac{1}{\left(\Delta{t}\,\sqrt{c_{3}}\right)^{3}}\left\langle(r^{n+1}-r^{n})^{3}\right\rangle$
$\displaystyle=\;\frac{3}{\left(\Delta{t}\,\sqrt{c_{3}}\right)^{3}}\left(\left\langle
r^{n}r^{n}r^{n+1}\right\rangle-\left\langle
r^{n}r^{n+1}r^{n+1}\right\rangle\right)\,,$ (55)

where the two correlations on the right hand side are part of the solution to
Eq. (45).

For any given choice of GJ method ($c_{2}$) the above system of equations can
be easily solved numerically for given system parameters and skewness
$\Gamma_{3}$. Examples of this will be discussed in Sec. IV, and are shown in
Fig. 6cd along with the results of corresponding numerical solution of Eqs.
(6) and (17) for a particular asymmetric noise distribution. Closed form
expressions are obtained in the following for a couple of special cases.

Special Case, $X=0$

The matrix ${\cal A}_{3}$ and the vector ${\cal R}$ simplify considerably with
easy accessible solutions

$\displaystyle\langle r^{n}r^{n}r^{n}\rangle$ $\displaystyle=$
$\displaystyle\frac{2}{1+c_{2}^{3}}\left(\frac{k_{B}T}{\kappa}\right)^{\frac{3}{2}}\,{\cal
R}_{4}$ (56a) $\displaystyle\langle r^{n}r^{n}r^{n+1}\rangle$ $\displaystyle=$
$\displaystyle\frac{1-c_{2}}{1+c_{2}^{3}}\left(\frac{k_{B}T}{\kappa}\right)^{\frac{3}{2}}\,{\cal
R}_{4}$ (56b) $\displaystyle\langle r^{n}r^{n+1}r^{n+1}\rangle$
$\displaystyle=$
$\displaystyle\frac{1+c_{2}^{2}}{1+c_{2}^{3}}\left(\frac{k_{B}T}{\kappa}\right)^{\frac{3}{2}}\,{\cal
R}_{4}\,,$ (56c)

where we have used that this special case for
$\Omega_{0}^{2}\Delta{t}^{2}=2\frac{c_{3}}{c_{1}}$ offers ${\cal R}_{6}=2{\cal
R}_{5}=2{\cal R}_{4}$, and where

$\displaystyle{\cal R}_{4}$ $\displaystyle=$
$\displaystyle\left(\frac{1-c_{2}^{2}}{2}\right)^{\frac{3}{2}}\,\Gamma_{3}\,.$
(57)

Equations (55), (56b), and (56c) give

$\displaystyle\left\langle\left(u^{n+\frac{1}{2}}\right)^{3}\right\rangle\;=\;-\frac{3}{\left(\Delta{t}\sqrt{c_{3}}\right)^{3}}\frac{c_{2}}{c_{2}^{2}-c_{2}+1}\left(\frac{k_{B}T}{\kappa}\right)^{\frac{3}{2}}{\cal
R}_{4}$
$\displaystyle=\;-\frac{3}{2\sqrt{2}}\,\left(\frac{1-c_{2}^{2}}{1+c_{2}}\right)^{\frac{3}{2}}\frac{c_{2}}{c_{2}^{2}-c_{2}+1}\,\left(\frac{k_{B}T}{m}\right)^{\frac{3}{2}}\,\Gamma_{3}\,.$
(58)

Examples of the use of these expressions are shown in Fig. 6cd as open markers
($\circ$) for a particular asymmetric noise distribution.

Special Case, $c_{2}=0$

In this case we have that $c_{1}=\frac{1}{2}$, $c_{3}=m/\alpha\Delta{t}$, and
$X=1-\Delta{t}\kappa/\alpha$. The simplified matrix ${\cal A}_{3}$ and vector
${\cal R}$ yield

$\displaystyle\langle r^{n}r^{n}r^{n}\rangle$ $\displaystyle=$
$\displaystyle\frac{3X^{2}+3X+2}{1-X^{3}}\left(\frac{k_{B}T}{\kappa}\right)^{\frac{3}{2}}\,{\cal
R}_{4}$ (59a) $\displaystyle\langle r^{n}r^{n}r^{n+1}\rangle$ $\displaystyle=$
$\displaystyle\frac{2X^{3}+3X^{2}+2X+1}{1-X^{3}}\left(\frac{k_{B}T}{\kappa}\right)^{\frac{3}{2}}{\cal
R}_{4}$ (59b) $\displaystyle\langle r^{n}r^{n+1}r^{n+1}\rangle$
$\displaystyle=$
$\displaystyle\frac{X^{4}+2X^{3}+2X^{2}+2X+1}{1-X^{3}}\left(\frac{k_{B}T}{\kappa}\right)^{\frac{3}{2}}{\cal
R}_{4}\,,$ (59c)

where

$\displaystyle{\cal R}_{4}$ $\displaystyle=$
$\displaystyle\left(\frac{\kappa\Delta{t}}{2\alpha}\right)^{\frac{3}{2}}\,\Gamma_{3}\,.$
(60)

Equations (55), (59b), and (59c) give

$\displaystyle\left\langle\left(u^{n+\frac{1}{2}}\right)^{3}\right\rangle$
$\displaystyle=$
$\displaystyle\frac{3X^{2}}{2\sqrt{2}}\frac{1-X^{2}}{1-X^{3}}\left(\frac{k_{B}T}{m}\right)^{\frac{3}{2}}\Gamma_{3}\,.$
(61)

Examples of the use of these expressions are shown in Fig. 6cd as closed
markers ($\bullet$) for a particular asymmetric noise distribution.

###  III.3 Fourth Moments

The fourth moments of $r^{n}$ and $u^{n+\frac{1}{2}}$ are needed for the
calculations of energy and pressure fluctuations, which are the basis of
several important thermodynamic quantities, such as heat capacity, thermal
expansion, compressibility, and the thermal pressure coefficient
AllenTildesley . Fourth moments are therefore critical for proper
thermodynamic characterization of a system. Following the same procedure that
was used for calculating the second and third moments, we form the statistical
averages of the product of Eq. (18) with the ten unique combinations of
$r^{\ell_{1}}r^{\ell_{2}}r^{\ell_{3}}$, where $\ell_{i}=n-1,n,n+1$. The result
is the linear system

$\displaystyle{\cal A}_{4}\left(\begin{array}[]{c}\langle
r^{n-1}r^{n-1}r^{n-1}r^{n+1}\rangle\\\ \langle
r^{n-1}r^{n-1}r^{n}r^{n+1}\rangle\\\ \langle
r^{n-1}r^{n-1}r^{n+1}r^{n+1}\rangle\\\ \langle
r^{n-1}r^{n}r^{n}r^{n+1}\rangle\\\ \langle
r^{n-1}r^{n}r^{n+1}r^{n+1}\rangle\\\ \langle
r^{n-1}r^{n+1}r^{n+1}r^{n+1}\rangle\\\ \langle
r^{n}r^{n}r^{n}r^{n+1}\rangle\\\ \langle r^{n}r^{n}r^{n+1}r^{n+1}\rangle\\\
\langle r^{n}r^{n+1}r^{n+1}r^{n+1}\rangle\\\ \langle
r^{n}r^{n}r^{n}r^{n}\rangle\end{array}\right)$ $\displaystyle=$
$\displaystyle\left(\frac{k_{B}T}{\kappa}\right)^{2}{\cal
R}\;=\;\left(\begin{array}[]{c}0\\\ \langle
r^{n-1}r^{n-1}r^{n}\beta^{n}\rangle\\\ \langle
r^{n-1}r^{n-1}r^{n+1}(\beta^{n}+\beta^{n+1})\rangle\\\ \langle
r^{n-1}r^{n}r^{n}\beta^{n}\rangle\\\ \langle
r^{n-1}r^{n}r^{n+1}(\beta^{n}+\beta^{n+1})\rangle\\\ \langle
r^{n-1}r^{n+1}r^{n+1}(\beta^{n}+\beta^{n+1})\rangle\\\ \langle
r^{n}r^{n}r^{n}\beta^{n}\rangle\\\ \langle
r^{n}r^{n}r^{n+1}(\beta^{n}+\beta^{n+1})\rangle\\\ \langle
r^{n}r^{n+1}r^{n+1}(\beta^{n}+\beta^{n+1})\rangle\\\ \langle
r^{n+1}r^{n+1}r^{n+1}(\beta^{n}+\beta^{n+1})\rangle\end{array}\right)\frac{c_{3}\Delta{t}}{2m}\,,$
(82)

where

$\displaystyle{\cal A}_{4}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccccccc}1&0&0&0&0&0&-2c_{1}X&0&0&c_{2}\\\
0&1&0&0&0&0&c_{2}&-2c_{1}X&0&0\\\ c_{2}&-2c_{1}X&1&0&0&0&0&0&0&0\\\
0&0&0&1&0&0&0&c_{2}&-2c_{1}X&0\\\ 0&c_{2}&0&-2c_{1}X&1&0&0&0&0&0\\\
0&0&c_{2}&0&-2c_{1}X&1&0&0&0&0\\\ 0&0&0&0&0&0&1&0&c_{2}&-2c_{1}X\\\
0&0&0&c_{2}&0&0&-2c_{1}X&1&0&0\\\ 0&0&0&0&c_{2}&0&0&-2c_{1}X&1&0\\\
0&0&0&0&0&c_{2}&0&0&-2c_{1}X&1\end{array}\right)\,.$ (93)

The right hand side ${\cal R}=\\{{\cal R}_{i}\\}$ can after some cumbersome
work be expressed as

$\displaystyle{\cal R}_{1}$ $\displaystyle=$ $\displaystyle 0$ (94a)
$\displaystyle{\cal R}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1-c_{2}^{2})(1-X)$ (94b) $\displaystyle{\cal R}_{3}$
$\displaystyle=$ $\displaystyle[2+2c_{1}X]{\cal R}_{2}$ (94c)
$\displaystyle{\cal R}_{4}$ $\displaystyle=$
$\displaystyle[2c_{1}X+1-c_{2}]{\cal R}_{2}$ (94d) $\displaystyle{\cal R}_{5}$
$\displaystyle=$ $\displaystyle[2c_{1}X(2c_{1}X+2-c_{2})+1-2c_{2}]{\cal
R}_{2}$ (94e) $\displaystyle{\cal R}_{6}$ $\displaystyle=$
$\displaystyle[(2c_{1}X)^{3}+(2c_{1}X)^{2}(3-c_{2})+4c_{1}X(1-2c_{2})-4c_{2}]{\cal
R}_{2}$ (94f) $\displaystyle{\cal R}_{7}$ $\displaystyle=$
$\displaystyle[3+(\Gamma_{4}-3){\cal R}_{2}]{\cal R}_{2}$ (94g)
$\displaystyle{\cal R}_{8}$ $\displaystyle=$
$\displaystyle[4c_{1}X+3-c_{2}+(2c_{1}X+1)(\Gamma_{4}-3){\cal R}_{2}]{\cal
R}_{2}$ (94h) $\displaystyle{\cal R}_{9}$ $\displaystyle=$
$\displaystyle[(2c_{1}X)^{2}+2c_{1}X(3-c_{2})+3-2c_{2}+(1+2c_{1}X)^{2}(\Gamma_{4}-3){\cal
R}_{2}]{\cal R}_{2}$ (94i) $\displaystyle{\cal R}_{10}$ $\displaystyle=$
$\displaystyle[3(2c_{1}X+2)+(1+(1+2c_{1}X)^{3})(\Gamma_{4}-3){\cal
R}_{2}]{\cal R}_{2}\,,$ (94j)

where we have used Eq. (10b). Notice that ${\cal R}_{7}$-${\cal R}_{10}$ in
Eqs. (94g)-(94j) depend on the kurtosis $\Gamma_{4}$

$\displaystyle\Gamma_{4}$ $\displaystyle=$
$\displaystyle\frac{\langle\beta^{n}\beta^{n}\beta^{n}\beta^{n}\rangle}{\langle\beta^{n}\beta^{n}\rangle^{2}}\,,$
(95)

which is not determined by the fluctuation-dissipation relationship in Eq.
(10).

While Eqs. (82)-(94) are not immediately informative, they provide a
straightforward path to calculating the sought-after fluctuation $\sigma_{p}$
in potential energy of Eq. (16)

$\displaystyle\sigma_{p}$ $\displaystyle=$ $\displaystyle\sqrt{\left\langle
E^{2}_{p}(r^{n})\right\rangle-\left\langle E_{p}(r^{n})\right\rangle^{2}}$
(96) $\displaystyle=$ $\displaystyle\frac{\kappa}{2}\sqrt{\langle
r^{n}r^{n}r^{n}r^{n}\rangle-\langle r^{n}r^{n}\rangle^{2}}\,.$

Equations (82)-(94) also lead directly to the fluctuations $\sigma_{k}$ in
kinetic energy of the velocity from Eq. (7a)

$\displaystyle\sigma_{k}$ $\displaystyle=$ $\displaystyle\sqrt{\left\langle
E^{2}_{k}(u^{n+\frac{1}{2}})\right\rangle-\left\langle
E_{k}(u^{n+\frac{1}{2}})\right\rangle^{2}}$ $\displaystyle=$
$\displaystyle\frac{m}{2}\sqrt{\langle
u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle-\langle
u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle^{2}}\,,$

where $\langle u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle$ is given by Eq.
(15), and

$\displaystyle\langle
u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle\;=\;\frac{1}{c_{3}^{2}\Delta{t}^{4}}\langle(r^{n+1}-r^{n})^{4}\rangle$
$\displaystyle=\;\frac{2}{c_{3}^{2}\Delta{t}^{4}}(\langle
r^{n}r^{n}r^{n}r^{n}\rangle-2\langle r^{n}r^{n}r^{n}r^{n+1}\rangle$
$\displaystyle+\;3\langle r^{n}r^{n}r^{n+1}r^{n+1}\rangle-2\langle
r^{n}r^{n+1}r^{n+1}r^{n+1}\rangle)\,.$ (98)

All the correlations on the right hand side are given from the solution to
Eqs. (82)-(94).

For any given choice of GJ method ($c_{2}$) the above system of equations can
be easily solved numerically for given system parameters and kurtosis
$\Gamma_{4}$. Examples of this will be discussed in Sec. IV, and are shown in
Figs. 2cd and 6ef along with the results of corresponding numerical solution
of Eqs. (6) and (17) for particular noise distributions. Closed form
expressions can be given in the special cases of $X=0$ and $c_{2}=0$.

Special Case, $X=0$

This particular time step is

$\displaystyle\Omega_{0}\Delta{t}$ $\displaystyle=$
$\displaystyle\sqrt{2}\sqrt{\frac{c_{1}}{c_{3}}}\,.$ (99)

Cumbersome, yet straightforward, algebra leads to the closed expressions for
$\langle r^{n}r^{n}r^{n}r^{n}\rangle$ and $\langle
u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle$
for $X=0$; i.e., for $\Omega_{0}\Delta{t}=\sqrt{2}\sqrt{c1/c3}$. The
corresponding fluctuations in potential and kinetic energies can be written

$\displaystyle\sigma_{p}$ $\displaystyle=$
$\displaystyle\frac{\kappa}{2}\sqrt{\langle
r^{n}r^{n}r^{n}r^{n}\rangle-\langle r^{n}r^{n}\rangle^{2}}$ (100)
$\displaystyle=$
$\displaystyle\frac{k_{B}T}{\sqrt{2}}\,\sqrt{\frac{1}{4}\frac{1+7c_{2}^{2}}{1+c_{2}^{2}}+\frac{1}{4}\Gamma_{4}\frac{1-c_{2}^{2}}{1+c_{2}^{2}}}$
$\displaystyle\sigma_{k}$ $\displaystyle=$
$\displaystyle\frac{m}{2}\sqrt{\langle
u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle-\langle
u^{n+\frac{1}{2}}u^{n+\frac{1}{2}}\rangle^{2}}$ (101) $\displaystyle=$
$\displaystyle\frac{k_{B}T}{\sqrt{2}}\,\sqrt{\frac{1}{2}\frac{c_{2}^{3}}{1+c_{2}^{2}}\left[3-\Gamma_{4}\right]+\frac{1}{4}\left[1+\Gamma_{4}\right]}\,.$

Notice that if $\beta^{n}$ is a Gaussian distributed variable, where the
kurtosis is $\Gamma_{4}=3$, the above energy fluctuations yield the expected
values, $\sigma_{p}=\sigma_{k}=k_{B}T/\sqrt{2}$.

Examples of the use of these expressions are shown in Figs. 2cd and 6ef as
open markers ($\circ$) for particular noise distributions.

Special Case, $c_{2}=0$

This special case (GJ-0, see Refs. GJ ) corresponds to the first order
difference equation

$\displaystyle r^{n+1}$ $\displaystyle=$ $\displaystyle
r^{n}+\frac{1}{\alpha}\left[\Delta{t}\,f^{n}+\frac{1}{2}(\beta^{n}+\beta^{n+1})\right]\,,$
(102)

which can approximate the solution to the Langevin equation with no inertia;
i.e., Brownian motion. The linearized, Hooke’s law system ($f=-\kappa r$)
reads

$\displaystyle r^{n+1}$ $\displaystyle=$ $\displaystyle
Xr^{n}+\frac{1}{2\alpha}(\beta^{n}+\beta^{n+1})\,,$ (103)

with $X=1-\kappa\Delta{t}/\alpha$, in correspondence with Eqs. (18) and (19).
This simplicity allows the last four of the ten equations in Eq. (82) to
decouple from the rest, and the fourth moments can be calculated to yield the
average and fluctuations in potential energy

$\displaystyle\left\langle E_{p}(r^{n})\right\rangle$ $\displaystyle=$
$\displaystyle\frac{\kappa}{2}\,\left\langle
r^{n}r^{n}\right\rangle\;=\;\frac{1}{2}k_{B}T$ (104) $\displaystyle\sigma_{p}$
$\displaystyle=$ $\displaystyle\frac{\kappa}{2}\,\sqrt{\left\langle
r^{n}r^{n}r^{n}r^{n}\right\rangle-\left\langle r^{n}r^{n}\right\rangle^{2}}$
(105) $\displaystyle=$
$\displaystyle\frac{k_{B}T}{\sqrt{2}}\,\sqrt{\frac{\frac{1}{2}X^{4}\left[3-\Gamma_{4}\right]+\frac{1}{4}X^{3}\left[7-\Gamma_{4}\right]+\frac{1}{4}(X^{2}+X+1)\left[1+\Gamma_{4}\right]}{X^{3}+X^{2}+X+1}}\,.$

We also here notice that if $\beta^{n}$ is a Gaussian distributed variable,
where the kurtosis is $\Gamma_{4}=3$, the potential energy fluctuations yield
the expected value, $\sigma_{p}=k_{B}T/\sqrt{2}$.

Similarly, we can find the fluctuations $\sigma_{k}$ in kinetic energy from
the correlation

$\displaystyle\left\langle\left(u^{n+\frac{1}{2}}\right)^{4}\right\rangle\;=\;(1-X)^{2}\frac{4\sigma_{p}^{2}+(k_{B}T)^{2}}{m^{2}}+6\left(\frac{k_{B}T}{m}\right)^{2}X^{2}(1-X)-\left(\frac{k_{B}T}{m}\right)^{2}(3+\Gamma_{4})X\left(1-\frac{3}{2}X+X^{2}\right)\,,$
(106)

such that

$\displaystyle\sigma_{k}$ $\displaystyle=$
$\displaystyle\frac{1}{2}m\,\sqrt{\left\langle\left(u^{n+\frac{1}{2}}\right)^{4}\right\rangle-\left(\frac{k_{B}T}{m}\right)^{2}}\,,$
(107)

where we have used Eq. (15). Thus, $\sigma_{k}$ can be easily computed also in
this special case. We again notice that for Gaussian fluctuations, where
$\Gamma_{4}=3$, the kinetic energy fluctuations yield the expected value
$\sigma_{k}=k_{B}T/\sqrt{2}$.

Examples of the use of these expressions are shown in Figs. 2cd and 6ef as
closed markers ($\bullet$) for particular noise distributions.

Figure 1: The tree applied distributions, $\rho_{G}(\beta^{n})$ (Gaussian),
$\rho_{U}(\beta^{n})$ (Uniform), and $\rho_{P}(\beta^{n})$ (Peaked), given in
Eqs. (Gaussian,), (Uniform,), and (Peaked,), respectively. All three
distributions satisfy the discrete-time fluctuation-dissipation relationship,
Eq. (10).

Figure 2: Simulations and analyses of average potential and kinetic energies,
$\langle E_{p}\rangle$ from Eq. (31) (a) and $\langle E_{k}\rangle$ from Eq.
(32) (b), and their respective fluctuations, $\sigma_{p}$ from Eq. (96) (c)
and $\sigma_{k}$ from Eq. (LABEL:eq:sigma_k) (d), as a function of reduced
time step for the harmonic system, Eqs. (6) and (17), with non-zero
temperature $T$. Results are shown for Gaussian $\beta_{G}^{n}$ (Eq. (108a))
and uniform $\beta_{U}^{n}$ (Eq. (LABEL:eq:Flatn)) noise distributions.
Several friction parameters $\alpha$ are used, as indicated in the figure.
Simulation results are indistinguishable from the analytical ones. Markers
indicate special cases for $c_{2}=0$ ($\bullet$ from Eqs. (105)-(107)) and
$X=0$ ($\circ$ from Eqs. (100) and (101)).

Figure 3: Configurational and kinetic distributions $\rho_{c}(r^{n})$ and
$\rho_{k}(u^{n+\frac{1}{2}})$ for simulations with Gaussian and uniform noise
applied to the linear system Eqs. (6) and (17) with $\alpha=1m\Omega_{0}$,
$\Omega_{0}\Delta{t}=0.5$, $k_{B}T=E_{0}$. Shown curves are Eq. (3) (solid,
labeled exact), simulations data using Gaussian noise $\beta_{G}^{n}$ (dashed,
labeled Gaussian), and simulation data using uniform noise $\beta_{U}^{n}$
(dash-dot, labeled uniform). Distributions from simulations using Gaussian
noise are indistinguishable from the exact references Eq. (3). (c) and (d)
show the distributions, (a) and (b) show the deviations from the exact
distributions Eq. (3), and (e) and (f) show the effective potential (PMF) and
kinetic (MB) potentials (Eq. (113)) derived from the distributions.

Figure 4: Configurational and kinetic distributions $\rho_{c}(r^{n})$ and
$\rho_{k}(u^{n+\frac{1}{2}})$ for simulations with Gaussian and uniform noise
applied to the linear system Eqs. (6) and (17) with $\alpha=1m\Omega_{0}$,
$\Omega_{0}\Delta{t}=1.5$, $k_{B}T=E_{0}$. Shown curves are Eq. (3) (solid,
labeled exact), simulations data using Gaussian noise $\beta_{G}^{n}$ (dashed,
labeled Gaussian), and simulation data using uniform noise $\beta_{U}^{n}$
(dash-dot, labeled uniform). Distributions from simulations using Gaussian
noise are indistinguishable from the exact references Eq. (3). (c) and (d)
show the distributions, (a) and (b) show the deviations from the exact
distributions Eq. (3), and (e) and (f) show the effective potential (PMF) and
kinetic (MB) potentials (Eq. (113)) derived from the distributions.

Figure 5: Configurational and kinetic distributions $\rho_{c}(r^{n})$ and
$\rho_{k}(u^{n+\frac{1}{2}})$ for simulations with Gaussian and uniform noise
applied to the linear system Eqs. (6) and (17) with $\alpha=1m\Omega_{0}$,
$\Omega_{0}\Delta{t}=1.975$, $k_{B}T=E_{0}$. Shown curves are Eq. (3) (solid,
labeled exact), simulations data using Gaussian noise $\beta_{G}^{n}$ (dashed,
labeled Gaussian), and simulation data using uniform noise $\beta_{U}^{n}$
(dash-dot, labeled uniform). Distributions from simulations using Gaussian
noise are indistinguishable from the exact references Eq. (3). (c) and (d)
show the distributions, (a) and (b) show the deviations from the exact
distributions Eq. (3), and (e) and (f) show the effective potential (PMF) and
kinetic (MB) potentials (Eq. (113)) derived from the distributions.

Figure 6: Simulations and analyses of average potential and kinetic energies,
$\langle E_{p}\rangle$ from Eq. (31) (a) and $\langle E_{k}\rangle$ from Eq.
(32) (b), the respective third moments $\langle(r^{n})^{3}\rangle$ from Eqs.
(45)-(53) (c) and $\langle(u^{n+\frac{1}{2}})^{3}\rangle$ from Eq. (55) (d),
and the respective energy fluctuations, $\sigma_{p}$ from Eq. (96) (e) and
$\sigma_{k}$ from Eq. (LABEL:eq:sigma_k) (f), as a function of reduced time
step for the harmonic system, Eqs. (6) and (17), with non-zero temperature
$T$. Results are shown for Gaussian $\beta_{G}^{n}$ (Eq. (108a)) and the
asymmetrically peaked $\beta_{P}^{n}$ (Eq. (110a)) noise distributions.
Several friction parameters $\alpha$ are used, as indicated in the figure.
Simulation results are indistinguishable from the analytical ones. Markers
indicate special cases for $c_{2}=0$ ($\bullet$ from Eqs. (59a) and (61), Eqs.
(105)-(107)) and $X=0$ ($\circ$ from Eqs. (56a) and (58), Eqs. (100) and
(101)).

Figure 7: Configurational and kinetic distributions $\rho_{c}(r^{n})$ and
$\rho_{k}(u^{n+\frac{1}{2}})$ for simulations with Gaussian and peaked noise
applied to the linear system Eqs. (6) and (17) with $\alpha=0.5m\Omega_{0}$,
$\Omega_{0}\Delta{t}=0.5$, $k_{B}T=E_{0}$. Shown curves are Eq. (3) (solid,
labeled exact), simulations data using Gaussian noise $\beta_{G}^{n}$ (dashed,
labeled Gaussian), and simulation data using peaked noise $\beta_{P}^{n}$
(dash-dot, labeled uniform). Distributions from simulations using Gaussian
noise are indistinguishable from the exact references Eq. (3). (c) and (d)
show the distributions, (a) and (b) show the deviations from the exact
distributions Eq. (3), and (e) and (f) show the effective potential (PMF) and
kinetic (MB) potentials (Eq. (113)) derived from the distributions.

Figure 8: Configurational and kinetic distributions $\rho_{c}(r^{n})$ and
$\rho_{k}(u^{n+\frac{1}{2}})$ for simulations with Gaussian and uniform noise
applied to the linear system Eqs. (6) and (17) with $\alpha=0.5m\Omega_{0}$,
$\Omega_{0}\Delta{t}=1.5$, $k_{B}T=E_{0}$. Shown curves are Eq. (3) (solid,
labeled exact), simulations data using Gaussian noise $\beta_{G}^{n}$ (dashed,
labeled Gaussian), and simulation data using peaked noise $\beta_{P}^{n}$
(dash-dot, labeled uniform). Distributions from simulations using Gaussian
noise are indistinguishable from the exact references Eq. (3). (c) and (d)
show the distributions, (a) and (b) show the deviations from the exact
distributions Eq. (3), and (e) and (f) show the effective potential (PMF) and
kinetic (MB) potentials (Eq. (113)) derived from the distributions.

Figure 9: Configurational and kinetic distributions $\rho_{c}(r^{n})$ and
$\rho_{k}(u^{n+\frac{1}{2}})$ for simulations with Gaussian and uniform noise
applied to the linear system Eqs. (6) and (17) with $\alpha=0.5m\Omega_{0}$,
$\Omega_{0}\Delta{t}=1.975$, $k_{B}T=E_{0}$. Shown curves are Eq. (3) (solid,
labeled exact), simulations data using Gaussian noise $\beta_{G}^{n}$ (dashed,
labeled Gaussian), and simulation data using peaked noise $\beta_{P}^{n}$
(dash-dot, labeled uniform). Distributions from simulations using Gaussian
noise are indistinguishable from the exact references Eq. (3). (c) and (d)
show the distributions, (a) and (b) show the deviations from the exact
distributions Eq. (3), and (e) and (f) show the effective potential (PMF) and
kinetic (MB) potentials (Eq. (113)) derived from the distributions.

##  IV Simulations

While the above expressions provide both direct and indirect means for
evaluating the thermodynamically important first four moments of the
configurational and kinetic variables, we test these results against numerical
simulations such that both simulations and analysis can be mutually verified.
Further, the simulations provide additional details on the actual
distributions $\rho_{c}(r^{n})$ and $\rho_{k}(u^{n+\frac{1}{2}})$, which in
some cases reveal significant deviations from Boltzmann statistics when the
fluctuations $\beta^{n}$ are non-Gaussian.

We show results of three different distributions for the fluctuations
$\beta^{n}$, all satisfying the discrete-time fluctuation-dissipation balance
that constrains the first two moments given in Eq. (10), and thus, all
resulting in correct first and second moments of $r^{n}$ and
$u^{n+\frac{1}{2}}$, leading to correct temperature as measured by the kinetic
and potential energies (see Sec. III.1):

Gaussian,

$\beta^{n}=\beta_{G}^{n}$ with distribution $\rho_{G}(\beta)$:

$\displaystyle\rho_{G}(\beta)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi\langle\beta^{n}\beta^{n}\rangle}}\,\exp\left(-\frac{\beta\beta}{2\langle\beta^{n}\beta^{n}\rangle}\right)\,,$
(108a)

with skewness and kurtosis

$\displaystyle\Gamma_{3}$ $\displaystyle=$ $\displaystyle 0$ (108b)
$\displaystyle\Gamma_{4}$ $\displaystyle=$ $\displaystyle 3\,.$ (108c)

Uniform,

$\beta^{n}=\beta_{U}^{n}$ with distribution $\rho_{U}(\beta)$:

$\displaystyle\rho_{U}(\beta)$ $\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{3\langle\beta^{n}\beta^{n}\rangle}}\times\left\\{\begin{array}[]{ccc}1&,&-\sqrt{3}\,\leq\,\frac{\beta}{\sqrt{\langle\beta^{n}\beta^{n}\rangle}}\,<\,\sqrt{3}\\\
0&,&{\rm otherwise}\end{array}\right.,$ (109c)

with skewness and kurtosis

$\displaystyle\Gamma_{3}$ $\displaystyle=$ $\displaystyle 0$ (109e)
$\displaystyle\Gamma_{4}$ $\displaystyle=$ $\displaystyle\frac{9}{5}\,.$
(109f)

Peaked,

$\beta^{n}=\beta_{P}^{n}$ with distribution $\rho_{P}(\beta)$:

$\displaystyle\rho_{P}(\beta)$ $\displaystyle=$
$\displaystyle\frac{2}{3}\,\delta\left(\frac{\beta}{\sqrt{\langle\beta^{n}\beta^{n}\rangle}}+\frac{1}{\sqrt{2}}\right)$
(110a) $\displaystyle+$
$\displaystyle\frac{1}{3}\,\delta\left(\frac{\beta}{\sqrt{\langle\beta^{n}\beta^{n}\rangle}}-\sqrt{2}\right)\,,$

with skewness and kurtosis

$\displaystyle\Gamma_{3}$ $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2}}$
(110b) $\displaystyle\Gamma_{4}$ $\displaystyle=$
$\displaystyle\frac{12}{5}\,.$ (110c)

The three distributions are sketched in Fig. 1, and the realizations of the
pseudo-random numbers $\beta^{n}$ are produced from the uniformly distributed
numbers in $[0,1)$ generated by the RANMAR algorithm (used in LAMMPS LAMMPS-
Manual ).

In order to validate the analysis of the previous section, we simulated the
noisy harmonic oscillator Eqs. (6) and (17) with several of the GJ methods,
which all give the same statistical results for Gaussian noise once the time
step is appropriately scaled according to the method. For simplicity we limit
the displayed results to those of the GJ-I (GJF-2GJ) method 2GJ ; GJ , for
which

$\displaystyle c_{2}$ $\displaystyle=$
$\displaystyle\frac{1-\frac{\alpha\Delta{t}}{2m}}{1+\frac{\alpha\Delta{t}}{2m}}$
(111a) $\displaystyle c_{1}\;=\;c_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\frac{\alpha\Delta{t}}{2m}}\,.$ (111b)

The GJF-2GJ method overlaps in the configurational coordinate with the GJF
method GJF1 . We note that the GJ-I (GJF-2GJ) method is available for use in
LAMMPS LAMMPS-Manual . With time normalized by the characteristic unit $t_{0}$
given by the inverse of the natural frequency $\Omega_{0}$, the reduced
damping is $\alpha/m\Omega_{0}$, and the characteristic energy is
$E_{0}=\kappa r_{0}^{2}$, where $r_{0}$ is a chosen characteristic length to
which $r$ is normalized. The simulations are conducted for $k_{B}T=E_{0}$ for
a variety of values of $\alpha/m\Omega_{0}$ in the entire stability range
$\Omega_{0}\Delta{t}<2$. Each statistical value is calculated from averages
over 1,000 independent simulations, each with 108 time steps.

Figure 2 shows the results of using the uniform noise distribution
$\beta_{U}^{n}$ (Eq. (LABEL:eq:Flatn)) alongside the results of using the
Gaussian $\beta_{G}^{n}$ (Eq. (108a)) as reference. Figures 2ab, respectively
displaying results for configurational and kinetic variables, verify the
general result from Sec. III.1; namely that a statistical property, which
depends only on the first and second moments of the noise distribution, will
be correctly evaluated for the GJ methods if the noise satisfies Eq. (10). The
shown data on each plot are for damping parameters $\alpha/m\Omega_{0}$=0.1,
0.25, 0.5, 1, 2.5, 5, and the results for both Gaussian and uniform noise
distributions are clearly in close agreement with the expected values given in
Eqs. (31) and (32). This has previously been extensively validated for
Gaussian noise 2GJ ; GJ .

Figures 2cd display simulation results and comparisons to the expectations
from the analysis in Sec. III.3 for the configurational and kinetic energy
fluctuations $\sigma_{p}$ and $\sigma_{k}$ found in Eqs. (96) and
(LABEL:eq:sigma_k). Since the Gaussian noise $\beta_{G}^{n}$ from Eq. (108a)
applied to the linear system produces Gaussian distributions for $r^{n}$ and
$u^{n+\frac{1}{2}}$ with correct, and time step independent, first and second
moments (see Sec. III.1), it follows that the fourth moments are also correct
and time step independent. This is clearly observed on the energy fluctuation
figures, where both simulation results and the results of solving Eqs.
(82)-(94) for $\Gamma_{4}=3$ are shown to be indistinguishable and constant at
the correct values $\sigma_{p}=\sigma_{k}=k_{B}T/\sqrt{2}$. However, the
uniformly distributed noise $\beta_{U}^{n}$ (Eq. (LABEL:eq:Flatn)) produces
neither Gaussian nor uniform distributions for $r^{n}$ and
$u^{n+\frac{1}{2}}$, but instead produces distributions that depend on the
time step and the damping parameter, even if Eq. (10) is satisfied. This can
be seen in the figures for the energy fluctuations, Figs. 2cd, where both
simulation results and the results of Eqs. (82)-(94) for $\Gamma_{4}=3$ and
$\Gamma_{4}=\frac{9}{5}$, the latter being the kurtosis for the uniform
distribution, are shown. The simulation results are indistinguishable from the
analysis. It is obvious that when using the uniform noise, the energy
fluctuations, which depend on the fourth moments of the variables, can deviate
significantly from the correct value, given by the Gaussian noise. We also
observe that the fluctuations approach the correct values for small time steps
or small damping parameters. This is understandable, since the damping per
time step in those limits is very small, and the fluctuations in $r^{n}$ and
$u^{n+\frac{1}{2}}$ therefore are composed of $\beta_{U}^{n}$ contributions
from many time steps, allowing the integrated noise in $r^{n}$ and
$u^{n+\frac{1}{2}}$ to become near-Gaussian by the central limit theorem.
Conversely, if the damping per time step becomes appreciable, then the
effective noise in $r^{n}$ and $u^{n+\frac{1}{2}}$ will be composed by only a
few $\beta_{U}^{n}$ contributions, which do not approximate a Gaussian outcome
very well.

In order to look at the details of the simulated distributions, we have
selected a few representative parameter values as illustrations. Figures 3-5
show the simulated configurational and kinetic distributions, $\rho_{c}(r)$
and $\rho_{k}(u)$, for $\alpha=1m\Omega_{0}$ and
$\Omega_{0}\Delta{t}=0.5,1.5,1.975$, respectively, for simulations using
Gaussian (dashed) and uniform (dash-dotted) noise distributions. Also shown
are the results from the exact Gaussian distributions, $\rho_{k,e}(u)$ and
$\rho_{c,e}(r)$, (solid, from Eqs. (3a) and (12a)) that are expected from
continuous-time Langevin dynamics. The distributions $\rho_{c}(r^{n})$ and
$\rho_{k}(u^{n+\frac{1}{2}})$, are shown in Figures 3-5cd, the deviations

$\displaystyle\Delta\rho_{c}(r^{n})$ $\displaystyle=$
$\displaystyle\rho_{c}(r^{n})-\rho_{c,e}(r^{n})$ (112a)
$\displaystyle\Delta\rho_{k}(u^{n+\frac{1}{2}})$ $\displaystyle=$
$\displaystyle\rho_{k}(u^{n+\frac{1}{2}})-\rho_{k,e}(u^{n+\frac{1}{2}})\,,$
(112b)

from the expected (correct) distributions are shown in Figures 3-5ab, and the
effective configurational and kinetic potentials,

$\displaystyle U_{PMF}(r^{n})$ $\displaystyle=$ $\displaystyle-
k_{B}T\ln\rho_{c}(r^{n})+\tilde{{\cal C}}_{c}$ (113a) $\displaystyle
U_{MB}(u^{n+\frac{1}{2}})$ $\displaystyle=$ $\displaystyle-
k_{B}T\ln\rho_{k}(u^{n+\frac{1}{2}})+\tilde{{\cal C}}_{k}\,,$ (113b)

where $\tilde{{\cal C}}_{c}$ and $\tilde{{\cal C}}_{k}$ are determined such
that ${\rm min}[U_{PMF}(r^{n})]={\rm min}[U_{MB}(u^{n+\frac{1}{2}})]=0$, are
shown in Figures 3-5ef. For the harmonic potential, the statistically correct
values of these potentials should be $\frac{1}{2}\kappa(r)^{2}$ and
$\frac{1}{2}m(u)^{2}$, which are indicated by thin solid curves. For
$\Omega_{0}\Delta{t}=0.5$, Figure 3 shows that even a seemingly modest
deviation from a pure Gaussian distribution can have rather large impacts on
fourth-moment thermodynamic measures, which for these parameters in Fig. 2cd
are seen to result in configurational and kinetic energy fluctuations being
depressed by about 7% and 14%, respectively. Increasing the time step to
$\Omega_{0}\Delta{t}=1.5$ amplifies the deformation of the
$\beta_{U}^{n}$-generated distributions away from the Gaussian. As seen in
Figs. 4ef, the increasingly noticeable difference is that the uniform noise
yields a more confined exploration of phase-space than the Gaussian
distribution does, consistent with the depression in the fourth moment seen in
Figs. 2cd. Finally, in Fig. 5, we show the results for a time step,
$\Omega_{0}\Delta{t}=1.975$, very close to the stability limit. This extreme
case shows that the uniform distribution in noise also can produce a near
uniform distribution for $r^{n}$, while the distribution for
$u^{n+\frac{1}{2}}$ is near triangular, consistent with a sum of two uniformly
distributed numbers contributing to the kinetic fluctuations. Again, we see
from Figs. 5ef that the sampling of the phase space is much more limited when
using uniform noise than when using Gaussian, even if the measured
temperature, configurational as well as kinetic, are measured to the correct
values.

Following the spirit of Refs. Greiner_1 ; Dunweg_1  we explore the application
of a more challenging noise distribution $\beta_{P}^{n}$ (the peaked
distribution defined in Eq. (110a)), which is both discrete and asymmetric
(see Fig. 1). As was the case for the uniform noise discussed above, we also
here reference the results next to those of Gaussian noise, which give time
step independent and correct behavior. Conducting numerical simulations of
Eqs. (6) and (17) for different time steps and for normalized friction
$\alpha/m\Omega_{0}=0.1,0.25,0.5,1,2.5,5$ as described above, and comparing to
the evaluation of the moments from the analyses of Sec. III, we obtain the
data shown in Fig. 6, where the simulation results are indistinguishable from
the results of the analyses for third and fourth moments in Secs. III.2 and
III.3. As expected from the analysis, the results of the second moments shown
in Fig. 6ab of $r^{n}$ and $u^{n+\frac{1}{2}}$, namely the potential and
kinetic energies, are perfectly aligned with statistical mechanics, since the
peaked distribution for $\beta_{P}^{n}$ satisfies the two moments in Eq. (10).
Since the applied noise is here asymmetric, the resulting third moments of
$r^{n}$ and $u^{n+\frac{1}{2}}$ may also become non-zero for non-zero time
steps. This is seen in Figures 6cd, where we observe rather complex behavior
as a function of the system parameters. We do, however, see that for
$\Omega_{0}\Delta{t}\rightarrow 0$ the third moments approach zero, consistent
with the central limit theorem that guarantees a Gaussian outcome when a very
large number of noise values contribute to the variables. Finally, in
comparison with Fig. 2cd for uniform noise, we see very similar behavior for
the energy fluctuations (fourth moments) for the application of the peaked
distribution in Fig. 6ef. Even if the uniform and peaked distributions are
very different in appearance, the similarities between their outcomes in their
fourth moments are not surprising, given that the analysis for the fourth
moments in Sec. III.3 shows that the difference between the two only depends
on the kurtoses, which have the values $\Gamma_{4}=\frac{9}{5}$ (uniform) and
$\Gamma_{4}=\frac{12}{5}$ (peaked). Again, we find that all signatures of non-
Gaussian noise vanish for $\Omega_{0}\Delta{t}\rightarrow 0$, and we find that
non-zero time steps result in depressions of the thermodynamically important
energy fluctuations, as derived from the fourth moments of the configurational
and kinetic coordinates.

The details of the simulated coordinate distributions arising from the peaked
noise are exemplified in Figs. 7-9 for $\alpha/m\Omega_{0}=0.5$ and
$\Omega_{0}\Delta{t}=0.5,1.5,1.975$, respectively. For the smaller of the time
steps, shown in Fig. 7, we see the seemingly modest skewness and deformations
of the coordinate distributions arising from the peaked noise. Yet, these
modest deformations are what provide the somewhat significant deviations from
Gaussian characteristics in third and fourth moments seen in Figs. 6cd and
Figs. 6ef for $\Omega_{0}\Delta{t}=0.5$. More dramatic deviations from
Gaussian/Boltzmann characteristics are found in a large range of
$\Omega_{0}\Delta{t}$, including the value $\Omega_{0}\Delta{t}=1.5$ shown in
Fig. 8. The seemingly discontinuous distribution is not a result of
insufficient statistics. Rather, it is the interference between the discrete
nature of the applied noise distribution with the discrete time step that
happens to distinctively select certain preferred values of $r^{n}$ and
$u^{n+\frac{1}{2}}$ over others. It is noticeable that these types of peculiar
distributions appear without any obvious or abrupt signatures in the first
four moments, as seen in Fig. 6. As the time step $\Omega_{0}\Delta{t}=1.975$
is pushed close to the stability limit, we again find smooth coordinate
distributions, seen in Fig. 9, where the applied peaked noise is visible
throughout the different displays. In analogy with the visualization of the
distributions from the uniform noise, this is the limit where the friction per
time step is relatively large, thereby making the behavior of the coordinates
$r^{n}$ and $u^{n+\frac{1}{2}}$ subject to only a few noise contributions at a
time.

##  V Discussion

In light of recent advances in stochastic thermostats for simulating Langevin
equations with accurate statistics across the stability range of the applied
time step when using Gaussian noise 2GJ ; GJ , we have revisited the
investigations of advantages and disadvantages of using non-Gaussian thermal
noise in discrete time. Given the systemic first and second order time step
errors of the traditional methods (e.g., Refs. BBK ; SS ), which necessitated
rather small time steps for accurate simulations, it was previously concluded
Greiner_1 ; Dunweg_1  that other distributions, including the desirable
uniform distribution, would be efficient substitutions since the central limit
theorem would ensure near-Gaussian outcomes for small time steps, thereby not
further significantly distort the simulation results due to the imperfect
noise. This result has been a very useful tool over the years when
computational efficiency could benefit from not converting stochastic
variables from uniform to Gaussian. However, as we have demonstrated in this
paper, when the time step becomes large, which is allowed by the modern GJ
methods when using Gaussian noise, the application of non-Gaussian noise does
not retain the time step independent benefits of these methods in
thermodynamic measures that involve moments higher than the second. We have
found that, given that the applied noise conforms to Eq. (10), the GJ methods
are invariant to the specific noise distributions in measures of first and
second moments, such as configurational and kinetic temperatures. Thus, it is
deceiving to judge the quality of the thermostat based on those moments alone.
As is evident from the third and fourth moments, as well as the visual
impressions of the actual distributions of $r^{n}$ and $u^{n+\frac{1}{2}}$,
the sampling of the phase-space can be quite distorted (non-Boltzmann) even if
the measured temperatures yield correct values. All of these results are
consistent with the significance of how the noise must be defined in discrete-
time; namely through the integral Eq. (9), which ensures that any underlying
distribution for $\beta(t)$ will yield a Gaussian outcome for $\beta^{n}$ due
to the central limit theorem. It follows that applying any other distribution
than Gaussian in discrete time is formally invalid, and certainly leads to
significant sampling errors unless $\alpha\Delta{t}/m$ is small enough that
the one-time-step attenuation parameter $c_{2}\approx 1$. It is therefore the
conclusion of this work that Gaussian noise must be applied to the modern
stochastic integrators if one wishes to take advantage of the large time step
benefits of their properties.

##  VI Acknowledgments

The author is grateful to Charlie Sievers for sharing LAMMPS simulations that
support the results of this paper in more complex MD simulations using the
GJF-2GJ method, and to Lorenzo Mambretti for assistance with the RANMAR random
number generator. The author is also grateful for initial discussions with
Chungho Cheng on non-Gaussian noise.

## References

  * (1) M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Inc., 1989.
  * (2) D. Frenkel and B. Smit, Understanding Molecular Simulations: From Algorithms to Applications, (Academic Press, San Diego, 2002).
  * (3) D. C. Rapaport, The Art of Molecular Dynamics Simulations, (Cambridge University Press, Cambridge, 2004).
  * (4) W. M. Hoover, Computational Statistical Mechanics, (Elsevier Science B.V., 1991).
  * (5) A. Leach, Molecular Modeling: Principles and Applications, 2nd ed., (Prentice Hall: Harlow, England, 2001).
  * (6) N. Grønbech-Jensen and O. Farago, Mol. Phys. 111, 983 (2013).
  * (7) L. F. G. Jensen and N. Grønbech-Jensen, Mol. Phys. 117, 2511 (2019).
  * (8) N. Grønbech-Jensen, Mol. Phys. 118, e1662506 (2020).
  * (9) P. Langevin, C. R. Acad. Sci. Paris 146, 530 (1908).
  * (10) W. T. Coffey and Y. P. Kalmykov, The Langevin Equation, 3rd ed., World Scientific, Singapore, 2012.
  * (11) G. Parisi, Statistical Field Theory, (Addison-Wesley, Menlo Park, 1988).
  * (12) Although not necessary, it is generally reasonable to expect that the distribution of $\beta(t)$ is symmetric, yielding all odd moments zero.
  * (13) A. Papoulis, Probability, Random Variables, and Stochastic Processes, (McGraw-Hill, London, 1965).
  * (14) For a general description of pseudo-random generators, see, e.g., W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes 3rd. ed. (Cambridge University Press, Cambridge, 2007).
  * (15) A. Greiner, W. Strittmatter, and J. Honerkamp, J. Stat. Phys. 51, 95 (1988).
  * (16) B. Dünweg and W. Paul, Int. J. Mod. Phys. C 2, 817 (1991).
  * (17) J. Finkelstein, G. Fiorin, and S. Seibold, Mol. Phys. 118, e1649493 (2020).
  * (18) A. Brünger, C. L. Brooks, and M. Karplus, Chem. Phys. Lett. 105, 495 (1984).
  * (19) R. W. Pastor, B. R. Brooks, and A. Szabo, Mol. Phys. 65, 1409 (1988).
  * (20) T. Schneider and E. Stoll, Phys. Rev. B 17, 1302 (1978).
  * (21) S. Plimpton, J. Comp. Phys. 117, 1 (1995).
  * (22) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ′t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. Stevens, J. Tranchida, C. Trott, S. J. Plimpton, Compute. Phys. Commun. 271, 108171 (2022).
  * (23) See the LAMMPS documentation LAMMPS-Manual  for the use of fix_langevin using the Schneider and Stoll thermostat of Ref. SS .
  * (24) F. Weik, R. Weeber, K. Szuttor, K Breitsprecher, J. de Graaf, M. Kuron, J. Landsgesell, H. Menke, D. Sean and C. Holm., Euro. Phys. J. Special Topics 227, 1789 (2019). See also https://espressomd.github.io/ documentation for ESPResSo-4.2.0 Sec. 6.3.1.
  * (25) L. Verlet, Phys. Rev. 159, 98 (1967).
  * (26) W. C. Swope, H. C. Andersen, P. H. Berens, K. R. Wilson, J. Chem. Phys. 76, 637 (1982).
  * (27) D. Beeman, J. Comp. Phys. 20, 130 (1976).
  * (28) O. Bumeman, J. Comp. Phys. 1, 517 (1967).
  * (29) R. W. Hockney, Methods Comput. Phys. 9, 136 (1970).
  * (30) See http://lammps.sandia.gov/doc/Manual.pdf, for the description of the “fix_langevin” command.