# Accurate and robust splitting methods for the generalized Langevin equation
with a positive Prony series memory kernel

Manh Hong Duong  School of Mathematics, University of Birmingham, Edgbaston,
Birmingham, B15 2TT, United Kingdom<EMAIL_ADDRESS>Xiaocheng Shang  The
Alan Turing Institute, British Library, 96 Euston Road, London, NW1 2DB,
United Kingdom<EMAIL_ADDRESS>
###### Abstract

We study numerical methods for the generalized Langevin equation (GLE) with a
positive Prony series memory kernel, in which case the GLE can be written in
an extended variable Markovian formalism. We propose a new splitting method
that is easy to implement and is able to substantially improve the accuracy
and robustness of GLE simulations in a wide range of the parameters. An error
analysis is performed in the case of a one-dimensional harmonic oscillator,
revealing that several averages are exact for the newly proposed method.
Various numerical experiments in both equilibrium and nonequilibrium
simulations are also conducted to compare the method with popular alternatives
in interacting multi-particle systems.

###### keywords:

Stochastic differential equations , Splitting methods , Generalized Langevin
equation , Memory kernel , Error analysis , Harmonic oscillator

##  1 Introduction

It is well known that the popular Langevin dynamics is based on the assumption
that there is a clear separation between the characteristic time scale of the
“massive” particles and that of those smaller solvent particles forming the
heat bath. Within this assumption, the effect of the solvent is simply reduced
to an instantaneous drag force and a delta-correlated random force, thereby
significantly reducing the computational cost of an “explicit solvent” model.
However, it has been well documented that this underlying assumption breaks
down in many physically compelling scenarios, in which case a generalized
version of the Langevin dynamics is more suitable—it is known as the
generalized Langevin equation (GLE) derived from the Mori–Zwanzig formalism
[42, 55]. The GLE has been widely applied in a large number of applications,
including molecular simulations [21, 30, 17], mesoscopic modeling [11, 22],
solids [23, 24, 52], nuclear quantum effects [8, 9], and various anomalous
diffusion phenomena [16, 40].

Unlike the Langevin dynamics, a temporally nonlocal drag force and a random
force with nontrivial correlations are typically associated with the GLE,
which make its numerical integration highly nontrivial [21, 6]. To be more
precise, the temporally nonlocal drag force, in the form of a convolution of
the momentum with a memory kernel, requires the storage of the momentum
history, whose numerical evaluation at each time step can be computationally
demanding. On the other hand, it can also be computationally expensive to
generate a random force with nontrivial correlations, which may require the
storage of a sequence of random numbers at each time step. Although
considerable effort has been devoted to developing accurate and robust
numerical methods that circumvent either or both of these challenges mentioned
above, it is very difficult to claim any individual method as being optimal,
especially given the fact that there are such a broad range of applications
where the GLE can be applied.

To this end, in this article we focus our attention on a general form of the
memory kernel, which is a sum of exponentials and also known as a positive
Prony series. It has been widely used in a large number of studies in the
literature [15, 41, 45, 6, 46, 20]. Moreover, as pointed out in [1, 41, 16],
with certain choice of the parameters, a Prony series can be seen as an
approximation of a power law [28, 40, 16], another primary example of the
memory kernel.

The rest of the article is organized as follows. In Section 2, we discuss the
mathematical formulations of the GLE, including not only the extended variable
Markovian formalism and its white noise limit but also the Markovian
approximations of the GLE with a more general memory kernel. Section 3 reviews
three popular integration methods for GLE, followed by the derivations of a
new promising scheme. Error analysis on the averages in the case of a one-
dimensional harmonic oscillator will also be performed. A variety of numerical
experiments are performed in Section 4 to compare all the schemes described in
the article. Our findings are summarized in Section 5.

##  2 Mathematical formulations

Consider an $N$-particle system evolving in dimension $d$ with position
${{\mathbf{q}}}_{i}\in\mathbb{R}^{d}$, momentum
${{\mathbf{p}}}_{i}\in\mathbb{R}^{d}$, and mass $m_{i}\in\mathbb{R}$ for
$i=1,\dots,N$, the equations of motion for the non-Markovian form of the GLE
are given by

$\displaystyle{\rm d}{{\mathbf{q}}}_{i}$
$\displaystyle=m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm d}t\,,$ (1a)
$\displaystyle{\rm d}{{\mathbf{p}}}_{i}$
$\displaystyle=-\nabla_{{{\mathbf{q}}}_{i}}U({{\mathbf{q}}}){\rm
d}t-\int^{t}_{0}\hat{K}(t-s)m_{i}^{-1}{{\mathbf{p}}}_{i}(s)\,{\rm d}s\,{\rm
d}t+{\boldsymbol{\eta}}_{i}{\rm d}t\,,$ (1b)

where $U({{\mathbf{q}}})$ is a smooth potential energy. Component-wise, the
random force $\eta^{x}_{i}(t)$, where the index $x$ indicates the appropriate
Cartesian components, is a mean zero stationary Gaussian process with an
autocorrelation function $\hat{K}$, satisfying the fluctuation-dissipation
relation [7, 27],

$\langle\eta^{x}_{i}(t+s)\eta^{y}_{j}(t)\rangle=\beta^{-1}\hat{K}(s)\delta_{ij}\delta_{xy}\,,\quad
s\geq 0\,,$ (2)

where $\beta$ denotes the inverse temperature, both $\delta_{ij}$ and
$\delta_{xy}$ are Kronecker delta functions. In this article, we focus our
attention on the following general form of the memory kernel:

$\hat{K}(t)=\sum^{M}_{k=1}\tilde{\lambda}^{2}_{k}\exp\left(-\tilde{\alpha}_{k}t\right)\,,\quad
t\geq 0\,,$ (3)

where $M$ is a positive integer that represent the number of modes,

$\tilde{\lambda}_{k}=\lambda_{k}/\sqrt{\epsilon}>0\,,\quad\tilde{\alpha}_{k}=\alpha_{k}/\epsilon>0\,,$
(4)

with $\lambda_{k}$ and $\alpha_{k}$ being two constant parameters and
$\epsilon>0$ being a rescaling parameter (see more discussions in [45]). Note
that for the sake of notational simplicity we do not explicitly express the
dependence of $\tilde{\lambda}_{k}$ and $\tilde{\alpha}_{k}$ on $\epsilon$.

###  2.1 Extended variable Markovian formalism

In order to avoid dealing with the integral form of (1b), it is desirable to
rewrite (1) as an extended variable Markovian formalism [28, 9, 45, 6, 46, 10,
36]. More specifically, following [6], we define the extended variable,
${\boldsymbol{\xi}}_{i,k}\in\mathbb{R}^{d}$, associated with the $k$-th Prony
mode’s action on the $i$-th component of ${{\mathbf{q}}}$ and
${{\mathbf{p}}}$:

${\boldsymbol{\xi}}_{i,k}=-\int^{t}_{0}\tilde{\lambda}_{k}\exp\left[-\tilde{\alpha}_{k}(t-s)\right]m_{i}^{-1}{{\mathbf{p}}}_{i}(s)\,{\rm
d}s\,.$ (5)

Subsequently, (1) may be rewritten as

$\displaystyle{\rm d}{{\mathbf{q}}}_{i}$
$\displaystyle=m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm d}t\,,$ (6a)
$\displaystyle{\rm d}{{\mathbf{p}}}_{i}$
$\displaystyle=-\nabla_{{{\mathbf{q}}}_{i}}U({{\mathbf{q}}}){\rm
d}t+\sum^{M}_{k=1}\tilde{\lambda}_{k}{\boldsymbol{\xi}}_{i,k}{\rm
d}t+{\boldsymbol{\eta}}_{i}{\rm d}t\,.$ (6b)

Differentiating (5) gives a simple stochastic differential equation (SDE):

${\rm
d}{\boldsymbol{\xi}}_{i,k}=-\tilde{\alpha}_{k}{\boldsymbol{\xi}}_{i,k}{\rm
d}t-\tilde{\lambda}_{k}m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm d}t\,.$ (7)

In order to construct a random force that satisfies the fluctuation-
dissipation relation (2), we consider the following SDE:

${\rm
d}{\boldsymbol{\eta}}_{i,k}=-\tilde{\alpha}_{k}{\boldsymbol{\eta}}_{i,k}{\rm
d}t+\sqrt{2\tilde{\alpha}_{k}\beta^{-1}}{\rm d}{\mathrm{\mathbf{W}}}_{i,k}\,,$
(8)

where
${\mathrm{\mathbf{W}}}_{i,k}={\mathrm{\mathbf{W}}}_{i,k}(t)\in\mathbb{R}^{d}$
is a vector of uncorrelated standard Wiener processes. It can be easily seen
that (8) is an Ornstein–Uhlenbeck (OU) process where, component-wise,
$\eta^{x}_{i,k}(t)$ has mean zero (i.e., $\langle\eta^{x}_{i,k}\rangle=0$) and
time correlation function [49] of

$\langle\eta^{x}_{i,k}(t+s)\eta^{x}_{i,k}(t)\rangle=\beta^{-1}\exp\left(-\tilde{\alpha}_{k}s\right)\,,\quad
s\geq 0\,.$ (9)

Therefore, the random force ${\boldsymbol{\eta}}_{i}$ in (1b) can be rewritten
as

${\boldsymbol{\eta}}_{i}=\sum^{M}_{k=1}\tilde{\lambda}_{k}{\boldsymbol{\eta}}_{i,k}\,.$
(10)

Moreover, combining the results of (3)–(10) and introducing a new variable of
${{\mathbf{z}}}_{i,k}={\boldsymbol{\xi}}_{i,k}+{\boldsymbol{\eta}}_{i,k}$, the
GLE (1) may be rewritten as the following extended variable Markovian form:

