A stochastic Hamiltonian formulation applied to dissipative particle dynamics

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.

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].

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)

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$.

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.

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.