# A stochastic Hamiltonian formulation applied to dissipative particle
dynamics

Linyu Peng<EMAIL_ADDRESS>Noriyoshi Arai<EMAIL_ADDRESS>Kenji
Yasuoka<EMAIL_ADDRESS>Department of Mechanical Engineering, Keio
University, Yokohama 223-8522, Japan

###### Abstract

In this paper, a stochastic Hamiltonian formulation (SHF) is proposed and
applied to dissipative particle dynamics (DPD) simulations. As an extension of
Hamiltonian dynamics to stochastic dissipative systems, the SHF provides
necessary foundations and great convenience for constructing efficient
numerical integrators. As a first attempt, we develop the Störmer–Verlet type
of schemes based on the SHF, which are structure-preserving for deterministic
Hamiltonian systems without external forces, the dissipative forces in DPD.
Long-time behaviour of the schemes is shown numerically by studying the damped
Kubo oscillator. In particular, the proposed schemes include the conventional
Groot–Warren’s modified velocity-Verlet method and a modified version of
Gibson–Chen–Chynoweth as special cases. The schemes are applied to DPD
simulations and analysed numerically.

Keywords: Dissipative particle dynamics; Hamiltonian mechanics; Stochastic
differential equations; Störmer–Verlet methods

##  1 Introduction

A dissipative particle dynamics (DPD) simulation [1, 2] method is a type of
coarse-grained molecular simulation method, which has proven to be a powerful
tool for investigating fluid events occurring on a wide range of spatio-
temporal scales compared to all-atom simulations. Using DPD method, many
studies have been conducted for both the statics and dynamics of complex
system at the mesoscopic level, such as unique self-assembled structures
formed by nanoparticles or polymers [3, 4, 5, 6], mechanical or rheological
properties of soft materials [7, 8, 9], medical materials and biological
functions [10, 11, 12, 13], and so forth. Huang $et~{}al.$ [4] proposed a
method to fabricate various two-dimensional nanostructures using self-assembly
of block copolymers and demonstrated it in DPD simulations. The simulations
showed that surface patterns of three-dimensional nanostructures could be
evolved to solve problems in lithography and transistors. In order to overcome
the problem of low toughness in the use of humanoid robotic hands, Pan
$et~{}al.$ [8] developed an ultra-tough electric tendon based on spider silk
toughened with single-wall carbon nanotubes (SWCNTs). In that study, DPD
simulations were performed to understand how SWCNTs improve the mechanical
properties of the fibers at a molecular level. Sicard and Toro-Mendoza [13]
reported on the computational design of soft nanocarriers using pickering
emulsions (nanoparticle armored droplet), able to selectively encapsulate or
release a probe load under specific flow conditions. They described in detail
the mechanisms at play in the formation of pocket-like structures and their
stability under external flow. Moreover, the rheological properties of the
designed nanocarriers were compared with those of delivery systems used in
pharmaceutical and cosmetic technologies.

On the other hand, during the last decades, a lot of efforts have been made
for proposing efficient simulation methods for DPD to achieve simultaneous
temperature control and momentum preservation. Examples include Groot–Warren’s
modified velocity-Verlet (GW) method [2], the method of Gibson–Chen–Chynoweth
(GCC) [14], and splitting methods [15, 16]; a review and comparison of
commonly used methods for DPD are available in [17]. In the current study, we
will show that various velocity-Verlet methods for DPD, including GW and GCC
methods, are actually special cases of the Störmer–Verlet (SV) schemes for a
novel stochastic Hamiltonian formulation (SHF) with dissipative forces which
are often called external forces in classical Hamiltonian mechanics; in DPD,
these dissipative forces are in fact internal forces (see Section 2). To be
consistent, they will be called external forces in the general setting but
dissipative forces in DPD. SV schemes are well-known symplectic-preserving
numerical methods for deterministic Hamiltonian systems without external
forces.

Symplecticity is a crucial feature of Hamiltonian systems. Geometrically, it
implies area or volume preservation of the corresponding phase flows due to
Liouville’s Theorem. Symplectic integrators are among the most important types
of geometric numerical integrators for Hamiltonian systems [18, 19].
Symplectic integrators for stochastic Hamiltonian systems with or without
external forces have received great attention as well, e.g., [20, 21, 22, 23,
24]. The SHF we propose in the current study can be viewed as a matrix
generalisation of stochastic forced Hamiltonian systems studied in [23]; see
also [22, 25]. The Hamiltonian structure brings us a convenient setting for
analysis of the underlying dynamical system; moreover, it allows the
systematic construction of structure-preserving integrators possible. In this
paper, we will mainly be focused on the extension of SV type of symplectic
schemes to systems of SHF and to DPD.

The paper is organised as follows. In Section 2, we propose the SHF and derive
the DPD by specifying the Hamiltonian functions and external/dissipative
forces properly. SV type of schemes for the SHF and the DPD are constructed in
Section 3 and in particular, we will be focused on several explicit schemes
that are applied to DPD simulations in Section 4. Finally, we conclude and
point out some future researches in Section 5.

##  2 The stochastic Hamiltonian formulation with external forces

Let $Q$ be an $n$-dimensional configuration space of a mechanical system with
$\bm{q}$ the generalised coordinates. Let $(\bm{q},\dot{\bm{q}})\in TQ$ and
$(\bm{q},\bm{p})\in T^{*}Q$ be coordinates of the tangent bundle and the
cotangent bundle, respectively. We propose a stochastic Hamiltonian
formulation (SHF) with external forces as a dynamical system in $T^{*}Q$ as
follows:

$\displaystyle\left(\begin{array}[]{c}\operatorname{d}\\!{\bm{q}}\\\
\operatorname{d}\\!{\bm{p}}\end{array}\right)=J\nabla
H(\bm{q},\bm{p})\operatorname{d}\\!t$
$\displaystyle+\left(\begin{array}[]{c}0\\\
\bm{F}^{\operatorname{D}}(\bm{q},\bm{p})\end{array}\right)\operatorname{d}\\!t$
(1) $\displaystyle+\sum_{i=1}^{K}\sum_{j=1}^{K}\left(J\nabla
h_{ij}(\bm{q},\bm{p})+\left(\begin{array}[]{c}0\\\
\bm{F}^{\operatorname{SD}}_{ij}(\bm{q},\bm{p})\end{array}\right)\right)\circ{\operatorname{d}\\!W_{ij}(t)},$

where $\circ$ denotes the Stratonovich integration, $J$ is the canonical
symplectic matrix

$J=\left(\begin{array}[]{cc}0&I_{n}\\\ -I_{n}&0\end{array}\right),$ (2)

$\bm{F}^{\operatorname{D}}:T^{*}Q\rightarrow T^{*}Q$ and
$\bm{F}_{ij}^{\operatorname{SD}}:T^{*}Q\rightarrow T^{*}Q$ are fibre-
preserving maps of the external forces leading to dissipation, the functions
$H:T^{*}Q\rightarrow\mathbb{R}$ and $h_{ij}:T^{*}Q\rightarrow\mathbb{R}$ are
the Hamiltonian functions, and components of the symmetric $K\times K$ random
matrix $W(t)$ are independent Wiener processes. Note that the indices $i,j$
are not necessary of the same dimension.The superindices $\operatorname{D}$
and $\operatorname{SD}$ are shorthand for ‘Dissipation’ and ‘Stochastic
Dissipation’, respectively. For more details on stochastic differential
equations, the reader may refer to [26, 27, 28].

The SHF (1) can be written in the Itô form as

$\operatorname{d}\\!\bm{z}=A(\bm{z})\operatorname{d}\\!t+\sum_{i=1}^{K}\sum_{j=1}^{K}B_{ij}(\bm{z})\operatorname{d}\\!W_{ij}(t),$
(3)

where $\bm{z}=(\bm{q},\bm{p})^{\operatorname{T}}$,