$\displaystyle{\rm d}{{\mathbf{q}}}_{i}$
$\displaystyle=m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm d}t\,,$ (11a)
$\displaystyle{\rm d}{{\mathbf{p}}}_{i}$
$\displaystyle=-\nabla_{{{\mathbf{q}}}_{i}}U({{\mathbf{q}}}){\rm
d}t+\sum^{M}_{k=1}\tilde{\lambda}_{k}{{\mathbf{z}}}_{i,k}{\rm d}t\,,$ (11b)
$\displaystyle{\rm d}{{\mathbf{z}}}_{i,k}$
$\displaystyle=-\tilde{\lambda}_{k}m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm
d}t-\tilde{\alpha}_{k}{{\mathbf{z}}}_{i,k}{\rm
d}t+\sqrt{2\tilde{\alpha}_{k}\beta^{-1}}{\rm
d}{\mathrm{\mathbf{W}}}_{i,k}\,,\quad k=1,\dots,M\,.$ (11c)

It is worth mentioning that, in the special case of $M=1$, (11) is very
similar to a third-order Langevin dynamics that has certain advantages in
controlling the discretization errors over its corresponding second-order
form, i.e., an underdamped Langevin dynamics (see more discussions in [43]).
Note also that, by introducing an additional variable, the extended variable
Markovian form (11) can be cast into the GENERIC (General Equation for Non-
Equilibrium Reversible–Irreversible Coupling) formalism (see more discussions
in [44, 14, 26]), while this cannot be easily done in the original non-
Markovian form (1). We can write down the Fokker–Planck (or forward
Kolmogorov) equation associated with (11)

$\displaystyle\frac{\partial\rho}{\partial t}$
$\displaystyle\,=\mathcal{L}^{{\dagger}}_{\mathrm{GLE}}\rho=\sum^{N}_{i=1}\left[-m_{i}^{-1}{{\mathbf{p}}}_{i}\cdot\nabla_{{{\mathbf{q}}}_{i}}\rho+\nabla_{{{\mathbf{q}}}_{i}}U({{\mathbf{q}}})\cdot\nabla_{{{\mathbf{p}}}_{i}}\rho\right]$
(12)
$\displaystyle+\sum^{N}_{i=1}\sum^{M}_{k=1}\left[\tilde{\lambda}_{k}\left(-{{\mathbf{z}}}_{i,k}\cdot\nabla_{{{\mathbf{p}}}_{i}}\rho+m_{i}^{-1}{{\mathbf{p}}}_{i}\cdot\nabla_{{{\mathbf{z}}}_{i,k}}\rho\right)+\tilde{\alpha}_{k}\left(\nabla_{{{\mathbf{z}}}_{i,k}}\cdot\left({{\mathbf{z}}}_{i,k}\rho\right)+\beta^{-1}\Delta_{{{\mathbf{z}}}_{i,k}}\rho\right)\right]\,.$

It can then be shown that there exists a unique invariant measure defined by
the density

$\rho_{\beta}({{\mathbf{q}}},{{\mathbf{p}}},{{\mathbf{z}}})=\frac{1}{\hat{Z}}\exp\left(-\beta\left[U({{\mathbf{q}}})+\sum^{N}_{i=1}\frac{{{\mathbf{p}}}_{i}\cdot{{\mathbf{p}}}_{i}}{2m_{i}}+\sum^{N}_{i=1}\sum^{M}_{k=1}\frac{{{\mathbf{z}}}_{i,k}\cdot{{\mathbf{z}}}_{i,k}}{2}\right]\right)\,,$
(13)

where $\hat{Z}$ is the partition function. That is, (13) is the unique
solution of the stationary Fokker–Planck equation (12), i.e.,
$\mathcal{L}^{{\dagger}}_{\mathrm{GLE}}\rho_{\beta}=0$.

###  2.2 White noise limit

From (11c) we have

${{\mathbf{z}}}_{i,k}{\rm d}t=\frac{1}{\tilde{\alpha}_{k}}\left(-{\rm
d}{{\mathbf{z}}}_{i,k}-\tilde{\lambda}_{k}m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm
d}t+\sqrt{2\tilde{\alpha}_{k}\beta^{-1}}{\rm
d}{\mathrm{\mathbf{W}}}_{i,k}\right)\,,$ (14)

and substituting it into (11b) gives

${\rm d}{{\mathbf{p}}}_{i}=-\nabla_{{{\mathbf{q}}}_{i}}U({{\mathbf{q}}}){\rm
d}t-\sum^{M}_{k=1}\frac{\tilde{\lambda}_{k}}{\tilde{\alpha}_{k}}{\rm
d}{{\mathbf{z}}}_{i,k}-\sum^{M}_{k=1}\frac{\tilde{\lambda}^{2}_{k}}{\tilde{\alpha}_{k}}m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm
d}t+\sum^{M}_{k=1}\sqrt{\frac{2\tilde{\lambda}^{2}_{k}\beta^{-1}}{\tilde{\alpha}_{k}}}{\rm
d}{\mathrm{\mathbf{W}}}_{i,k}\,,$ (15)

which, in the white noise limit of $\epsilon\rightarrow 0$ (4), becomes

${\rm d}{{\mathbf{p}}}_{i}=-\nabla_{{{\mathbf{q}}}_{i}}U({{\mathbf{q}}}){\rm
d}t-\sum^{M}_{k=1}\frac{\lambda^{2}_{k}}{\alpha_{k}}m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm
d}t+\sum^{M}_{k=1}\sqrt{\frac{2\lambda^{2}_{k}\beta^{-1}}{\alpha_{k}}}{\rm
d}{\mathrm{\mathbf{W}}}_{i,k}\,.$ (16)

Rewriting the summation of ${\rm d}{\mathrm{\mathbf{W}}}_{i,k}$ in the
equation above as

$\sqrt{\sum^{M}_{k=1}\frac{2\lambda^{2}_{k}\beta^{-1}}{\alpha_{k}}}{\rm
d}{\mathrm{\mathbf{W}}}_{i}\,,$ (17)

we have

${\rm d}{{\mathbf{p}}}_{i}=-\nabla_{{{\mathbf{q}}}_{i}}U({{\mathbf{q}}}){\rm
d}t-\gamma m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm d}t+\sqrt{2\gamma\beta^{-1}}{\rm
d}{\mathrm{\mathbf{W}}}_{i}\,,$ (18)

where the friction coefficient is given by

$\gamma=\sum^{M}_{k=1}\frac{\lambda^{2}_{k}}{\alpha_{k}}\,,$ (19)

which is precisely the integral of the memory kernel defined in (3),

$\gamma=\int^{\infty}_{0}\hat{K}(t)\,{\rm d}t\,.$ (20)

Thus, we have heuristically verified that, in the white noise limit of
$\epsilon\rightarrow 0$, the solution of the GLE (11) converges weakly to that
of the Langevin dynamics (see a rigorous proof in [45]).

###  2.3 Markovian approximations of the GLE with a more general memory kernel

We have demonstrated in Section 2.1 that, for the particular form of the
memory kernel (3), the GLE (1) may be rewritten as the extended variable
Markovian form (11). We consider here the case of $N=d=1$ while dropping the
subscripts for simplicity. Moreover, for a more general form of the memory
kernel, we may be able to approximate the trajectories of $(q,p)$ in (1) by
trajectories of $(\tilde{q},\tilde{p})$ that solve a Markovian system with a
vector of $M$ auxiliary variables,
${\boldsymbol{\xi}}=\left(\xi_{1},\xi_{2},\ldots,\xi_{M}\right)^{\mathrm{T}}\in\mathbb{R}^{M}$:

$\displaystyle{\rm d}\tilde{q}$ $\displaystyle=m^{-1}\tilde{p}{\rm d}t\,,$
(21a) $\displaystyle{\rm d}\tilde{p}$ $\displaystyle=-\nabla U(\tilde{q}){\rm
d}t+{{\mathbf{g}}}^{\mathrm{T}}{\boldsymbol{\xi}}{\rm d}t\,,$ (21b)
$\displaystyle{\rm d}{\boldsymbol{\xi}}$
$\displaystyle=-m^{-1}\tilde{p}{{\mathbf{g}}}{\rm
d}t-{{\mathbf{A}}}{\boldsymbol{\xi}}{\rm d}t+{{\mathbf{C}}}{\rm
d}{\mathrm{\mathbf{W}}}\,,\quad{\boldsymbol{\xi}}(0)\sim\mathcal{N}\left(\mathbf{0},{\boldsymbol{\Sigma}}\right)\,,$
(21c)

where ${{\mathbf{g}}}\in\mathbb{R}^{M}$ is a constant vector,
$\mathbf{A},\mathbf{C}\in\mathbb{R}^{M\times M}$ are constant matrices, and
${\mathrm{\mathbf{W}}}={\mathrm{\mathbf{W}}}(t)\in\mathbb{R}^{M}$ is a vector
of uncorrelated standard Wiener processes. As demonstrated in [28], (21) would
be a Markovian approximation of (1) if
${{\mathbf{g}}},{{\mathbf{A}}},{{\mathbf{C}}}$ and ${\boldsymbol{\Sigma}}$
satisfy the following relations:

${\boldsymbol{\Sigma}}=\beta^{-1}{{\mathbf{I}}}\,,\quad{{\mathbf{C}}}{{\mathbf{C}}}^{\mathrm{T}}=\beta^{-1}\left({{\mathbf{A}}}+{{\mathbf{A}}}^{\mathrm{T}}\right)\,,$
(22)

with the approximated memory kernel being

$\tilde{K}(t)={{\mathbf{g}}}^{\mathrm{T}}e^{-{{\mathbf{A}}}t}{{\mathbf{g}}}\,.$
(23)

We can also write down the Laplace transform of the approximated memory kernel
as

$\mathcal{L}\\{\tilde{K}\\}(s)={{\mathbf{g}}}^{\mathrm{T}}\left({{\mathbf{A}}}+s{{\mathbf{I}}}\right)^{-1}{{\mathbf{g}}}\,,$
(24)

where the right-hand side is a rational function of $s$. Thus, a Markovian
approximation is established in two steps: first we approximate the Laplace
transform of the memory kernel by a rational function; then we construct a
matrix ${{\mathbf{A}}}$ and a vector ${{\mathbf{g}}}$ so that (24) holds.
Subsequently, ${\boldsymbol{\Sigma}}$ and ${{\mathbf{C}}}$ are determined from
(22).

For a general memory kernel $\hat{K}$, these tasks are nontrivial. The
simplest case is when the memory kernel can be approximated by a sum of
exponentials, e.g., the positive Prony series memory kernel (3), whose Laplace
transform can be easily computed as

$\quad\mathcal{L}\\{\tilde{K}\\}(s)=\sum_{k=1}^{M}\frac{\tilde{\lambda}^{2}_{k}}{s+\tilde{\alpha}_{k}}\,.$
(25)

In this case, ${{\mathbf{A}}}$ and ${{\mathbf{g}}}$ are chosen as

${{\mathbf{A}}}=\mathrm{diag}\left(\tilde{\alpha}_{1},\ldots,\tilde{\alpha}_{M}\right)\,,\quad{{\mathbf{g}}}=\left(\tilde{\lambda}_{1},\ldots,\tilde{\lambda}_{M}\right)^{\mathrm{T}}\,.$
(26)

Note also that in cases where ${{\mathbf{A}}}$ is diagonalizable they can be
reduced to the case above after an appropriate orthogonal transformation [28].

##  3 Numerical methods

Splitting methods have been widely used in a range of systems, including
Hamiltonian dynamics [35, 19], dissipative systems [51], and various
stochastic dynamics [31, 32, 34, 37, 38, 39, 50, 48]. The techniques have also
been adopted in the construction of numerical methods for the GLE (e.g., [9,
6, 36]). In what follows we will first review three splitting methods proposed
by Baczewski and Bond [6]. We will then propose a new method based on an
alternative splitting of the vector field. Error analysis on the averages in
the case of a one-dimensional harmonic oscillator will also be performed.

###  3.1 The BACSCAB method

The vector field of the GLE (11) can be decomposed into pieces, for instance,
A, B, C, and S:

$\displaystyle{\rm d}\left[\begin{array}[]{c}{{\mathbf{q}}}_{i}\\\
{{\mathbf{p}}}_{i}\\\ {{\mathbf{z}}}_{i,1}\\\ {{\mathbf{z}}}_{i,2}\\\
\vdots\\\ {{\mathbf{z}}}_{i,M}\end{array}\right]=$
$\displaystyle\,\underbrace{\left[\begin{array}[]{c}m_{i}^{-1}{{\mathbf{p}}}_{i}\\\
\mathbf{0}\\\ \mathbf{0}\\\ \mathbf{0}\\\ \vdots\\\
\mathbf{0}\end{array}\right]{\rm
d}t}_{\mathrm{A}}+\underbrace{\left[\begin{array}[]{c}\mathbf{0}\\\
-\nabla_{{{\mathbf{q}}}_{i}}U({{\mathbf{q}}})\\\ \mathbf{0}\\\ \mathbf{0}\\\
\vdots\\\ \mathbf{0}\end{array}\right]{\rm
d}t}_{\mathrm{B}}+\underbrace{\left[\begin{array}[]{c}\mathbf{0}\\\
\sum^{M}_{k=1}\tilde{\lambda}_{k}{{\mathbf{z}}}_{i,k}\\\ \mathbf{0}\\\
\mathbf{0}\\\ \vdots\\\ \mathbf{0}\end{array}\right]{\rm d}t}_{\mathrm{C}}$
(27) $\displaystyle+\underbrace{\left[\begin{array}[]{c}\mathbf{0}\\\
\mathbf{0}\\\ -\tilde{\lambda}_{1}m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm
d}t-\tilde{\alpha}_{1}{{\mathbf{z}}}_{i,1}{\rm
d}t+\sqrt{2\tilde{\alpha}_{1}\beta^{-1}}{\rm d}{\mathrm{\mathbf{W}}}_{i,1}\\\
-\tilde{\lambda}_{2}m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm
d}t-\tilde{\alpha}_{2}{{\mathbf{z}}}_{i,2}{\rm
d}t+\sqrt{2\tilde{\alpha}_{2}\beta^{-1}}{\rm d}{\mathrm{\mathbf{W}}}_{i,2}\\\
\vdots\\\ -\tilde{\lambda}_{M}m_{i}^{-1}{{\mathbf{p}}}_{i}{\rm
d}t-\tilde{\alpha}_{M}{{\mathbf{z}}}_{i,M}{\rm
d}t+\sqrt{2\tilde{\alpha}_{M}\beta^{-1}}{\rm
d}{\mathrm{\mathbf{W}}}_{i,M}\end{array}\right]}_{\mathrm{S}}\,,$

in such a way that each subsystem can be solved “exactly”. It is worth
mentioning that the S part consists of a vector of uncorrelated OU processes,
each of which has an exact (in the sense of distributional fidelity) solution
[25].

In describing splitting methods, we use the formal notation of the generator
of the diffusion as in [12, 47, 54]. The generators for each part of the
system may be written out as follows:

$\displaystyle\mathcal{L}_{\mathrm{A}}$
$\displaystyle=\sum^{N}_{i=1}m_{i}^{-1}{{\mathbf{p}}}_{i}\cdot\nabla_{\mathbf{q}_{i}}\,,$
(28a) $\displaystyle\mathcal{L}_{\mathrm{B}}$
$\displaystyle=-\sum^{N}_{i=1}\nabla_{\mathbf{q}_{i}}U(\mathbf{q})\cdot\nabla_{\mathbf{p}_{i}}\,,$
(28b) $\displaystyle\mathcal{L}_{\mathrm{C}}$
$\displaystyle=\sum^{N}_{i=1}\sum^{M}_{k=1}\tilde{\lambda}_{k}{{\mathbf{z}}}_{i,k}\cdot\nabla_{\mathbf{p}_{i}}\,,$
(28c) $\displaystyle\mathcal{L}_{\mathrm{S}}$
$\displaystyle=\sum^{N}_{i=1}\sum^{M}_{k=1}\left[-\tilde{\lambda}_{k}m_{i}^{-1}{{\mathbf{p}}}_{i}-\tilde{\alpha}_{k}{{\mathbf{z}}}_{i,k}+\tilde{\alpha}_{k}\beta^{-1}\nabla_{{{\mathbf{z}}}_{i,k}}\right]\cdot\nabla_{{{\mathbf{z}}}_{i,k}}\,.$
(28d)

The generator for the GLE thus can be written as

$\mathcal{L}_{\mathrm{GLE}}=\mathcal{L}_{\mathrm{A}}+\mathcal{L}_{\mathrm{B}}+\mathcal{L}_{\mathrm{C}}+\mathcal{L}_{\mathrm{S}}\,.$
(29)

Moreover, the flow map (or phase space propagator) of the system may be given
by the shorthand notation

$\mathcal{F}_{t}=e^{t\mathcal{L}_{\mathrm{GLE}}}\,,$ (30)

where the exponential map is used to formally denote the solution operator to
the equation $\partial_{t}u=\mathcal{L}_{\mathrm{GLE}}u$. Furthermore,
approximations of $\mathcal{F}_{t}$ may be obtained as products (taken in
different arrangements) of exponentials of the various splitting terms. For
instance, the phase space propagation of a splitting method in Section II E of
[6], termed BACSCAB, can be written as

$e^{\Delta t\hat{\mathcal{L}}_{\mathrm{BACSCAB}}}=e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{B}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{A}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{C}}}e^{\Delta
t\mathcal{L}_{\mathrm{S}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{C}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{A}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{B}}}\,,$ (31)

where $\exp\left(\Delta t\mathcal{L}_{f}\right)$ denotes the phase space
propagator associated with the corresponding vector field $f$. Note that the
steplengths associated with various operations are uniform and span the
interval $\Delta t$. Therefore, each of the A, B, and C steps in (31) is taken
with a steplength of $\Delta t/2$, while a steplength of $\Delta t$ is
associated with the S step. Moreover, it is worth mentioning that, while each
of the OU processes in the S step is exactly solvable, known as the “method 2”
in [6],

${{\mathbf{z}}}_{i,k}(\Delta
t)=\theta_{k}{{\mathbf{z}}}_{i,k}(0)-\left(1-\theta_{k}\right)\frac{\tilde{\lambda}_{k}m_{i}^{-1}{{\mathbf{p}}}_{i}}{\tilde{\alpha}_{k}}+\eta_{k}\sqrt{\beta^{-1}}{\mathrm{\mathbf{R}}}_{i,k}\,,$
(32)

where

$\theta_{k}=\exp\left(-\tilde{\alpha}_{k}\Delta
t\right)\,,\quad\eta_{k}=\sqrt{1-\theta^{2}_{k}}\,,$ (33)

${{\mathbf{z}}}_{i,k}(0)$ is the initial value of ${{\mathbf{z}}}_{i,k}$, and
${\mathrm{\mathbf{R}}}_{i,k}$ is a vector of independent and identically
distributed (i.i.d.) standard normal random variables. Note that an
alternative approach was used in the BACSCAB method, known as the “method 3”
in [6], by slightly modifying the $\eta_{k}$ defined above,

$\tilde{\eta}_{k}=\sqrt{\frac{2\left(1-\theta_{k}\right)^{2}}{\Delta
t\tilde{\alpha}_{k}}}\,,$ (34)

in order to improve the stability [6]. However, we would like to point that in
the current form there does not appear to have stability issues in the exact
solutions (32)–(33) as $\tilde{\alpha}_{k}\rightarrow\infty$. (The memory
kernel was written in a slightly different way in [6].) The Euler–Maruyama
method can also be used for solving the OU process, known as the “method 1” in
[6]. However, it will not be included for comparisons since it has been
observed that its performance is not as good as the other two alternative
methods. The integration steps of the BACSCAB method (31) read