$A(\bm{z})=\left(\begin{array}[]{c}\nabla_{\bm{p}}H+\frac{1}{2}\sum\limits_{i=1}^{K}\sum\limits_{j=1}^{K}\left(\frac{\partial^{2}h_{ij}}{\partial\bm{p}\partial\bm{q}}\left(\nabla_{\bm{p}}h_{ij}\right)+\frac{\partial^{2}h_{ij}}{\partial\bm{p}^{2}}\left(\bm{F}_{ij}^{\operatorname{SD}}-\nabla_{\bm{q}}h_{ij}\right)\right)\\\
-\nabla_{\bm{q}}H+\bm{F}^{\operatorname{D}}+\frac{1}{2}\sum\limits_{i=1}^{K}\sum\limits_{j=1}^{K}\left(\left(\frac{\partial^{2}h_{ij}}{\partial\bm{q}\partial\bm{p}}-\nabla_{\bm{\bm{p}}}\bm{F}_{ij}^{\operatorname{SD}}\right)\left(\nabla_{\bm{p}}h_{ij}-\bm{F}_{ij}^{\operatorname{SD}}\right)\right.\\\
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\left.-\left(\frac{\partial^{2}h_{ij}}{\partial\bm{q}^{2}}-\nabla_{\bm{q}}\bm{F}_{ij}^{\operatorname{SD}}\right)\left(\nabla_{\bm{p}}h_{ij}\right)\right)\end{array}\right)$
(4)

and

$B_{ij}(\bm{z})=\left(\begin{array}[]{c}\nabla_{\bm{p}}h_{ij}\\\
-\nabla_{\bm{q}}h_{ij}+\bm{F}_{ij}^{\operatorname{SD}}\end{array}\right).$ (5)

Here, ${\partial^{2}h_{ij}}/{\partial\bm{p}\partial\bm{q}}$,
${\partial^{2}h_{ij}}/{\partial\bm{q}}^{2}$ and
${\partial^{2}h_{ij}}/{\partial\bm{p}}^{2}$ denote the Hessian matrices of
$h_{ij}$, and $\nabla$ denotes the gradient of functions. Throughout the
paper, we will employ the conventional assumptions that the Hamiltonians $H$
and $h_{ij}$ are all $C^{2}$ functions and $A$ and $B_{ij}$ are globally
Lipschitz [26, 28].

######  Remark 2.1.

The SHF can be derived through variational calculus. It will be called a
stochastic Lagrange–d’Alembert principle in the phase space $T^{*}Q$, reading

$\displaystyle\delta$
$\displaystyle\int_{t_{a}}^{t_{b}}\left(\bm{p}\circ\operatorname{d}\\!{\bm{q}}-H(\bm{q},\bm{p})\operatorname{d}\\!t\right)+\int_{t_{a}}^{t_{b}}\bm{F}^{\operatorname{D}}(\bm{q},\bm{p})\cdot\delta\bm{q}\operatorname{d}\\!t$
(6)
$\displaystyle~{}~{}+\sum_{i=1}^{K}\sum_{j=1}^{K}\left(\delta\int_{t_{a}}^{t_{b}}-h_{ij}(\bm{q},\bm{p})\circ\operatorname{d}\\!W_{ij}(t)+\int_{t_{a}}^{t_{b}}\left(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p})\cdot\delta\bm{q}\right)\circ\operatorname{d}\\!W_{ij}(t)\right)=0.$

The time interval is $[t_{a},t_{b}]$ ($t_{a}<t_{b}$). The first row denotes
all deterministic terms, while the second row includes all stochastic terms.

Solutions of the SHF (1) satisfies the stochastic Lagrange–d’Alembert
principle (6); see, e.g., [23]. The converse is also true providing the
regularity of $\bm{q}$ and $\bm{p}$ [29]. In particular if
$h_{ij}=h_{ij}(\bm{q})$ are all independent of $\bm{p}$, which is exactly the
case for DPD, $(\bm{q},\bm{p})$ is a solution of the SHF (1) if and only if it
satisfies the stochastic Lagrange–d’Alembert principle (6) [30].

DPD derived from the SHF. To derive the DPD system of $N$ particles, we assume
that there exist no stochastic dissipative forces, meaning that

$\bm{F}^{\operatorname{SD}}_{ij}(\bm{q},\bm{p})\equiv 0,\quad\forall
i,j=1,2,\ldots,N.$ (7)

In the general SHF formulation (1), introduce the local coordinates for the
cotangent bundle of $N$ copies of $Q$ as

$\bm{q}=(\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{N}),\quad\bm{p}=(\bm{p}_{1},\bm{p}_{2},\ldots,\bm{p}_{N}),$
(8)

where $(\bm{q}_{i},\bm{p}_{i})$ are the coordinates of the phase space
$T^{*}Q$ for the $i$-th particle. As commonly considered in DPD, we will be
focused on the three-dimensional Euclidean space, i.e., $Q=\mathbb{R}^{3}$, in
the current study. Define the Hamiltonian $H(\bm{q},\bm{p})$ as the total
energy:

$H(\bm{q},\bm{p})=\sum_{i=1}^{N}\frac{1}{2m_{i}}|\bm{p}_{i}|^{2}+V(\bm{q}),$
(9)

where the potential energy $V(\bm{q})$ is given by

$V(\bm{q})=\sum_{i=1}^{N}\sum_{j=1}^{N}\frac{a_{ij}}{4}q_{\mathrm{c}}\left(1-\frac{q_{ij}}{q_{\mathrm{c}}}\right)^{2}\delta_{ij}.$
(10)

Here $m_{i}$ is mass of the $i$-th particle, $q_{\mathrm{c}}$ is a constant,
$a_{N\times N}$ is a constant symmetric matrix,
$q_{ij}=|\bm{q}_{i}-\bm{q}_{j}|$ is the distance of the $i$-th and the $j$-th
particles, and $\delta_{ij}$ is given by

$\delta_{ij}=\left\\{\begin{array}[]{cl}1,&q_{ij}<q_{\mathrm{c}},\vspace{0.2cm}\\\
0,&q_{ij}\geq q_{\mathrm{c}}.\end{array}\right.$ (11)

######  Remark 2.2.

Direct computation gives gradient of the Hamiltonian $H$ as follows

$\nabla H(\bm{q},\bm{p})=\left(-\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}),\frac{\bm{p}_{i}}{m_{i}}\right)^{\operatorname{T}},$
(12)

where the conservative force reads

$\bm{F}^{\operatorname{C}}_{ij}(\bm{q})=a_{ij}\left(1-\frac{q_{ij}}{q_{\mathrm{c}}}\right)\delta_{ij}\bm{\widehat{q}}_{ij},\quad
i,j=1,2,\ldots,N,\quad i\neq j,$ (13)

with

$\bm{\widehat{q}}_{ij}=\frac{\bm{q}_{i}-\bm{q}_{j}}{q_{ij}}=\frac{\bm{q}_{i}-\bm{q}_{j}}{|\bm{q}_{i}-\bm{q}_{j}|}$
(14)

and the superindex $\operatorname{C}$ meaning ‘Conservation’. Obviously, the
conservative force $\bm{F}^{\operatorname{C}}_{ij}(\bm{q})$ between the $i$-th
and the $j$-th particles only depends on their relative distance
$\bm{q}_{i}-\bm{q}_{j}$.

Furthermore, the (deterministic) dissipative force is defined by

$\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p})=-\gamma\sum_{j\neq
i}\omega^{\operatorname{D}}(q_{ij})\left(\bm{\widehat{q}}_{ij}\cdot\bm{v}_{ij}\right)\bm{\widehat{q}}_{ij},\quad
i=1,2,\ldots,N,$ (15)

where $\gamma$ is a constant friction parameter,

$\bm{v}_{ij}=\frac{\bm{p}_{i}}{m_{i}}-\frac{\bm{p}_{j}}{m_{j}}$ (16)

and
$\omega^{\operatorname{D}}(q_{ij})=\left(\omega^{\operatorname{R}}(q_{ij})\right)^{2}$
with

$\omega^{\operatorname{R}}(q_{ij})=\left(1-\frac{q_{ij}}{q_{\mathrm{c}}}\right)\delta_{ij}.$
(17)

Here, the superindex $\operatorname{R}$ means ‘Randomness’.

Let $k=N$ and define the Hamiltonian functions $h_{ij}(\bm{q},\bm{p})$
($i,j,=1,2,\ldots,N$) by

$h_{ij}(\bm{q})=\frac{\sigma}{4}q_{\mathrm{c}}\left(1-\frac{q_{ij}}{q_{\mathrm{c}}}\right)^{2}\delta_{ij},$
(18)

where $\sigma$ is a constant noise parameter. Obviously,
$h_{ii}(\bm{q})\equiv\text{const}$ for all $i=1,2,\ldots,N$ and hence $\nabla
h_{ii}(\bm{q})\equiv 0$.

######  Remark 2.3.

When $i\neq j$, since the Hamiltonian function $h_{ij}(\bm{q})$ only depends
on $\bm{q}_{i}$ and $\bm{q}_{j}$, nonzero components of its gradient are given
by

$\displaystyle\nabla_{\bm{q}_{i}}h_{ij}(\bm{q})$
$\displaystyle=-\frac{\sigma}{2}\left(1-\frac{q_{ij}}{q_{\mathrm{c}}}\right)\delta_{ij}\bm{\widehat{q}}_{ij}=-\frac{\sigma}{2}\omega^{\operatorname{R}}(q_{ij})\bm{\widehat{q}}_{ij},$
(19) $\displaystyle\nabla_{\bm{q}_{j}}h_{ij}(\bm{q})$
$\displaystyle=-\frac{\sigma}{2}\left(1-\frac{q_{ij}}{q_{\mathrm{c}}}\right)\delta_{ij}\bm{\widehat{q}}_{ji}=-\frac{\sigma}{2}\omega^{\operatorname{R}}(q_{ij})\bm{\widehat{q}}_{ji}.$

Substituting the functions specified above to the SHF (1), we obtain the
system of DPD as follows

$\left\\{\begin{aligned} \dot{\bm{q}}_{i}&=\frac{\bm{p}_{i}}{m_{i}},\\\
\dot{\bm{p}}_{i}&=\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q})+\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p})+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij})\bm{\widehat{q}}_{ij}\circ\frac{\operatorname{d}\\!W_{ij}(t)}{\operatorname{d}\\!t},\end{aligned}\right.$
(20)

for $i=1,2,\ldots,N$, in which the dissipative force
$\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p})$ is given by (15) while the
conservative force $\bm{F}^{\operatorname{C}}_{i}(\bm{q},\bm{p})$ and
randomness contribution are respectively derived from the Hamiltonians
$H(\bm{q},\bm{p})$, i.e., the total energy (9), and $h_{ij}(\bm{q},\bm{p})$
defined in (18). It is obvious that the DPD system (20) can also be obtained
via the stochastic Lagrange–d’Alembert principle (6). Note that the SHF (1) is
formally divided by $\operatorname{d}\\!t$ on both sides to obtain the system
(20), which has been the conventional form of DPD.

######  Remark 2.4.

Since no stochastic dissipative forces exist and $h_{ij}=h_{ij}(\bm{q})$ are
independent from $\bm{p}$, SHF’s Itô form (3), in particular the coefficient
matrix $A(\bm{z})$ given by (4), yields that the DPD (20) takes the same form
in both the Itô framework and the Stratonovich framework.

##  3 Störmer–Verlet schemes for the SHF and the DPD

In this section, we propose the Störmer–Verlet (SV) type of symplectic schemes
for the DPD based on the SHF (1). That is, when no external forces and
randomness are imposed, the corresponding discrete ‘flow’ shall be symplectic
as well. In other words, dissipation in the numerical schemes is only
contributed by the external forces, same as what occurs in the continuous
counterpart.

###  3.1 SV type of schemes for the SHF

Discretize the time interval $[t_{a},t_{b}]$ as a series
$t_{a}=t_{0},t_{1},t_{2},\ldots,t_{K}=t_{b}$ and denote

$\Delta t=t_{k+1}-t_{k}=\frac{t_{b}-t_{a}}{K}$

as the time step. The space $TT^{*}Q$ where SHF systems (and the corresponding
variational structure) are defined is discretized into two copies of the
cotangent bundle, i.e., $T^{*}Q\times T^{*}Q$, with local coordinates
$(\bm{q}^{k},\bm{p}^{k},\bm{q}^{k+1},\bm{p}^{k+1})$ where
$\bm{q}^{k}=\bm{q}(t_{k})$, $\bm{p}^{k}=\bm{p}(t_{k})$ and so forth. In the
current paper, we will mainly be focused on extensions of the SV schemes for
systems of SHF (1), which are symplectic schemes of second order accuracy for
conservative Hamiltonian systems. The SV schemes arise as the composite of
Euler-A and Euler-B methods which are both symplectic, implicit and of first
order accuracy for conservative Hamiltonian systems. We will follow a similar
approach to introduce SV schemes for the SHF.

For SHF (1), we propose a family of Euler-A methods:

$\displaystyle{\bm{q}^{k+1}-\bm{q}^{k}}$ $\displaystyle=\Delta
t*\nabla_{\bm{p}}H(\bm{q}^{k+1},\bm{p}^{k})+\sum_{i,j}\nabla_{\bm{p}}h_{ij}(\bm{q}^{k+1},\bm{p}^{k})\circ\Delta
W_{ij}(t_{k}),$ (21) $\displaystyle{\bm{p}^{k+1}-\bm{p}^{k}}$
$\displaystyle=-\Delta
t\left(\nabla_{\bm{q}}H(\bm{q}^{k+1},\bm{p}^{k})+(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k}\right)$
$\displaystyle~{}~{}~{}~{}+\sum_{i,j}\left(-\nabla_{\bm{q}}h_{ij}(\bm{q}^{k+1},\bm{p}^{k})+(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k}\right)\circ\Delta
W_{ij}(t_{k}),$

and a family of Euler-B methods:

$\displaystyle{\bm{q}^{k+1}-\bm{q}^{k}}$ $\displaystyle=\Delta
t*\nabla_{\bm{p}}H(\bm{q}^{k},\bm{p}^{k+1})+\sum_{i,j}\nabla_{\bm{p}}h_{ij}(\bm{q}^{k},\bm{p}^{k+1})\circ\Delta
W_{ij}(t_{k}),$ (22) $\displaystyle{\bm{p}^{k+1}-\bm{p}^{k}}$
$\displaystyle=-\Delta
t\left(\nabla_{\bm{q}}H(\bm{q}^{k},\bm{p}^{k+1})+(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k}\right)$
$\displaystyle~{}~{}~{}~{}+\sum_{i,j}\left(-\nabla_{\bm{q}}h_{ij}(\bm{q}^{k},\bm{p}^{k+1})+(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k}\right)\circ\Delta
W_{ij}(t_{k}),$

where $(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k}$ and
$(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k}$ denote discretisations
of the external forces, and

$\Delta W_{ij}(t_{k})=W_{ij}(t_{k+1})-W_{ij}(t_{k})\sim\mathcal{N}(0,\Delta
t).$ (23)

Here $\mathcal{N}(0,\Delta t)$ denotes the normal distribution with mean $0$
and standard deviation $\sqrt{\Delta t}$.

Two types of SV schemes can be defined as composites of the Euler methods with
time step $\Delta t/2$, namely $\text{(Euler-A)}\circ\text{(Euler-B)}$ and
$\text{(Euler-B)}\circ\text{(Euler-A)}$, which will be called SV-AB schemes
and SV-BA schemes, respectively.

The family of SV-AB schemes, namely $\text{(Euler-A)}\circ\text{(Euler-B)}$,
reads

$\displaystyle\bm{p}^{k+1/2}$ $\displaystyle-\bm{p}^{k}=\frac{\Delta
t}{2}\left[-\nabla_{\bm{q}}H(\bm{q}^{k},\bm{p}^{k+1/2})+(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{1}}\right]$
(24)
$\displaystyle\quad\quad\quad+\sum_{i,j}\left(-\nabla_{\bm{q}}h_{ij}(\bm{q}^{k},\bm{p}^{k+1/2})+(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k_{1}}\right)\circ\overline{\Delta}W_{ij}(t_{k}),$
$\displaystyle\bm{q}^{k+1}$ $\displaystyle-\bm{q}^{k}=\frac{\Delta
t}{2}\left[\nabla_{\bm{p}}H(\bm{q}^{k},\bm{p}^{k+1/2})+\nabla_{\bm{p}}H(\bm{q}^{k+1},\bm{p}^{k+1/2})\right]$
$\displaystyle+\sum_{i,j}\nabla_{\bm{p}}h_{ij}(\bm{q}^{k},\bm{p}^{k+1/2})\circ\overline{\Delta}W_{ij}(t_{k})+\sum_{i,j}\nabla_{\bm{p}}h_{ij}(\bm{q}^{k+1},\bm{p}^{k+1/2})\circ\overline{\Delta}W_{ij}(t_{k+1/2}),$
$\displaystyle\bm{p}^{k+1}$ $\displaystyle-\bm{p}^{k+1/2}=\frac{\Delta
t}{2}\left[-\nabla_{\bm{q}}H(\bm{q}^{k+1},\bm{p}^{k+1/2})+(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{2}}\right]$
$\displaystyle\quad\quad\quad+\sum_{i,j}\left(-\nabla_{\bm{q}}h_{ij}(\bm{q}^{k+1},\bm{p}^{k+1/2})+(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k_{2}}\right)\circ\overline{\Delta}W_{ij}(t_{k+1/2}),$