$\displaystyle{{\mathbf{p}}}_{i}^{n+1/4}$
$\displaystyle={{\mathbf{p}}}_{i}^{n}-(\Delta
t/2)\nabla_{{{\mathbf{q}}}_{i}^{n}}U({{\mathbf{q}}}^{n})\,,$ (35a)
$\displaystyle{{\mathbf{q}}}_{i}^{n+1/2}$
$\displaystyle={{\mathbf{q}}}_{i}^{n}+(\Delta
t/2)m_{i}^{-1}{{\mathbf{p}}}_{i}^{n+1/4}\,,$ (35b)
$\displaystyle{{\mathbf{p}}}_{i}^{n+2/4}$
$\displaystyle={{\mathbf{p}}}_{i}^{n+1/4}+(\Delta
t/2)\sum^{M}_{k=1}\tilde{\lambda}_{k}{{\mathbf{z}}}_{i,k}^{n}\,,$ (35c)
$\displaystyle{{\mathbf{z}}}_{i,k}^{n+1}$
$\displaystyle=\theta_{k}{{\mathbf{z}}}_{i,k}^{n}-\left(1-\theta_{k}\right)\frac{\tilde{\lambda}_{k}m_{i}^{-1}{{\mathbf{p}}}_{i}^{n+2/4}}{\tilde{\alpha}_{k}}+\tilde{\eta}_{k}\sqrt{\beta^{-1}}{\mathrm{\mathbf{R}}}_{i,k}^{n}\,,$
(35d) $\displaystyle{{\mathbf{p}}}_{i}^{n+3/4}$
$\displaystyle={{\mathbf{p}}}_{i}^{n+2/4}+(\Delta
t/2)\sum^{M}_{k=1}\tilde{\lambda}_{k}{{\mathbf{z}}}_{i,k}^{n+1}\,,$ (35e)
$\displaystyle{{\mathbf{q}}}_{i}^{n+1}$
$\displaystyle={{\mathbf{q}}}_{i}^{n+1/2}+(\Delta
t/2)m_{i}^{-1}{{\mathbf{p}}}_{i}^{n+3/4}\,,$ (35f)
$\displaystyle{{\mathbf{p}}}_{i}^{n+1}$
$\displaystyle={{\mathbf{p}}}_{i}^{n+3/4}-(\Delta
t/2)\nabla_{{{\mathbf{q}}}_{i}^{n+1}}U({{\mathbf{q}}}^{n+1})\,.$ (35g)

###  3.2 The PASP method

An alternative splitting method was proposed in Section II C of [6], termed
PASP, whose phase space propagation can be written as

$\exp\left(\Delta
t\hat{\mathcal{L}}_{\mathrm{PASP}}\right)=\exp\left(\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{P}}\right)\exp\left(\Delta
t\mathcal{L}_{\mathrm{A}}\right)\exp\left(\Delta
t\mathcal{L}_{\mathrm{S}}\right)\exp\left(\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{P}}\right)\,,$ (36)

where

$\mathcal{L}_{\mathrm{P}}=\mathcal{L}_{\mathrm{B}}+\mathcal{L}_{\mathrm{C}}\,.$
(37)

Depending on how the S part is solved, there are three variants of the PASP
method. For example, if the S part is solved exactly as in (32)–(33), we
arrive the PASP-2 method, whose integration steps read

$\displaystyle{{\mathbf{p}}}_{i}^{n+1/2}$
$\displaystyle={{\mathbf{p}}}_{i}^{n}-(\Delta
t/2)\nabla_{{{\mathbf{q}}}_{i}^{n}}U({{\mathbf{q}}}^{n})+(\Delta
t/2)\sum^{M}_{k=1}\tilde{\lambda}_{k}{{\mathbf{z}}}_{i,k}^{n}\,,$ (38a)
$\displaystyle{{\mathbf{q}}}_{i}^{n+1}$
$\displaystyle={{\mathbf{q}}}_{i}^{n}+\Delta
tm_{i}^{-1}{{\mathbf{p}}}_{i}^{n+1/2}\,,$ (38b)
$\displaystyle{{\mathbf{z}}}_{i,k}^{n+1}$
$\displaystyle=\theta_{k}{{\mathbf{z}}}_{i,k}^{n}-\left(1-\theta_{k}\right)\frac{\tilde{\lambda}_{k}m_{i}^{-1}{{\mathbf{p}}}_{i}^{n+1/2}}{\tilde{\alpha}_{k}}+\eta_{k}\sqrt{\beta^{-1}}{\mathrm{\mathbf{R}}}_{i,k}^{n}\,,$
(38c) $\displaystyle{{\mathbf{p}}}_{i}^{n+1}$
$\displaystyle={{\mathbf{p}}}_{i}^{n+1/2}-(\Delta
t/2)\nabla_{{{\mathbf{q}}}_{i}^{n+1}}U({{\mathbf{q}}}^{n+1})+(\Delta
t/2)\sum^{M}_{k=1}\tilde{\lambda}_{k}{{\mathbf{z}}}_{i,k}^{n+1}\,.$ (38d)

Note that the integration steps for the PASP-3 method [6] are exactly the same
as PASP-2 except the $\eta_{k}$ in (38c) is replaced by $\tilde{\eta}_{k}$ as
in (34).

###  3.3 The BAEOEAB method

BACSCAB and both PASP methods mentioned above rely on solving the S part
exactly as in (32)–(33) or with a slight modification (34). However, one
potential drawback of that approach is, as
$\tilde{\alpha}_{k}\rightarrow\infty$, the update on ${{\mathbf{z}}}_{i,k}$
becomes increasingly dominated by the noise in (32), thereby shrinking the
contribution from the “force term” (i.e.,
$-\tilde{\lambda}_{k}m_{k}^{-1}{{\mathbf{p}}}_{k}$), which is likely to cause
potential accuracy and/or stability issues in sampling the invariant measure
(13). To this end, we propose to combine the “force term” with the original C
part in (27) to form the E part as follows:

$\displaystyle{\rm d}\left[\begin{array}[]{c}{{\mathbf{q}}}_{i}\\\
{{\mathbf{p}}}_{i}\\\ {{\mathbf{z}}}_{i,1}\\\ {{\mathbf{z}}}_{i,2}\\\
\vdots\\\ {{\mathbf{z}}}_{i,M}\end{array}\right]=$
$\displaystyle\,\underbrace{\left[\begin{array}[]{c}m_{i}^{-1}{{\mathbf{p}}}_{i}\\\
\mathbf{0}\\\ \mathbf{0}\\\ \mathbf{0}\\\ \vdots\\\
\mathbf{0}\end{array}\right]{\rm
d}t}_{\mathrm{A}}+\underbrace{\left[\begin{array}[]{c}\mathbf{0}\\\
-\nabla_{{{\mathbf{q}}}_{i}}U({{\mathbf{q}}})\\\ \mathbf{0}\\\ \mathbf{0}\\\
\vdots\\\ \mathbf{0}\end{array}\right]{\rm
d}t}_{\mathrm{B}}+\underbrace{\left[\begin{array}[]{c}\mathbf{0}\\\
\sum^{M}_{k=1}\tilde{\lambda}_{k}{{\mathbf{z}}}_{i,k}\\\
-\tilde{\lambda}_{1}m_{i}^{-1}{{\mathbf{p}}}_{i}\\\
-\tilde{\lambda}_{2}m_{i}^{-1}{{\mathbf{p}}}_{i}\\\ \vdots\\\
-\tilde{\lambda}_{M}m_{i}^{-1}{{\mathbf{p}}}_{i}\end{array}\right]{\rm
d}t}_{\mathrm{E}}$ (39)
$\displaystyle+\underbrace{\left[\begin{array}[]{c}\mathbf{0}\\\ \mathbf{0}\\\
-\tilde{\alpha}_{1}{{\mathbf{z}}}_{i,1}{\rm
d}t+\sqrt{2\tilde{\alpha}_{1}\beta^{-1}}{\rm d}{\mathrm{\mathbf{W}}}_{i,1}\\\
-\tilde{\alpha}_{2}{{\mathbf{z}}}_{i,2}{\rm
d}t+\sqrt{2\tilde{\alpha}_{2}\beta^{-1}}{\rm d}{\mathrm{\mathbf{W}}}_{i,2}\\\
\vdots\\\ -\tilde{\alpha}_{M}{{\mathbf{z}}}_{i,M}{\rm
d}t+\sqrt{2\tilde{\alpha}_{M}\beta^{-1}}{\rm
d}{\mathrm{\mathbf{W}}}_{i,M}\end{array}\right]}_{\mathrm{O}}\,.$

In this case, each of the OU processes in the O step is still exactly
solvable,

${{\mathbf{z}}}_{i,k}(\Delta
t)=\theta_{k}{{\mathbf{z}}}_{i,k}(0)+\sqrt{\beta^{-1}\left(1-\theta^{2}_{k}\right)}{\mathrm{\mathbf{R}}}_{i,k}\,.$
(40)

Moreover, we can further split the E part in such a way that each subsystem,
$\mathrm{E^{x}_{i,k}}$, component-wise,

${\rm d}\left[\begin{array}[]{c}q_{i}^{x}\\\ p_{i}^{x}\\\
z_{i,k}^{x}\end{array}\right]=\left[\begin{array}[]{c}0\\\
\tilde{\lambda}_{k}z_{i,k}^{x}\\\
-\tilde{\lambda}_{k}m_{i}^{-1}p_{i}^{x}\end{array}\right]{\rm d}t\,,$ (41)

with $q_{i}^{x}$ remaining fixed, corresponds to a harmonic oscillator with
the exact solution:

$\displaystyle p_{i}^{x}(\Delta t)$
$\displaystyle=\cos(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t)p_{i}^{x}(0)+\sin(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t)m^{1/2}_{i}z_{i,k}^{x}(0)\,,$ (42a) $\displaystyle z_{i,k}^{x}(\Delta t)$
$\displaystyle=-\sin(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t)m^{-1/2}_{i}p_{i}^{x}(0)+\cos(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t)z_{i,k}^{x}(0)\,,$ (42b)

where the superscript $x\leq d$ is a positive integer that represents a
specific dimension, $p_{i}^{x}(0)$ and $z_{i,k}^{x}(0)$ are the initial values
of $p_{i}^{x}$ and $z_{i,k}^{x}$, respectively. We further propose a new
splitting method, termed BAEOEAB, whose phase space propagation can be written
as

$e^{\Delta t\hat{\mathcal{L}}_{\mathrm{BAEOEAB}}}=e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{B}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{A}}}e^{\frac{\Delta
t}{2}\hat{\mathcal{L}}_{\mathrm{E}}}e^{\Delta
t\mathcal{L}_{\mathrm{O}}}e^{\frac{\Delta
t}{2}\hat{\mathcal{L}}_{\mathrm{E}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{A}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{B}}}\,.$ (43)

Given that the E part has been further split, the propagation of either of the
E parts in (43) should be more explicitly defined as

$e^{\frac{\Delta t}{2}\hat{\mathcal{L}}_{\mathrm{E}}}=e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{E}^{z}_{N,M}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{E}^{y}_{N,M}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{E}^{x}_{N,M}}}\dots e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{E}^{z}_{1,1}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{E}^{y}_{1,1}}}e^{\frac{\Delta
t}{2}\mathcal{L}_{\mathrm{E}^{x}_{1,1}}}\,.$ (44)

Note that one may wish to reverse the order in either of the E parts in (43),
which would affect neither its overall performance nor the order of
convergence to the invariant measure. The corresponding integration steps may
be written out as follows:

$\displaystyle{{\mathbf{p}}}_{i}^{n+1/4}$
$\displaystyle={{\mathbf{p}}}_{i}^{n}-(\Delta
t/2)\nabla_{{{\mathbf{q}}}_{i}^{n}}U({{\mathbf{q}}}^{n})\,,$ (45a)
$\displaystyle{{\mathbf{q}}}_{i}^{n+1/2}$
$\displaystyle={{\mathbf{q}}}_{i}^{n}+(\Delta
t/2)m_{i}^{-1}{{\mathbf{p}}}_{i}^{n+1/4}\,,$ (45b) $\displaystyle
p_{i}^{x,n+2/4}$ $\displaystyle=\cos(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t/2)p_{i}^{x,n+1/4}+\sin(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t/2)m^{1/2}_{i}z_{i,k}^{x,n}\,,$ (45c) $\displaystyle z_{i,k}^{x,n+1/3}$
$\displaystyle=-\sin(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t/2)m^{-1/2}_{i}p_{i}^{x,n+1/4}+\cos(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t/2)z_{i,k}^{x,n}\,,$ (45d) $\displaystyle{{\mathbf{z}}}_{i,k}^{n+2/3}$
$\displaystyle=\theta_{k}{{\mathbf{z}}}_{i,k}^{n+1/3}+\sqrt{\beta^{-1}\left(1-\theta^{2}_{k}\right)}{\mathrm{\mathbf{R}}}_{i,k}\,,$
(45e) $\displaystyle p_{i}^{x,n+3/4}$
$\displaystyle=\cos(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t/2)p_{i}^{x,n+2/4}+\sin(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t/2)m^{1/2}_{i}z_{i,k}^{x,n+2/3}\,,$ (45f) $\displaystyle z_{i,k}^{x,n+1}$
$\displaystyle=-\sin(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t/2)m^{-1/2}_{i}p_{i}^{x,n+2/4}+\cos(\tilde{\lambda}_{k}m^{-1/2}_{i}\Delta
t/2)z_{i,k}^{x,n+2/3}\,,$ (45g) $\displaystyle{{\mathbf{q}}}_{i}^{n+1}$
$\displaystyle={{\mathbf{q}}}_{i}^{n+1/2}+(\Delta
t/2)m_{i}^{-1}{{\mathbf{p}}}_{i}^{n+3/4}\,,$ (45h)
$\displaystyle{{\mathbf{p}}}_{i}^{n+1}$
$\displaystyle={{\mathbf{p}}}_{i}^{n+3/4}-(\Delta
t/2)\nabla_{{{\mathbf{q}}}_{i}^{n+1}}U({{\mathbf{q}}}^{n+1})\,.$ (45i)

As indicated in (44), both (45c)–(45d) and (45f)–(45g) loop over not only each
pair of $i=1\dots N$ and $k=1\dots M$, but also each Cartesian component
(i.e., $x,y,$ or $z$). Note that we can combine the E and O parts in (39)
together to form an OU process that can also be solved exactly as in [9, 36].
However, this approach relies on operations involving potentially very large
matrices, which could be computationally very demanding (especially when the
number of modes, $M$, is large and/or some of the parameters,
$\tilde{\lambda}_{k}$ and $\tilde{\alpha}_{k}$, are, for instance, position-
dependent) and thus we would like to avoid. It is also worth mentioning that
we can easily adopt the procedures based on the framework of the long-time
Talay–Tubaro expansion [53, 13, 31, 32, 34, 3, 2, 38, 33, 48] to analyze the
accuracy of ergodic averages (with respect to the invariant measure) in the
stochastic numerical methods previously mentioned in this section (with a
general nonlinear force), and conclude that they all have second order
convergence to the invariant measure.

###  3.4 Error analysis

As in [6], we adopt a standard test case typically used in classical Langevin
dynamics [32] by applying a one-dimensional harmonic oscillator
$U(q)=Kq^{2}/2$ ($K>0$) with a single mode (i.e., $M=1$ in (3) while dropping
the subscripts for simplicity to have $m$, $\tilde{\lambda}$, and
$\tilde{\alpha}$). In such a simple case, we can explicitly write down one
iteration of a general numerical method evolving the dynamics as

$\left[\begin{array}[]{c}q_{n+1}\\\ p_{n+1}\\\ z_{n+1}\\\
\end{array}\right]=\Psi\left[\begin{array}[]{c}q_{n}\\\ p_{n}\\\ z_{n}\\\
\end{array}\right]+\mu_{n}\,,$ (46)

where $\Psi=(\psi_{ij})$, $i,j\in\\{1,2,3\\}$, is a constant matrix and
$\mu_{n}$ is a vector of stochastic processes, whose components can be denoted
as $\mu_{n,i}$. Taking products of $q_{n+1}$, $p_{n+1}$, and $z_{n+1}$ in the
update equations above and then taking expectations on both sides of the
equations with respect to the invariant measure of the numerical method yields

$\displaystyle\langle q^{2}\rangle=$ $\displaystyle\,\psi^{2}_{11}\langle
q^{2}\rangle+\psi^{2}_{12}\langle p^{2}\rangle+\psi^{2}_{13}\langle
z^{2}\rangle+\langle\mu^{2}_{1}\rangle+2\psi_{11}\psi_{12}\langle
qp\rangle+2\psi_{11}\psi_{13}\langle qz\rangle+2\psi_{12}\psi_{13}\langle
pz\rangle\,,$ $\displaystyle\langle p^{2}\rangle=$
$\displaystyle\,\psi^{2}_{21}\langle q^{2}\rangle+\psi^{2}_{22}\langle
p^{2}\rangle+\psi^{2}_{23}\langle
z^{2}\rangle+\langle\mu^{2}_{2}\rangle+2\psi_{21}\psi_{22}\langle
qp\rangle+2\psi_{21}\psi_{23}\langle qz\rangle+2\psi_{22}\psi_{23}\langle
pz\rangle\,,$ $\displaystyle\langle z^{2}\rangle=$
$\displaystyle\,\psi^{2}_{31}\langle q^{2}\rangle+\psi^{2}_{32}\langle
p^{2}\rangle+\psi^{2}_{33}\langle
z^{2}\rangle+\langle\mu^{2}_{3}\rangle+2\psi_{31}\psi_{32}\langle
qp\rangle+2\psi_{31}\psi_{33}\langle qz\rangle+2\psi_{32}\psi_{33}\langle
pz\rangle\,,$ $\displaystyle\langle qp\rangle=$
$\displaystyle\,\psi_{11}\psi_{21}\langle
q^{2}\rangle+\psi_{12}\psi_{22}\langle p^{2}\rangle+\psi_{13}\psi_{23}\langle
z^{2}\rangle+\langle\mu_{1}\mu_{2}\rangle+\left(\psi_{11}\psi_{22}+\psi_{12}\psi_{21}\right)\langle
qp\rangle$
$\displaystyle+\left(\psi_{11}\psi_{23}+\psi_{13}\psi_{21}\right)\langle
qz\rangle+\left(\psi_{12}\psi_{23}+\psi_{13}\psi_{22}\right)\langle
pz\rangle\,,$ $\displaystyle\langle qz\rangle=$
$\displaystyle\,\psi_{11}\psi_{31}\langle
q^{2}\rangle+\psi_{12}\psi_{32}\langle p^{2}\rangle+\psi_{13}\psi_{33}\langle
z^{2}\rangle+\langle\mu_{1}\mu_{3}\rangle+\left(\psi_{11}\psi_{32}+\psi_{12}\psi_{31}\right)\langle
qp\rangle$
$\displaystyle+\left(\psi_{11}\psi_{33}+\psi_{13}\psi_{31}\right)\langle
qz\rangle+\left(\psi_{12}\psi_{33}+\psi_{13}\psi_{32}\right)\langle
pz\rangle\,,$ $\displaystyle\langle pz\rangle=$
$\displaystyle\,\psi_{21}\psi_{31}\langle
q^{2}\rangle+\psi_{22}\psi_{32}\langle p^{2}\rangle+\psi_{23}\psi_{33}\langle
z^{2}\rangle+\langle\mu_{2}\mu_{3}\rangle+\left(\psi_{21}\psi_{32}+\psi_{22}\psi_{31}\right)\langle
qp\rangle$
$\displaystyle+\left(\psi_{21}\psi_{33}+\psi_{23}\psi_{31}\right)\langle
qz\rangle+\left(\psi_{22}\psi_{33}+\psi_{23}\psi_{32}\right)\langle
pz\rangle\,,$

where $\langle\cdot\rangle$ denotes an ensemble average with respect to the
invariant measure of the numerical method. See examples in the case of
Langevin dynamics in Chapter 7 of [33]. While the analytical expressions to
the averages on the left-hand side above are given respectively by