while the family of SV-BA schemes, namely
$\text{(Euler-B)}\circ\text{(Euler-A)}$, reads

$\displaystyle\bm{q}^{k+1/2}-\bm{q}^{k}$ $\displaystyle=\frac{\Delta
t}{2}*\nabla_{\bm{p}}H(\bm{q}^{k+1/2},\bm{p}^{k})+\sum_{i,j}\nabla_{\bm{p}}h_{ij}(\bm{q}^{k+1/2},\bm{p}^{k})\circ\overline{\Delta}W_{ij}(t_{k}),$
(25) $\displaystyle\bm{p}^{k+1}-\bm{p}^{k}$ $\displaystyle=\frac{\Delta
t}{2}\left[-\nabla_{\bm{q}}H(\bm{q}^{k+1/2},\bm{p}^{k})-\nabla_{\bm{q}}H(\bm{q}^{k+1/2},\bm{p}^{k+1})\right]+\Delta
t*(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k}$
$\displaystyle~{}~{}~{}~{}+\sum_{i,j}\left(-\nabla_{\bm{q}}h_{ij}(\bm{q}^{k+1/2},\bm{p}^{k})+(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k_{1}}\right)\circ\overline{\Delta}W_{ij}(t_{k})$
$\displaystyle~{}~{}~{}~{}+\sum_{i,j}\left(-\nabla_{\bm{q}}h_{ij}(\bm{q}^{k+1/2},\bm{p}^{k+1})+(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k_{2}}\right)\circ\overline{\Delta}W_{ij}(t_{k+1/2}),$
$\displaystyle\bm{q}^{k+1}-\bm{q}^{k+1/2}$ $\displaystyle=\frac{\Delta
t}{2}*\nabla_{\bm{p}}H(\bm{q}^{k+1/2},\bm{p}^{k+1})+\sum_{i,j}\nabla_{\bm{p}}h_{ij}(\bm{q}^{k+1/2},\bm{p}^{k+1})\circ\overline{\Delta}W_{ij}(t_{k+1/2}),$

where $(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k}$,
$(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{1}}$ and
$(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{2}}$ denote three independent
discretisations of the force $\bm{F}^{\operatorname{D}}(\bm{q},\bm{p})$,
$(\bm{F}^{\operatorname{SD}}_{ij}(\bm{q},\bm{p}))^{k_{1}}$ and
$(\bm{F}^{\operatorname{SD}}_{ij}(\bm{q},\bm{p}))^{k_{2}}$ denote two
independent discretisations of the force
$\bm{F}^{\operatorname{SD}}_{ij}(\bm{q},\bm{p})$, and

$\overline{\Delta}W_{ij}(t_{k})=W_{ij}(t_{k+1/2})-W_{ij}(t_{k})\sim\mathcal{N}(0,\Delta
t/2).$ (26)

######  Remark 3.5.

It is obvious that the SV schemes (24) and (25) reduce to the ordinary SV
schemes for conservative Hamiltonian systems, assuming the absence of external
forces and stochastic terms. Consequently, discretisations of the external
forces can, in principle, be chosen arbitrarily, providing the resulting
schemes are stable and convergent. Only when discretisations of the external
forces are specified properly, they are a 2-stage stochastic partitioned
Runge–Kutta method given in [23]; however, in DPD simulations, for instance
the GW and GCC methods, these discretisations are often chosen very
differently as we will find out below.

Separable Hamiltonians. Assuming the Hamiltonians can be separated as

$H(\bm{q},\bm{p})=T(\bm{p})+V(\bm{q})\text{ and
}h_{ij}(\bm{q},\bm{p})=S_{ij}(\bm{p})+U_{ij}(\bm{q}),$ (27)

the SV-AB schemes (24) and SV-BA schemes (25) become

$\displaystyle\bm{p}^{k+1/2}-\bm{p}^{k}$ $\displaystyle=\frac{\Delta
t}{2}\left[-\nabla_{\bm{q}}V(\bm{q}^{k})+(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{1}}\right]$
(28)
$\displaystyle\quad\quad+\sum_{i,j}\left(-\nabla_{\bm{q}}U_{ij}(\bm{q}^{k})+(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k_{1}}\right)\circ\overline{\Delta}W_{ij}(t_{k}),$
$\displaystyle\bm{q}^{k+1}-\bm{q}^{k}$ $\displaystyle=\Delta
t*\nabla_{\bm{p}}T(\bm{p}^{k+1/2})+\sum_{i,j}\nabla_{\bm{p}}S_{ij}(\bm{p}^{k+1/2})\circ\Delta
W_{ij}(t_{k}),$ $\displaystyle\bm{p}^{k+1}-\bm{p}^{k+1/2}$
$\displaystyle=\frac{\Delta
t}{2}\left[-\nabla_{\bm{q}}V(\bm{q}^{k+1})+(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{2}}\right]$
$\displaystyle\quad\quad+\sum_{i,j}\left(-\nabla_{\bm{q}}U_{ij}(\bm{q}^{k+1})+(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k_{2}}\right)\circ\overline{\Delta}W_{ij}(t_{k+1/2}),$

and

$\displaystyle\bm{q}^{k+1/2}-\bm{q}^{k}$ $\displaystyle=\frac{\Delta
t}{2}*\nabla_{\bm{p}}T(\bm{p}^{k})+\sum_{i,j}\nabla_{\bm{p}}S_{ij}(\bm{p}^{k})\circ\overline{\Delta}W_{ij}(t_{k}),$
(29) $\displaystyle\bm{p}^{k+1}-\bm{p}^{k}$ $\displaystyle=\Delta
t*\left(-\nabla_{\bm{q}}V(\bm{q}^{k+1/2})+(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k}\right)-\sum_{ij}\nabla_{\bm{q}}U_{ij}(\bm{q}^{k+1/2})\circ\Delta
W_{ij}(t_{k})$
$\displaystyle~{}~{}~{}~{}+\sum_{i,j}\left((\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k_{1}}\circ\overline{\Delta}W_{ij}(t_{k})+(\bm{F}_{ij}^{\operatorname{SD}}(\bm{q},\bm{p}))^{k_{2}}\circ\overline{\Delta}W_{ij}(t_{k+1/2})\right),$
$\displaystyle\bm{q}^{k+1}-\bm{q}^{k+1/2}$ $\displaystyle=\frac{\Delta
t}{2}*\nabla_{\bm{p}}T(\bm{p}^{k+1})+\sum_{i,j}\nabla_{\bm{p}}S_{ij}(\bm{p}^{k+1})\circ\overline{\Delta}W_{ij}(t_{k+1/2}).$

###  3.2 SV schemes for the DPD

The Hamiltonians (9) and (18) of the DPD (20) can obviously be separated with
respect to their position and momentum coordinates. Substituting them into
(28) and (29), the SV-AB schemes and SV-BA schemes for DPD turn out to be

$\displaystyle{\bm{p}^{k+1/2}_{i}-\bm{p}^{k}_{i}}$ $\displaystyle=\frac{\Delta
t}{2}\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{1}}\right]+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\overline{\Delta}W_{ij}(t_{k}),$
(30) $\displaystyle{\bm{q}^{k+1}_{i}-\bm{q}^{k}_{i}}$ $\displaystyle=\Delta
t*\frac{\bm{p}_{i}^{k+1/2}}{m_{i}},$
$\displaystyle{\bm{p}^{k+1}_{i}-\bm{p}^{k+1/2}_{i}}$
$\displaystyle=\frac{\Delta t}{2}\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k+1})+(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{2}}\right]+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k+1})\bm{\widehat{q}}_{ij}^{k+1}\circ\overline{\Delta}W_{ij}(t_{k+1/2}),$

and

$\displaystyle\bm{q}^{k+1/2}_{i}-\bm{q}^{k}_{i}$ $\displaystyle=\frac{\Delta
t}{2}*\frac{\bm{p}_{i}^{k}}{m_{i}},$ (31)
$\displaystyle\bm{p}^{k+1}_{i}-\bm{p}^{k}_{i}$ $\displaystyle=\Delta
t\left(\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k+1/2})+(\bm{F}_{i}^{\operatorname{D}}(\bm{q},\bm{p}))^{k}\right)+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k+1/2})\bm{\widehat{q}}_{ij}^{k+1/2}\circ\Delta
W_{ij}(t_{k}),$ $\displaystyle\bm{q}^{k+1}_{i}-\bm{q}^{k+1/2}_{i}$
$\displaystyle=\frac{\Delta t}{2}*\frac{\bm{p}_{i}^{k+1}}{m_{i}}.$