$\langle q^{2}\rangle=\frac{1}{K\beta}\,,\quad\langle
p^{2}\rangle=\frac{m}{\beta}\,,\quad\langle
z^{2}\rangle=\frac{1}{\beta}\,,\quad\langle qp\rangle=0\,,\quad\langle
qz\rangle=0\,,\quad\langle pz\rangle=0\,,$ (47)

we can also solve the above linear system to obtain the averages associated
with each numerical method. For brevity, some results are shown as leading
order series in $\Delta t$.

####  3.4.1 PASP-2

In the PASP-2 method (36), we have the following averages:

$\displaystyle\langle q^{2}\rangle$
$\displaystyle=\frac{1}{K\beta}\left[1+\frac{\Delta
t^{2}\left(m\tilde{\alpha}^{2}+3K\right)}{12m}\right]+O(\Delta t^{4})\,,$
(48a) $\displaystyle\langle p^{2}\rangle$
$\displaystyle=\frac{m}{\beta}\left[1+\frac{\Delta
t^{2}\tilde{\alpha}^{2}}{12}\right]+O(\Delta t^{4})\,,$ (48b)
$\displaystyle\langle z^{2}\rangle$
$\displaystyle=\frac{1}{\beta}\left[1+\frac{\Delta
t^{2}\tilde{\lambda}^{2}}{4m}\right]+O(\Delta t^{4})\,,$ (48c)
$\displaystyle\langle qp\rangle$ $\displaystyle=\langle pz\rangle=0\,,$ (48d)
$\displaystyle\langle qz\rangle$ $\displaystyle=-\frac{\Delta
t^{2}\tilde{\lambda}}{4m\beta}+O(\Delta t^{4})\,.$ (48e)

While the averages of $\langle qp\rangle$ and $\langle pz\rangle$ are exact,
there exist second order errors in the averages of $\langle q^{2}\rangle$,
$\langle p^{2}\rangle$, $\langle z^{2}\rangle$, and $\langle qz\rangle$. Note
also that the minus sign of the leading order term in $\langle qz\rangle$
appeared to be missing in [6].

####  3.4.2 PASP-3

Similarly in the PASP-3 method, we have the following averages:

$\displaystyle\langle q^{2}\rangle$
$\displaystyle=\frac{1}{K\beta}\left[1+\frac{\Delta
t^{2}K}{4m}\right]+O(\Delta t^{4})\,,$ (49a) $\displaystyle\langle
p^{2}\rangle$ $\displaystyle=\frac{m}{\beta}\,,$ (49b) $\displaystyle\langle
z^{2}\rangle$ $\displaystyle=\frac{1}{\beta}\left[1+\frac{\Delta
t^{2}\left(3\tilde{\lambda}^{2}-m\tilde{\alpha}^{2}\right)}{12m}\right]+O(\Delta
t^{4})\,,$ (49c) $\displaystyle\langle qp\rangle$ $\displaystyle=\langle
pz\rangle=0\,,$ (49d) $\displaystyle\langle qz\rangle$
$\displaystyle=-\frac{\Delta t^{2}\tilde{\lambda}}{4m\beta}+O(\Delta
t^{4})\,.$ (49e)

Similar to the case of the PASP-2 method, the averages of $\langle qp\rangle$
and $\langle pz\rangle$ are exact, while there exist second order errors in
the averages of $\langle q^{2}\rangle$, $\langle p^{2}\rangle$, $\langle
z^{2}\rangle$, and $\langle qz\rangle$. Again, the minus sign of the leading
order term in $\langle qz\rangle$ appeared to be missing in [6].

####  3.4.3 BACSCAB

In the BACSCAB method (31), we have the following averages:

$\displaystyle\langle q^{2}\rangle$ $\displaystyle=\frac{1}{K\beta}\,,$ (50a)
$\displaystyle\langle p^{2}\rangle$
$\displaystyle=\frac{m}{\beta}\left[1-\frac{\Delta t^{2}K}{4m}\right]\,,$
(50b) $\displaystyle\langle z^{2}\rangle$
$\displaystyle=\frac{1}{\beta}\left[1+\frac{\Delta
t^{2}\left(3\tilde{\lambda}^{2}-m\tilde{\alpha}^{2}\right)}{12m}\right]+O(\Delta
t^{4})\,,$ (50c) $\displaystyle\langle qp\rangle$ $\displaystyle=\langle
qz\rangle=\langle pz\rangle=0\,.$ (50d)

While the averages of $\langle q^{2}\rangle$, $\langle qp\rangle$, $\langle
qz\rangle$, and $\langle pz\rangle$ are exact, there exist second order errors
in the averages of $\langle p^{2}\rangle$ and $\langle z^{2}\rangle$.

####  3.4.4 BAEOEAB

In the BAEOEAB method (43), we have the following averages:

$\langle q^{2}\rangle=\frac{1}{K\beta}\,,\quad\langle
p^{2}\rangle=\frac{m}{\beta}\left[1-\frac{\Delta
t^{2}K}{4m}\right]\,,\quad\langle z^{2}\rangle=\frac{1}{\beta}\,,\quad\langle
qp\rangle=\langle qz\rangle=\langle pz\rangle=0\,.$ (51)

That is, all of the averages are exact, with the only exception of $\langle
p^{2}\rangle$ that has second order errors, which yields a friction-
independent upper bound of the stepsize $\Delta t_{\mathrm{max}}=2\sqrt{m/K}$,
coinciding with the (deterministic) Verlet stability threshold [32]. Note that
a similar method to BAEOEAB (i.e., BAOEOAB) yields exactly the same averages,
but requires additional generations of random numbers. It is also worth
mentioning that the averages of $\langle q^{2}\rangle$, $\langle
p^{2}\rangle$, and $\langle qp\rangle$ match perfectly with their counterparts
in the BAOAB method of the underdamped Langevin dynamics [32].

##  4 Numerical experiments

In this section, a variety of numerical experiments are conducted to
systematically compare the newly proposed BAEOEAB method (43) with alternative
methods described in Section 3. While all the methods presented can be used in
multi-mode systems, we focus our attention on a single mode (i.e., $M=1$) for
simplicity where the subscripts in $\tilde{\lambda}$ and $\tilde{\alpha}$ (3)
can be dropped.

###  4.1 Harmonic oscillator

In order to verify the error analysis results in Section 3.4, we compare the
performance of various methods with a one-dimensional harmonic oscillator
$U(q)=Kq^{2}/2$ ($K>0$). The following set of parameters were used:
$m=K=\beta=\tilde{\alpha}=1$ and $\tilde{\lambda}=2$.

Figure 1 compares the control of the relative errors in the computed averages
of $\langle q^{2}\rangle$ and $\langle z^{2}\rangle$. According to the dashed
order lines, both PSAP-2 and PASP-3 methods exhibit second order convergence
in both averages whereas the BACSCAB method appears to be second order only in
the average of $\langle z^{2}\rangle$ on the right panel. Moreover, the
BAEOEAB method is exact in both averages while the BACSCAB method is exact
only in the average of $\langle q^{2}\rangle$ on the left panel. It is worth
mentioning that the error in cases where the corresponding averages are exact
comes solely from the sampling error, rather than the discretization error.
Note also that all the observations are consistent with our findings in
Section 3.4.

More precisely, for the average of $\langle q^{2}\rangle$ on the left panel of
Figure 1, PASP-3 appears to be only slightly more accurate than PASP-2 while
both methods are significantly outperformed by either BACSCAB or BAEOEAB in
terms of accuracy. It is worth pointing out that the BAEOEAB method also
appears to be more robust than alternative methods in a wide range of
stepsizes—while both PASP-2 and PASP-3 methods exhibit around 50% relative
errors at a stepsize of around $\Delta t=0.75$ and the BACSCAB method became
unstable just over $\Delta t=1$, the BAEOEAB method is still extremely
accurate (i.e., up to sampling error) when the stepsize is around $\Delta
t=1.9$, close to its stability threshold $\Delta t_{\mathrm{max}}=2$. The
behavior is largely similar for the average of $\langle z^{2}\rangle$ on the
right panel of Figure 1 except the performance of PASP-2, PASP-3, and BACSCAB
are almost indistinguishable. In this case, the BAEOEAB method substantially
outperforms alternative methods in terms of not only accuracy but also
robustness.

The distributions of positions associated with different stepsizes were
compared in Figure 2. The behavior is consistent with our findings on the left
panel of Figure 1. That is, the distributions of both PASP-2 and PASP-3
methods visibly deviate from the reference solutions with a stepsize of
$\Delta t=0.7$ and substantial deviations can be observed with a stepsize of
$\Delta t=0.9$. In contrast, the distributions of both BACSCAB and BAEOEAB
methods are both indistinguishable from the reference solutions; however the
largest stepsizes used in the BACSCAB and BAEOEAB methods were $\Delta t=1$
and $\Delta t=1.9$, respectively. This again demonstrates the superior
accuracy and robustness of the BAEOEAB method over alternative methods.

###  4.2 Multi-particle system

Figure 1: Double logarithmic plot of the relative error in the computed
averages of $\langle q^{2}\rangle$ (left) and $\langle z^{2}\rangle$ (right)
against stepsize by using various numerical methods for the GLE described in
Section 3 with parameters of $\tilde{\lambda}=2$ and $\tilde{\alpha}=1$. The
system was simulated for $10^{7}$ reduced time units but only the last 80% of
the data were collected to calculate the static quantity in order to make sure
the system was well equilibrated. 1000 different runs were averaged to reduce
the sampling errors. The stepsizes tested began at $\Delta t=0.25$ and were
increased incrementally by 20% until all methods either started to show
significant relative errors or became unstable. The dashed black line
represents the second order convergence to the invariant measure.

Figure 2: Comparisons of the distributions obtained from various methods for
the GLE with parameters of $\tilde{\lambda}=2$ and $\tilde{\alpha}=1$. The
solid black line is the reference solution obtained by using the BACSCAB
method with a very small stepsize of $\Delta t=0.01$, while the colored lines
correspond to different stepsizes as indicated. (For interpretation of the
colors in the figure(s), the reader is referred to the web version of this
article.)