######  Remark 3.6.

If we (partially) eliminate the half values in the SV-AB schemes (30) and SV-
BA schemes (31), we can rewrite them in the following equivalent
representatives

$\displaystyle\bm{q}^{k+1}_{i}-\bm{q}^{k}_{i}$ $\displaystyle=\frac{\Delta
t}{m_{i}}\left\\{\bm{p}_{i}^{k}+\frac{\Delta t}{2}\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{1}}\right]+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\overline{\Delta}W_{ij}(t_{k})\right\\},$
(32) $\displaystyle\bm{p}^{k+1}_{i}-\bm{p}^{k}_{i}$
$\displaystyle=\frac{\Delta t}{2}\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{1}}\right]+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\overline{\Delta}W_{ij}(t_{k})$
$\displaystyle\quad\quad+\frac{\Delta t}{2}\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k+1})+(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{2}}\right]+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k+1})\bm{\widehat{q}}_{ij}^{k+1}\circ\overline{\Delta}W_{ij}(t_{k+1/2})$

and

$\displaystyle\bm{q}^{k+1/2}_{i}-\bm{q}^{k}_{i}$ $\displaystyle=\frac{\Delta
t}{2}*\frac{\bm{p}_{i}^{k}}{m_{i}},$ (33)
$\displaystyle\bm{p}^{k+1}_{i}-\bm{p}^{k}_{i}$ $\displaystyle=\Delta
t\left(\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k+1/2})+(\bm{F}_{i}^{\operatorname{D}}(\bm{q},\bm{p}))^{k}\right)+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k+1/2})\bm{\widehat{q}}_{ij}^{k+1/2}\circ\Delta
W_{ij}(t_{k}),$ $\displaystyle\bm{q}^{k+1}_{i}-\bm{q}^{k}_{i}$
$\displaystyle=\frac{\Delta
t}{m_{i}}*\frac{\bm{p}_{i}^{k}+\bm{p}_{i}^{k+1}}{2}.$

Note that in the latter, $\bm{q}^{k+1/2}$ can also be totally eliminated. We
keep it to avoid heavy arguments for the functions.

In the rest of the paper, we will be focused on the SV-AB schemes (30) (or
(32)) for DPD, which include the GW and GCC methods as special cases. Further
studies on the SV-BA schemes and other symplectic methods will be conducted in
our future work. We need only specify the force discretisations
$(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{1}}$ and
$(\bm{F}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{2}}$, respectively. There are
certainly many other choices expect for what we introduce below.

To recover the conventional GW and GCC methods, the approximation
$\overline{\Delta}W_{ij}(t_{k})\approx\overline{\Delta}W_{ij}(t_{k+1/2})$ will
have to be employed, and hence

$\overline{\Delta}W_{ij}\approx\Delta W_{ij}/2,$ (34)

as $\overline{\Delta}W_{ij}(t_{k})+\overline{\Delta}W_{ij}(t_{k+1/2})=\Delta
W_{ij}(t_{k})$. However, it should be noted that this approximation will
change the nature of the schemes in the sense that the increments
$\overline{\Delta}W_{ij}(t_{k})$ and $\overline{\Delta}W_{ij}(t_{k+1/2})$ are
not longer independent; in fact, this approximation is not necessary in
practical applications.

  * 1.

SV-AB-1 is an implicit scheme by choosing

$(\bm{F}_{i}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{1}}=\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\bm{p}^{k}),\quad(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{2}}=\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k+1},\bm{p}^{k+1}).$
(35)

For DPD, the dissipative force $\bm{F}^{\operatorname{D}}$ is linear in
$\bm{p}$, so the scheme can be written explicitly, in principle. However, one
may need to solve a linear system with a sparse coefficient matrix.

  * 2.

SV-AB-2 is an explicit scheme by defining

$(\bm{F}_{i}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{1}}=\bm{F}_{i}^{\operatorname{D}}(\bm{q}^{k},\bm{p}^{k}),\quad(\bm{F}_{i}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{2}}=\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k+1},\bm{p}^{k+\lambda}),$
(36)

where $\bm{p}^{k+\lambda}$ ($\lambda\in[0,1]$) is defined by

$\frac{\bm{p}^{k+\lambda}_{i}-\bm{p}^{k}_{i}}{\Delta
t}=\lambda\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\bm{p}^{k})+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\frac{\Delta
W_{ij}(t_{k})}{\Delta t}\right].$ (37)

This is exactly the GCC method [14].

  * 3.

SV-AB-3 is explicit by specifying

$(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{1}}=\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\bm{p}^{k-1+\lambda}),\quad(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{2}}=\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k+1},\bm{p}^{k+\lambda}),$
(38)

where $\bm{p}^{k+\lambda}$ ($\lambda\in[0,1]$) is defined by

$\frac{\bm{p}^{k+\lambda}_{i}-\bm{p}^{k}_{i}}{\Delta
t}=\lambda\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\bm{p}^{k-1+\lambda})+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\frac{\Delta
W_{ij}(t_{k})}{\Delta t}\right]$ (39)

and the initial value of $\bm{p}^{k+\lambda}$ is $\bm{p}^{\lambda}=\bm{p}^{1}$
when $k=1$. This is exactly the GW method [2].

  * 4.

SV-AB-4 is a generalisation of the three methods above, which can, in
principle, be expressed explicitly for the DPD:

$\displaystyle(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{1}}$
$\displaystyle=\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\alpha\bm{p}^{k}+(1-\alpha)\bm{p}^{k-1+\lambda}),\quad\alpha\in[0,1],$
(40) $\displaystyle(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{2}}$
$\displaystyle=\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k+1},\beta\bm{p}^{k+\lambda}+(1-\beta)\bm{p}^{k+1}),\quad\beta\in[0,1],$

where $\bm{p}^{k+\lambda}$ ($\lambda\in[0,1]$) is defined by

$\displaystyle\frac{\bm{p}^{k+\lambda}_{i}-\bm{p}^{k}_{i}}{\Delta t}$
$\displaystyle=\lambda\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\alpha\bm{p}^{k}+(1-\alpha)\bm{p}^{k-1+\lambda})\right.$
(41)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\left.+~{}\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\frac{\Delta
W_{ij}(t_{k})}{{\Delta t}}\right].$

It reduces to the SV-AB-1 method for $\alpha=1,\beta=0$, to the SV-AB-2 (GCC)
method for $\alpha=1,\beta=1$ and to the SV-AB-3 (GW) method for
$\alpha=0,\beta=1$.

  * 5.

SV-AB-5 is explicit by choosing

$\displaystyle(\bm{F}_{i}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{1}}$
$\displaystyle=\bm{F}_{i}^{\operatorname{D}}(\bm{q}^{k},\bm{p}^{k-1+\lambda_{1}}),\quad\lambda_{1}\in[0,1],$
(42) $\displaystyle(\bm{F}_{i}^{\operatorname{D}}(\bm{q},\bm{p}))^{k_{2}}$
$\displaystyle=\bm{F}_{i}^{\operatorname{D}}(\bm{q}^{k+1},\bm{p}^{k+\lambda_{2}}),\quad\lambda_{2}\in[0,1],$

where

$\displaystyle\frac{\bm{p}^{k+\lambda_{1}}_{i}-\bm{p}^{k}_{i}}{\Delta t}$
$\displaystyle=\lambda_{1}\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\bm{p}^{k-1+\lambda_{1}})+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\frac{\Delta
W_{ij}(t_{k})}{\Delta t}\right],$ (43)
$\displaystyle\frac{\bm{p}^{k+\lambda_{2}}_{i}-\bm{p}^{k}_{i}}{\Delta t}$
$\displaystyle=\lambda_{2}\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\bm{p}^{k-1+\lambda_{1}})+\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\frac{\Delta
W_{ij}(t_{k})}{\Delta t}\right].$

When $\lambda_{1}=\lambda_{2}$, it reduces to the SV-AB-3 (GW) method.

  * 6.

SV-AB-6 is a simultaneous generalisation of SV-AB-4 and SV-AB-5, which can be
written in an explicit form for the DPD:

$\displaystyle(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{1}}$
$\displaystyle=\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\alpha\bm{p}^{k}+(1-\alpha)\bm{p}^{k-1+\lambda_{1}}),\quad\alpha\in[0,1],\lambda_{1}\in[0,1],$
(44) $\displaystyle(\bm{F}^{\operatorname{D}}_{i}(\bm{q},\bm{p}))^{k_{2}}$
$\displaystyle=\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k+1},\beta\bm{p}^{k+\lambda_{2}}+(1-\beta)\bm{p}^{k+1}),\quad\beta\in[0,1],\lambda_{2}\in[0,1],$

where $\bm{q}^{k+\lambda_{1}}$ and $\bm{q}^{k+\lambda_{2}}$ are given by

$\displaystyle\frac{\bm{p}^{k+\lambda_{1}}_{i}-\bm{p}^{k}_{i}}{\Delta t}$
$\displaystyle=\lambda_{1}\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\alpha\bm{p}^{k}+(1-\alpha)\bm{p}^{k-1+\lambda_{1}})\right.$
(45)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\left.+~{}\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\frac{\Delta
W_{ij}(t_{k})}{\Delta t}\right],$
$\displaystyle\frac{\bm{p}^{k+\lambda_{2}}_{i}-\bm{p}^{k}_{i}}{\Delta t}$
$\displaystyle=\lambda_{2}\left[\sum_{j\neq
i}\bm{F}_{ij}^{\operatorname{C}}(\bm{q}^{k})+\bm{F}^{\operatorname{D}}_{i}(\bm{q}^{k},\alpha\bm{p}^{k}+(1-\alpha)\bm{p}^{k-1+\lambda_{1}})\right.$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\left.+~{}\sigma\sum_{j\neq
i}\omega^{\operatorname{R}}(q_{ij}^{k})\bm{\widehat{q}}_{ij}^{k}\circ\frac{\Delta
W_{ij}(t_{k})}{\Delta t}\right].$

When $\lambda_{1}=\lambda_{2}$, it becomes SV-AB-4, while when
$\alpha=0,\beta=1$, it becomes SV-AB-5.

######  Remark 3.7.

The schemes SV-AB-1$\sim$6 are related through the following diagram:

SV-AB-1SV-AB-2 (GCC)SV-AB-3 (GW)SV-AB-4SV-AB-5SV-
AB-6$\alpha=1$$\beta=1$$\alpha=1$$\beta=0$$\alpha=0$$\beta=1$$\lambda_{1}=\lambda_{2}$$\lambda_{1}=\lambda_{2}$$\alpha=0$$\beta=1$

Figure 1: Relations of the schemes SV-AB-1$\sim$6.

We summarize some general features of the schemes in Fig. 1 as follows.

  * 1.

Explicit schemes: SV-AB-2 (GCC), SV-AB-3 (GW), SV-AB-5. The other schemes are
implicit but can be written explicitly for the DPD by solving a sparse linear
system.

  * 2.

Number of independent parameters: SV-AB-1: 0, SV-AB-2 (GCC): 1, SV-AB-3 (GW):
1, SV-AB-4: 3, SV-AB-5: 2, SV-AB-6: 4.

###  3.3 Long-time behaviour of the SV-AB methods

When no external forces are involved, it was noticed that the Euler-A method
(21) and the Euler-B method (22) are not convergent in the mean-square sense
when the Hamiltonian functions $h_{ij}=h_{ij}(\bm{q},\bm{p})$ depend on both
the positions and momenta [22]. If we further assume that
$h_{ij}=h_{ij}(\bm{q})$ only depend on the positions, the Euler-A and Euler-B
methods are both convergent and hence are the SV methods.

In this subsection, we will numerically show the convergence of the SV-
AB-1$\sim$6 methods by studying the damped Kubo oscillator, a stochastic
Hamiltonian system whose Hamiltonians are separable given by

$H(q,p)=\frac{p^{2}}{2}+\frac{q^{2}}{2},\quad
h(q,p)=\sigma\left(\frac{p^{2}}{2}+\frac{q^{2}}{2}\right).$ (46)

Here $\sigma$ is the noise intensity. As its solution can be calculated
analytically, it has often been used for the validation of numerical methods
(e.g., [22, 20]). By employing the forces

$F^{\operatorname{D}}=-\varepsilon p,\quad
F^{\operatorname{SD}}=-\varepsilon\sigma p$ (47)

with $\varepsilon$ the nonnegative damping coefficient, the damped Kubo
oscillator has the following exact solution [23]

$\displaystyle\overline{q}(t)$
$\displaystyle=q_{0}\exp\left(-\frac{\varepsilon}{2}(t+\sigma
W(t))\right)\cos\omega\left(t+\sigma W(t)\right)$ (48)
$\displaystyle\quad+\frac{1}{\omega}\left(p_{0}+\frac{\varepsilon}{2}q_{0}\right)\exp\left(-\frac{\varepsilon}{2}(t+\sigma
W(t))\right)\sin\omega\left(t+\sigma W(t)\right),$
$\displaystyle\overline{p}(t)$
$\displaystyle=p_{0}\exp\left(-\frac{\varepsilon}{2}(t+\sigma
W(t))\right)\cos\omega\left(t+\sigma W(t)\right)$
$\displaystyle\quad-\frac{1}{\omega}\left(q_{0}+\frac{\varepsilon}{2}p_{0}\right)\exp\left(-\frac{\varepsilon}{2}(t+\sigma
W(t))\right)\sin\omega\left(t+\sigma W(t)\right),$

where $(q_{0},p_{0})$ are the initial conditions, the angular frequency is
$\omega={\sqrt{4-\varepsilon^{2}}}/{2}$ by assuming $\varepsilon<2$. The
expected value of the Hamiltonian $H$ is given by

$\displaystyle E(H(\overline{q}(t),\overline{p}(t)))$
$\displaystyle=a\exp\left(-\frac{\varepsilon(2-\varepsilon\sigma^{2})t}{2}\right)$
(49)
$\displaystyle\quad+\exp\left(-\left((2-\varepsilon^{2})\sigma^{2}+\varepsilon\right)t\right)\left(b\cos\left(2(1-\varepsilon\sigma^{2})\omega
t\right)+c\sin\left(2(1-\varepsilon\sigma^{2})\omega t\right)\right),$

where

$a=\frac{2(q_{0}^{2}+p_{0}^{2}+\varepsilon
q_{0}p_{0})}{4-\varepsilon^{2}},\quad
b=-\frac{\varepsilon^{2}(q_{0}^{2}+p_{0}^{2})+4\varepsilon
q_{0}p_{0}}{2(4-\varepsilon^{2})},\quad
c=\frac{\varepsilon(q_{0}^{2}-p_{0}^{2})}{2\sqrt{4-\varepsilon^{2}}}.$ (50)

In the simulations, the initial conditions are $q_{0}=0,p_{0}=1$, the noise
intensity is $\sigma=0.2$ and the damping coefficient is $\varepsilon=0.001$.
Time step is $\Delta t=0.1$ for a time span $[0,2000]$. For simplicity,
discretisation of $F^{\operatorname{SD}}(q,p)$ is chosen as
$F^{\operatorname{SD}}(q^{k},p^{k})$ at each step $k$ for all numerical
methods. Furthermore, we pick one special choice of the parameters
$\alpha,\beta,\lambda$ for each method as shown in the figures and in each
case $2,000$ sample paths are generated. Figs. 2 and 3 show the mean
Hamiltonians of the SV-AB methods and their differences with respect to the
exact Hamiltonian (49). Fluctuating behaviour of the energy can be noticed. In
particular, Fig. 3 shows that order of the error is approximately $10^{-3}$
and it tends to become smaller in a long time after a relatively stronger
vibration at the beginning.

Figure 2: Mean Hamiltonians of the SV-AB-1, SV-AB-2 ($\lambda=0.7$), SV-AB-3
($\lambda=0.3$), SV-AB-4 ($\alpha=0.5,\beta=1,\lambda=0.6$), SV-AB-5
($\lambda_{1}=0.3,\lambda_{2}=0.5$), and SV-AB-6
($\lambda_{1}=0.3,\lambda_{2}=0.4,\alpha=0.4,\beta=1$) methods. To clearly
show the tendency of the time evolution and the fluctuating behaviour, the
first 200 seconds are plotted here. Figure 3: The difference between numerical
mean Hamiltonians and the exact Hamiltonian for the SV-AB-1, SV-AB-2
($\lambda=0.7$), SV-AB-3 ($\lambda=0.3$), SV-AB-4
($\alpha=0.5,\beta=1,\lambda=0.6$), SV-AB-5
($\lambda_{1}=0.3,\lambda_{2}=0.5$), and SV-AB-6
($\lambda_{1}=0.3,\lambda_{2}=0.4,\alpha=0.4,\beta=1$) methods.

##  4 Applications to DPD simulations

Although all the SV-AB schemes proposed above can be made explicit for the
DPD, further efforts may be needed to achieve the corresponding explicit
representatives, in particular, by solving a huge sparse linear system. For
simplicity, we will be focused on the explicit SV-AB-2 (GCC) and SV-AB-4
($\beta=1$) methods in comparison with the SV-AB-3 (GW) method. Recall that
SV-AB-4 ($\beta=1$) reduces to the SV-AB-3 (GW) with $\alpha=0$ and reduces to
SV-AB-2 (GCC) with $\alpha=1$ (see Fig. 1).

In our simulations, the total number of fluid particles of the same mass $m$
is set to $3,000$ with $a=25k_{\mathrm{B}}T^{*}$, where $a$ is the repulsive
parameter (i.e., $a_{ij}=a$ for all $i\neq j$) to determine the magnitude of
the conservative force $\bm{F}^{\operatorname{C}}$, $T^{*}$ is the set
temperature and $k_{\mathrm{B}}$ is the Boltzmann constant. The noise
parameter $\sigma$ and the friction parameter $\gamma$ are set to $3.0$ and
$4.5$, respectively. All simulations are performed under the condition of
constant-volume and constant-temperature, i.e., the canonical ensemble is
generated. The size of simulation box is 10 $\times 10\times
10q_{\mathrm{c}}^{3}$. The periodic boundary condition is applied in all three
dimensions. Here, $q_{\mathrm{c}}$ is the cutoff distance, which is the unit
length in the DPD simulation. The initial configuration is random, and the
initial momentum is set appropriately so that the temperature would satisfy
the Boltzmann distribution for the set temperature satisfying
$k_{\mathrm{B}}T^{*}=1.0$. This gives the repulsion parameter $a=25$, yielding
the compressibility of water. Although Groot and Warren reported that there
was no statistical difference between simulations using uniform random numbers
and those using Gaussian random numbers [2], we use a Gaussian distribution to
generate the random numbers in the current simulations.

We examined twenty cases with the time step size $\Delta t$ ranging from
$0.001$ to $0.16\tau$. Here, we use reduced units for the cutoff radius
$q_{\mathrm{c}}$, the particle mass $m$, and the energy $k_{\mathrm{B}}T$.
Hence, the time unit is defined as
$\tau=\sqrt{mq_{\mathrm{c}}^{2}/k_{\mathrm{B}}T}$. All cases were simulated
for at least $1,000\tau$, and the last $16\%$ were used as statistical data.
Note that we were not able to calculate exactly $1,000\tau$ for all $\Delta t$
and the first 84 % of the data was discarded to equilibrate the system
sufficiently. As a comparison of the accuracy of the formulations, the kinetic
temperature $k_{\mathrm{B}}T=\left\langle\bm{v}^{2}\right\rangle/3$ was
calculated and its difference from the set temperature
$k_{\mathrm{B}}T^{*}=1.0$ was examined, where $\langle\cdot\rangle$ is the
average over all particles in the simulations and $\bm{v}=\bm{p}/m$. Since the
simulation was performed with a canonical ensemble, temperature of the system
will fluctuate around a certain average value after reaching the equilibrium
state. In the simulations, the average value is the set temperature, which
satisfies $k_{\mathrm{B}}T^{*}=1.0$.

Fig. 4 plots the artificial kinetic temperature increase of the SV-AB-2 (GCC),
SV-AB-3 (GW), and SV-AB-4 ($\beta=1$) schemes with representative parameters.
For results for all parameters, please refer to Figs. S1–S3 in Supporting
Information. It is confirmed that the statistical error of the temperature is
less than $1\%$, i.e., $k_{\mathrm{B}}T-1<10^{-2}$, for all schemes when
$\Delta t$ is less than $0.01$.

Figure 4: Kinetic temperature versus time step. Curves represent
representative results for the SV-AB-2 (GCC), SV-AB-3 (GW), and SV-AB-4
($\beta=1$) schemes. Note that the kinetic temperature is averaged over time
after equilibration.

Let us firstly compare the SV-AB-3 (GW) scheme with the SV-AB-2 (GCC) scheme.
We consider that $\lambda=0.5$ and $0.65$ are the representative parameters of
the SV-AB-2 (GCC) and SV-AB-3 (GW) schemes, respectively. When $\Delta
t<2\times 10^{-2}$, in several cases the error of SV-AB-2 (GCC)
($\lambda=0.5$) is smaller than that of the SV-AB-3 (GW) with $\lambda=0.65$.
When $\Delta t>3\times 10^{-2}$, error of SV-AB-2 (GCC) ($\lambda=0.5$) jumps
to bigger than $0.1\%$. However, when the time step becomes ever bigger, for
instance $\Delta t>10^{-1}$, error of SV-AB-2 (GCC) ($\lambda=0.5$) is
smaller; one should be noted that error for these cases is probably too big
for practical simulations.

Now consider the SV-AB-4 ($\beta=1$) scheme. For all $\alpha$s, the error
tends to be the smallest around $\lambda=0.6$. When $\Delta t<10^{-2}$, the
error of the SV-AB-4 ($\beta=1$) is similar to that of the SV-AB-3 (GW)
($\lambda=0.65$). As $\Delta t$ increases, the error also increases. The
accuracy of schemes with $\alpha=0.8$ and $\alpha=0.5$ interchanges at some
point as $\Delta t$ increases: for smaller $\Delta t$, the error of
$\alpha=0.5$ case is smaller, while for larger $\Delta t$, the error of
$\alpha=0.8$ case becomes smaller. The maximum $\Delta t$ that shows an
accuracy of less than $1\%$ error is $0.06$, which is the same as the SV-AB-3
(GW) ($\lambda=0.65$), but for $\Delta t=0.04$, its accuracy is higher than
the SV-AB-3 (GW) ($\lambda=0.65$). On the other hand, when $\Delta t<7\times
10^{-2}$ and $\lambda=1.0$, the error of SV-AB-4 ($\beta=1$) is larger than
that of the SV-AB-3 (GW) ($\lambda=0.65$) for all $\alpha$s. However, when
$\Delta t=0.1$, the error is approximately $0.5\%$, which is highly accurate.
Unfortunately, further studies are needed before this can be applied in
practical simulations easily.

Simulations of the SV-AB-4 ($\beta=1$) for $\alpha=0.9$ and $\lambda=1.0$ are
shown in Fig. 5 with the vertical axis illustrated in linear scale. Blue and
red curves show the error and the absolute error respectively. Note that the
green curve in Fig. 4 and the red curve in Fig. 5 coincide. As shown in the
zoomed-up view inserted in Fig. 5, starting at around $\Delta t=2\times
10^{-2}$, the statical error becomes bigger and bigger in negative values as
$\Delta t$ increases. It attains $-0.05$ at $\Delta t=7\times 10^{-2}$ and
then becomes smaller towards the positive direction as $\Delta t$ increases
further. Finally at $\Delta t=0.1$, the error shifts from negative to positive
and it is expected that the error is not small for all schemes. This
phenomenon that the error is shifting between negative and positive is
observed in all three methods including the GW and the GCC. As a consequence,
in practical applications, it is necessary to adopt the value of time step
size $\Delta t$ within the permissible error during when error becomes larger
as $\Delta t$ becomes larger, instead of the value of $\Delta t$ with the
smallest error.

Figure 5: Kinetic temperature versus time step on linear $y$-axis. SV-AB-4
($\beta=1$) scheme with $\alpha=0.9$ and $\lambda=1.0$. Red and blue curves
represent statical error of kinetic temperatures before and after taking
absolute values. Note that the kinetic temperature is averaged over time after
equilibration.

##  5 Conclusions and outlook