Figure 3: Double logarithmic plot of the relative error in the computed
configurational temperature against stepsize by using various numerical
methods for the GLE described in Section 3 with parameters of
$\tilde{\lambda}=1$ and $\tilde{\alpha}=16$ (top left), $\tilde{\alpha}=32$
(top right), $\tilde{\alpha}=64$ (bottom left), and $\tilde{\alpha}=128$
(bottom right). The system was simulated for $1000$ reduced time units but
only the last 80% of the data were collected to calculate the static quantity
in order to make sure the system was well equilibrated. Ten different runs
were averaged to reduce the sampling errors. The stepsizes tested began at
$\Delta t=0.03$ and were increased incrementally by 15% until all methods
either started to show significant relative errors or became unstable.

Figure 4: Double logarithmic plot of the relative error in the computed
configurational temperature against stepsize by using various numerical
methods for the GLE described in Section 3 with parameters of
$\tilde{\lambda}=1$ and $\tilde{\alpha}=16$ (top left), $\tilde{\alpha}=32$
(top right), $\tilde{\alpha}=64$ (bottom left), and $\tilde{\alpha}=128$
(bottom right). The format of the plots is the same as in Figure 3, except
with a shear rate of $\kappa=0.1$ generated by the Lees–Edwards boundary
conditions.

We also compare the performance of various methods in systems where multiple
particles are interacting with each other. To this end, we chose a soft pair
potential energy that has been widely used in the so-called dissipative
particle dynamics [37, 39, 48],

$U(\mathbf{q})=\sum_{i}\sum_{j>i}\varphi(r_{ij})\,,$ (52)

where

$\varphi(r_{ij})=\begin{cases}a_{ij}r_{\mathrm{c}}\left(1-r_{ij}/r_{\mathrm{c}}\right)^{2}/2\,,&r_{ij}<r_{\mathrm{c}}\,,\\\
\quad\quad\quad\quad 0\,,&r_{ij}\geq r_{\mathrm{c}}\,,\end{cases}$ (53)

where parameter $a_{ij}$ denotes the maximum repulsion strength between
particles $i$ and $j$,
$r_{ij}=|\mathbf{q}_{ij}|=|\mathbf{q}_{i}-\mathbf{q}_{j}|$ is the distance,
and $r_{\mathrm{c}}$ represents a certain cutoff radius, beyond which there is
no interaction between particles. We adopted a standard set of parameters
commonly used in algorithms tests as in [48]: $m_{i}=r_{\mathrm{c}}=\beta=1$.
Moreover, a particle density of $\rho_{\rm d}=3$ was used throughout the
current article, which also determines the repulsion parameter of
$a_{ij}=75/(\beta\rho_{\rm d})=25$ in order to match the compressibility of
water [18]. A system of $N=500$ identical particles was simulated in a cubic
box with periodic boundary conditions [5], unless otherwise stated. The
initial positions of the particles were i.i.d. with a uniform distribution
over the box, while the initial momenta were i.i.d. normal random variables
with mean zero and variance $1/\beta$.

The average of the computed configurational temperature in the canonical
ensemble is expected to be precisely the target temperature:

$\beta^{-1}=\frac{\left\langle\nabla_{i}U(\mathbf{q})\cdot\nabla_{i}U(\mathbf{q})\right\rangle}{\left\langle\nabla^{2}_{i}U(\mathbf{q})\right\rangle}\,,$
(54)

where $\nabla_{i}U$ and $\nabla^{2}_{i}U$, respectively, represent the
gradient and Laplacian of the potential energy $U$ with respect to the
position of particle $i$. The control of the configurational temperature has
been recommended in [4] as a verification of equilibrium; importantly, good
control of the configurational temperature also appears to imply good
performance in other configuration-based physical quantities (see more
discussions on the configurational temperature in [37, 39, 48]).

The control of the configurational temperature was compared in Figure 3 with
$\tilde{\lambda}=1$ and varying $\tilde{\alpha}$. All the methods appear to
have second order convergence to the invariant measure (dashed order lines not
shown). On the top left panel where $\tilde{\alpha}=16$, the PASP-3 method
appears to be more accurate than the PASP-2 method; however it is outperformed
by either BACSCAB or BAEOEAB methods, whose curves are almost
indistinguishable. As we increase the value of $\tilde{\alpha}$, the behavior
is largely similar except (a) PASP-2 becomes less and less accurate than
alternative methods and displays substantial relative error at a stepsize of
around $\Delta t=0.05$ on the bottom left panel where $\tilde{\alpha}=64$; (b)
BAEOEAB becomes increasingly better than alternative methods in terms of not
only accuracy but also robustness—while, on the bottom right panel where
$\tilde{\alpha}=128$, both PASP-3 and BACSCAB show an around 25% relative
error with a stepsize of slightly over $\Delta t=0.05$, the relative error is
around 36% for BAEOEAB with a stepsize of around $\Delta t=0.08$.

Figure 4 also compares the control of the configurational temperature when the
well-known Lees–Edwards boundary conditions [29] were applied in order to
generate a simple and steady shear flow (shear rate $\kappa=0.1$), again with
$\tilde{\lambda}=1$ and varying $\tilde{\alpha}$. The behavior is very similar
to Figure 3 except the superiority of the BAEOEAB method is even more evident
particularly in the bottom row where $\tilde{\alpha}$ is relatively large. For
instance, on the bottom right panel where $\tilde{\alpha}=128$, BAEOEAB
achieves around an order of magnitude improvement over both PASP-3 and BACSCAB
in terms of accuracy. This indicates that BAEOEAB has a even better
configurational temperature control over alternative methods in nonequilibrium
simulations than in equilibrium when $\tilde{\alpha}$ is relatively large. It
is worth mentioning that, although BAEOEAB has been designed in equilibrium
settings where the invariant measure is preserved, it appears that its
performance is also very good in nonequilibrium, particularly when the shear
rate is relatively small, which may be thought of as in “near-equilibrium”.

##  5 Conclusions

We have reviewed the constructions of popular splitting methods proposed for
the GLE. Having also pointed the potential drawback of those existing methods,
we have proposed an alternative method based on a different splitting of the
vector field of the GLE. We have also analyzed the errors on the averages
associated with different methods in the case of a one-dimensional harmonic
oscillator. We have demonstrated that all of the averages are exact for the
newly proposed BAEOEAB method, with the only exception being the average of
$\langle p^{2}\rangle$ that has second order errors, leading to a friction-
independent upper bound of the stepsize $\Delta t_{\mathrm{max}}=2\sqrt{m/K}$
that coincides with the (deterministic) Verlet stability threshold. In
contrast, there are at least two averages that are not exact for all the other
splitting methods examined.

Restricting our attention to a single mode (i.e., $M=1$) for simplicity, we
have systematically compared the BAEOEAB method with alternative methods in a
variety of numerical experiments. In the case of a one-dimensional harmonic
oscillator, the BAEOEAB method is exact in both averages of $\langle
q^{2}\rangle$ and $\langle z^{2}\rangle$ and substantially outperforms
alternative methods in terms of not only accuracy but also robustness. The
BACSCAB method is also exact in the average of $\langle q^{2}\rangle$, and
appears to be as accurate as the BAEOEAB method. However, the former became
unstable just over $\Delta t=1$, while the latter is still extremely accurate
(i.e., up to sampling error) when the stepsize is around $\Delta t=1.9$, close
to its stability threshold $\Delta t_{\mathrm{max}}=2$. In the average of
$\langle z^{2}\rangle$, the BACSCAB method becomes second order as both PASP
methods. Moreover, the performance of the three methods is very similar to
each other, all outperformed by the BAEOEAB method.

We have also examined the case of multiple particles that are interacting with
each other—a soft pair potential energy widely used in the so-called
dissipative particle dynamics. Fixing $\tilde{\lambda}$ to be unity, as we
increase the value of $\tilde{\alpha}$, the BAEOEAB method becomes
increasingly better than alternative methods in terms of not only accuracy but
also robustness. Moreover, the superiority of the BAEOEAB method is even more
evident when the well-known Lees–Edwards boundary conditions were applied in
order to generate a simple and steady shear flow, a technique commonly used in
nonequilibrium simulations.

Although we have focused our attention on a single mode in this article, it is
expected that the newly proposed BAEOEAB method will perform well in multi-
mode cases (i.e., a general Prony series) that can be viewed as a summation of
single modes. However, we leave a detailed examination (including more general
forms of the memory kernel, for instance, in Section 2.3) for future work.

Following similar results in Langevin dynamics [31, 38], a fourth-order
convergence to the invariant measure has been demonstrated recently for a
particular splitting method for the GLE with certain choices of the parameters
[36]. However, such superconvergence results have not been generally observed
in our numerical experiments. Therefore, we leave further exploration of this
observation for future work.

## CRediT authorship contribution statement

Manh Hong Duong: Conceptualization, Investigation, Methodology, Writing -
review & editing. Xiaocheng Shang: Conceptualization, Formal analysis,
Investigation, Methodology, Software, Validation, Visualization, Writing -
original draft, Writing - review & editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

The authors thank Nawaf Bou-Rabee, Benedict Leimkuhler, Matthias Sachs, and
the anonymous referees for their valuable suggestions and comments. XS
acknowledges the support of the London Mathematical Society through the
Research Reboot Grant (reference number 42025), the Royal Society through the
International Exchanges Scheme (reference number
IES$\backslash$R3$\backslash$203007), and the Institute of Mathematics and its
Applications through the QJMAM Fund for Applied Mathematics. The research of
MHD was supported by EPSRC Grants EP/W008041/1 and EP/V038516/1.

## References

  * [1] J. Abate and W. Whitt.  Infinite-series representations of Laplace transforms of probability density functions for numerical inversion.  J. Oper. Res. Soc. Jpn., 42(3):268–285, 1999\.
  * [2] A. Abdulle, G. Vilmart, and K. C. Zygalakis.  High order numerical approximation of the invariant measure of ergodic SDEs.  SIAM J. Numer. Anal., 52(4):1600–1622, 2014\.
  * [3] A. Abdulle, G. Vilmart, and K. C. Zygalakis.  Long time accuracy of Lie–Trotter splitting methods for Langevin dynamics.  SIAM J. Numer. Anal., 53(1):1–16, 2015.
  * [4] M. P. Allen.  Configurational temperature in membrane simulations using dissipative particle dynamics.  J. Phys. Chem. B, 110(8):3823–3830, 2006.
  * [5] M. P. Allen and D. J. Tildesley.  Computer Simulation of Liquids, Second Edition.  Oxford University Press, Oxford, 2017.
  * [6] A. D. Baczewski and S. D. Bond.  Numerical integration of the extended variable generalized Langevin equation with a positive Prony representable memory kernel.  J. Chem. Phys., 139(4):044107, 2013.
  * [7] H. B. Callen and T. A. Welton.  Irreversibility and generalized noise.  Phys. Rev., 83(1):34–40, 1951.
  * [8] M. Ceriotti, G. Bussi, and M. Parrinello.  Nuclear quantum effects in solids using a colored-noise thermostat.  Phys. Rev. Lett., 103(3):030603, 2009.
  * [9] M. Ceriotti, G. Bussi, and M. Parrinello.  Colored-noise thermostats à la carte.  J. Chem. Theory Comput., 6(4):1170–1180, 2010\.
  * [10] M. Chak, N. Kantas, and G. A. Pavliotis.  On the generalised Langevin equation for simulated annealing.  arXiv:2003.06448, 2020.
  * [11] A. Córdoba, T. Indei, and J. D. Schieber.  Elimination of inertia from a generalized Langevin equation: Applications to microbead rheology modeling and data analysis.  J. Rheol., 56(1):185–212, 2012.
  * [12] G. De Fabritiis, M. Serrano, P. Español, and P. V. Coveney.  Efficient numerical integrators for stochastic models.  Phys. A, 361(2):429–440, 2006.
  * [13] A. Debussche and E. Faou.  Weak backward error analysis for SDEs.  SIAM J. Numer. Anal., 50(3):1735–1752, 2012\.
  * [14] M. H. Duong, M. A. Peletier, and J. Zimmer.  GENERIC formalism of a Vlasov–Fokker–Planck equation and connection to large-deviation principles.  Nonlinearity, 26(11):2951, 2013.
  * [15] J. Fricks, L. Yao, T. C. Elston, and M. G. Forest.  Time-domain methods for diffusive transport in soft matter.  SIAM J. Appl. Math., 69(5):1277–1308, 2009\.
  * [16] N. E. Glatt-Holtz, D. P. Herzog, S. A. McKinley, and H. D. Nguyen.  The generalized Langevin equation with power-law memory in a nonlinear potential well.  Nonlinearity, 33(6):2820, 2020.
  * [17] F. Grogan, H. Lei, X. Li, and N. A. Baker.  Data-driven molecular modeling with the generalized Langevin equation.  J. Comput. Phys., 418:109633, 2020.
  * [18] R. D. Groot and P. B. Warren.  Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation.  J. Chem. Phys., 107(11):4423–4435, 1997.
  * [19] E. Hairer, C. Lubich, and G. Wanner.  Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations.  Springer, Berlin, 2006.
  * [20] E. J. Hall, M. A. Katsoulakis, and L. Rey-Bellet.  Uncertainty quantification for generalized Langevin dynamics.  J. Chem. Phys., 145(22):224108, 2016.
  * [21] C. Hijón, P. Español, E. Vanden-Eijnden, and R. Delgado-Buscalioni.  Mori–Zwanzig formalism as a practical computational tool.  Faraday Discuss., 144:301–322, 2010.
  * [22] G. Jung, M. Hanke, and F. Schmid.  Generalized Langevin dynamics: Construction and numerical integration of non-Markovian particle-based models.  Soft Matter, 14:9368–9382, 2018.
  * [23] L. Kantorovich.  Generalized Langevin equation for solids. I. Rigorous derivation and main properties.  Phys. Rev. B, 78(9):094304, 2008.
  * [24] L. Kantorovich and N. Rompotis.  Generalized Langevin equation for solids. II. Stochastic boundary conditions for nonequilibrium molecular dynamics simulations.  Phys. Rev. B, 78(9):094305, 2008.
  * [25] P. E. Kloeden and E. Platen.  Numerical Solution of Stochastic Differential Equations.  Springer, Berlin, 1992.
  * [26] R. C. Kraaij, A. Lazarescu, C. Maes, and M. Peletier.  Fluctuation symmetry leads to GENERIC equations with non-quadratic dissipation.  Stochastic Process. Appl., 130(1):139–170, 2020\.
  * [27] R. Kubo.  The fluctuation-dissipation theorem.  Rep. Progr. Phys., 29(1):255, 1966.
  * [28] R. Kupferman.  Fractional kinetics in kac–zwanzig heat bath models.  Journal of Statistical Physics, 114(1):291–326, 2004.
  * [29] A. W. Lees and S. F. Edwards.  The computer study of transport processes under extreme conditions.  J. Phys. C: Solid State Phys., 5(15):1921–1929, 1972.
  * [30] H. Lei, N. A. Baker, and X. Li.  Data-driven parameterization of the generalized Langevin equation.  Proc. Natl. Acad. Sci. USA, 113(50):14183–14188, 2016.
  * [31] B. Leimkuhler and C. Matthews.  Rational construction of stochastic numerical methods for molecular sampling.  Appl. Math. Res. Express, 2013(1):34–56, 2013\.
  * [32] B. Leimkuhler and C. Matthews.  Robust and efficient configurational molecular sampling via Langevin dynamics.  J. Chem. Phys., 138:174102, 2013.
  * [33] B. Leimkuhler and C. Matthews.  Molecular Dynamics: With Deterministic and Stochastic Numerical Methods.  Springer, Cham, Switzerland, 2015.
  * [34] B. Leimkuhler, C. Matthews, and G. Stoltz.  The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics.  IMA J. Numer. Anal., 36(1):13–79, 2016.
  * [35] B. Leimkuhler and S. Reich.  Simulating Hamiltonian Dynamics.  Cambridge University Press, Cambridge, 2005.
  * [36] B. Leimkuhler and M. Sachs.  Efficient numerical algorithms for the generalized Langevin equation.  SIAM J. Sci. Comput., 44(1):A364–A388, 2022\.
  * [37] B. Leimkuhler and X. Shang.  On the numerical treatment of dissipative particle dynamics and related systems.  J. Comput. Phys., 280:72–95, 2015.
  * [38] B. Leimkuhler and X. Shang.  Adaptive thermostats for noisy gradient systems.  SIAM J. Sci. Comput., 38(2):A712–A736, 2016\.
  * [39] B. Leimkuhler and X. Shang.  Pairwise adaptive thermostats for improved accuracy and stability in dissipative particle dynamics.  J. Comput. Phys., 324:174–193, 2016.
  * [40] S. A. McKinley and H. D. Nguyen.  Anomalous diffusion and the generalized Langevin equation.  SIAM J. Math. Anal., 50(5):5119–5160, 2018\.
  * [41] S. A. McKinley, L. Yao, and M. G. Forest.  Transient anomalous diffusion of tracer particles in soft matter.  J. Rheol., 53(6):1487–1506, 2009.
  * [42] H. Mori.  Transport, collective motion, and Brownian motion.  Progr. Theoret. Phys., 33(3):423–455, 1965\.
  * [43] W. Mou, Y.-A. Ma, M. J. Wainwright, P. L. Bartlett, and M. I. Jordan.  High-order Langevin diffusion yields an accelerated MCMC algorithm.  J. Mach. Learn. Res., 22(42):1–41, 2021.
  * [44] H. C. Öttinger.  Beyond Equilibrium Thermodynamics.  John Wiley & Sons, 2005.
  * [45] M. Ottobre and G. A. Pavliotis.  Asymptotic analysis for the generalized Langevin equations.  Nonlinearity, 24:1629–1653, 2011.
  * [46] G. A. Pavliotis.  Stochastic Processes and Applications: Diffusion Processes, the Fokker–Planck and Langevin Equations.  Springer, 2014.
  * [47] M. Serrano, G. De Fabritiis, P. Español, and P. V. Coveney.  A stochastic Trotter integration scheme for dissipative particle dynamics.  Math. Comput. Simulation, 72(2):190–194, 2006\.
  * [48] X. Shang.  Accurate and efficient splitting methods for dissipative particle dynamics.  SIAM J. Sci. Comput., 43(3):A1929–A1949, 2021\.
  * [49] X. Shang and M. Kröger.  Time correlation functions of equilibrium and nonequilibrium Langevin dynamics: Derivations and numerics using random numbers.  SIAM Rev., 62(4):901–935, 2020.
  * [50] X. Shang, M. Kröger, and B. Leimkuhler.  Assessing numerical methods for molecular and particle simulation.  Soft Matter, 13:8565–8578, 2017.
  * [51] X. Shang and H. C. Öttinger.  Structure-preserving integrators for dissipative systems based on reversible–irreversible splitting.  Proc. R. Soc. A, 476(2234):20190446, 2020.
  * [52] L. Stella, C. D. Lorenz, and L. Kantorovich.  Generalized Langevin equation: An efficient approach to nonequilibrium molecular dynamics of open systems.  Phys. Rev. B, 89(13):134303, 2014.
  * [53] D. Talay and L. Tubaro.  Expansion of the global error for numerical schemes solving stochastic differential equations.  Stoch. Anal. Appl., 8(4):483–509, 1990.
  * [54] F. Thalmann and J. Farago.  Trotter derivation of algorithms for Brownian and dissipative particle dynamics.  J. Chem. Phys., 127:124109, 2007.
  * [55] R. Zwanzig.  Nonequilibrium Statistical Mechanics.  Oxford University Press, 2001.