In conclusion, we proposed a novel stochastic Hamiltonian formulation (SHF)
with matrix noise and subject to external forces, which was found applicable
to DPD simulations as the DPD system could be obtained from the corresponding
stochastic Lagrange–d’Alembert principle by introducing proper Hamiltonian
functions and dissipative forces. In particular, we extended the well-known
symplectic SV scheme for conservative Hamiltonian systems to the SHF as
composites of the Euler-A and Euler-B methods. By discretising the dissipative
forces properly, several simple families of SV methods were constructed and
especially the SV-AB methods were focused which were derived as the composite
$\text{(Euler-A)}\circ\text{(Euler-B)}$. By studying the damped Kubo
oscillator, the fluctuating behaviour and damping energy/Hamiltonian
dissipation were realised with order of error approximately $10^{-3}$ between
the numerical and exact Hamiltonians. For DPD simulations, the SV-AB methods
include the conventional GW and GCC methods as special cases. Simulations of a
novel two-parametric explicit schemes were conducted and compared with the GW
and GCC methods. As time step varies, some of the novel schemes were
advantageous over the GW method but unfortunately no global advantage was
realised. It was also observed that for all schemes as the time step increases
the error can shift between positive and negative values, that requires one to
choose a time step in practical applications more carefully.

Beside the SV methods proposed in the current study, thanks to the SHF a
variety of other effective structure-preserving methods may be extended as
well, for instance, symplectic partitioned Runge–Kutta methods and variational
integrators. These are part of our current and future studies including their
applications to the DPD and other relevant stochastic physical systems. From
the theoretical viewpoint, it is worthwhile to study further the geometric and
algebraic structures of the SHF, for instance, conformal symplectic
structures, generating functions, symmetries and Noether’s conserved
quantities.

## 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

This work was partially supported by JSPS KAKENHI Grant Number JP20K14365,
JST-CREST Grant Number JPMJCR1914, and Keio Gijuku Fukuzawa Memorial Fund. We
thank the anonymous referees for their constructive comments.

## References

  * [1] P. J. Hoogerbrugge, J. M. V. A. Koelman, Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics, EPL 19 (3) (1992) 155–160, https://doi.org/10.1209/0295-5075/19/3/001.
  * [2] R. D. Groot, P. B. Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys. 107 (11) (1997) 4423–4435, https://doi.org/10.1063/1.474784.
  * [3] S. Chen, E. Olson, S. Jiang, X. Yong, Nanoparticle assembly modulated by polymer chain conformation in composite materials, Nanoscale 12 (27) (2020) 14560–14572, http://doi.org/10.1039/d0nr01740j.
  * [4] H. Huang, R. Liu, C. A. Ross, A. Alexander-Katz, Self-directed self-assembly of 3D tailored block copolymer nanostructures, ACS Nano 14 (11) (2020) 15182–15192, http://doi.org/10.1021/acsnano.0c05417.
  * [5] N. Arai, Y. Kobayashi, K. Yasuoka, A biointerface effect on the self-assembly of ribonucleic acids: a possible mechanism of RNA polymerisation in the self-replication cycle, Nanoscale 12 (2020) 6691–6698, http://doi.org/10.1039/c9nr09537c.
  * [6] Y.-T. Cheng, H.-K. Tsao, Y.-J. Sheng, Interfacial assembly of nanorods: smectic alignment and multilayer stacking, Nanoscale 13 (33) (2021) 14236–14244, http://doi.org/10.1039/d1nr03784f.
  * [7] S. V. Nikolov, A. Fernandez-Nieves, A. Alexeev, Behavior and mechanics of dense microgel suspensions, Proc. Natl. Acad. Sci. U.S.A. 117 (44) (2020) 27096–27103, http://doi.org/10.1073/pnas.2008076117.
  * [8] L. Pan, F. Wang, Y. Cheng, W. R. Leow, Y. W. Zhang, M. Wang, P. Cai, B. Ji, D. Li, X. Chen, A supertough electro-tendon based on spider silk composites, Nat. Commun. 11 (1) (2020) 1–9, http://doi.org/10.1038/s41467-020-14988-5.
  * [9] Y. Kobayashi, H. Gomyo, N. Arai, Molecular Insight into the Possible Mechanism of Drag Reduction of Surfactant Aqueous Solution in Pipe Flow, Int. J. Mol. Sci. 22 (14) (2021) 7573, http://doi.org/10.3390/ijms22147573.
  * [10] D. P. Papageorgiou, S. Z. Abidi, H. Y. Chang, X. Li, G. J. Kato, G. E. Karniadakis, S. Suresh, M. Dao, Simultaneous polymerization and adhesion under hypoxia in sickle cell disease, Proc. Natl. Acad. Sci. U.S.A. 115 (38) (2018) 9473–9478, http://doi.org/10.1073/pnas.1807405115.
  * [11] P. Chen, H. Yue, X. Zhai, Z. Huang, G. H. Ma, W. Wei, L. T. Yan, Transport of a graphene nanosheet sandwiched inside cell membranes, Sci. Adv. 5 (6) (2019) eaaw3192, http://doi.org/10.1126/sciadv.aaw3192.
  * [12] N. Arai, T. Koishi, T. Ebisuzaki, Nanotube Active Water Pump Driven by Alternating Hydrophobicity, ACS Nano 15 (2) (2021) 2481–2489, http://doi.org/10.1021/acsnano.0c06493.
  * [13] F. Sicard, J. Toro-Mendoza, Armored Droplets as Soft Nanocarriers for Encapsulation and Release under Flow Conditions, ACS Nano 15 (7) (2021) 11406–11416, http://doi.org/10.1021/acsnano.1c00955.
  * [14] J. B. Gibson, K. Chen, S. Chynoweth, The equilibrium of a velocity-Verlet type algorithm for DPD with finite time steps, Int. J. Mod. Phys. C 10 (01) (1999) 241–261, https://doi.org/10.1142/S0129183199000176.
  * [15] T. Shardlow, Splitting for Dissipative Particle Dynamics, SIAM J. Sci. Comput. 24 (2003) 1267–1282, https://doi.org/10.1137/S1064827501392879.
  * [16] X. Shang, Accurate and efficient splitting methods for dissipative particle dynamics, SIAM J. Sci. Comput. 43 (2021) A1929–A1949, https://doi.org/10.1137/20M1336230.
  * [17] B. Leimkuhler, X. Shang, On the numerical treatment of dissipative particle dynamics and related systems, J. Comput. Phys. 280 (2015) 72–95, https://doi.org/10.1016/j.jcp.2014.09.008.
  * [18] B. Leimkuhler, S. Reich, Simulating Hamiltonian Dynamics, Cambridge University Press, Cambridge, 2004.
  * [19] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integrators: Structure-Preserving Algorithms for Ordinary Differential Equations, Springer, Berlin, 2006.
  * [20] G. Milstein, Y. M. Repin, M. V. Tretyakov, Symplectic Integration of Hamiltonian Systems with Additive Noise, SIAM J. Numer. Anal. 39 (2002) 2066–2088, https://doi.org/10.1137/S0036142901387440.
  * [21] L. Wang, J. Hong, R. Scherer, F. Bai, Dynamics and variational integrators of stochastic Hamiltonian systems, Int. J. Numer. Anal. Model. 6 (2009) 586–602, http://global-sci.org/intro/article_detail/ijnam/785.html.
  * [22] D. D. Holm, T. Tyranowski, Stochastic discrete Hamiltonian variational integrators, BIT Numer. Math. 58 (2018) 1009–1048, https://doi.org/10.1007/s10543-018-0720-2.
  * [23] M. Kraus, T. M. Tyranowski, Variational integrators for stochastic dissipative Hamiltonian systems, IMA J. Numer. Anal. 41 (2) (2021) 1318–1367, https://doi.org/10.1093/imanum/draa022.
  * [24] J. E. Marsden, M. West, Discrete mechanics and variational integrators, Acta Numer. 10 (2001) 357–514, https://doi.org/10.1017/S096249290100006X.
  * [25] J.-A. Lázaro-Camí, J. Ortega, Stochastic Hamiltonian dynamical systems, Rep. Math. Phys. 61 (2008) 65–122, https://doi.org/10.1016/S0034-4877(08)80003-1.
  * [26] L. Arnold, Stochastic Differential Equations: Theory and Applications, John Willy & Sons, New York, 1994.
  * [27] L. C. Evans, An Introduction to Stochastic Differential Equations, AMS, 2013.
  * [28] P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer-Verlag, Berlin, 1995.
  * [29] O. D. Street, D. Crisan, Semi-martingale driven variational principles, Proc. R. Soc. A. 477 (2021) 20200957, http://doi.org/10.1098/rspa.2020.0957.
  * [30] T. M. Tyranowski, Stochastic variational principles for the collisional Vlasov–Maxwell and Vlasov–Poisson equations, Proc. R. Soc. A. 477 (2021) 20210167, http://doi.org/10.1098/rspa.2021.0167.