Investigation of dissipative particle dynamics with colored noise

TL;DR

This paper explores a modified dissipative particle dynamics model incorporating time-correlated noise, revealing unique equilibrium properties and potential for applications beyond traditional molecular simulations.

Contribution

It introduces a new stochastic force with algebraic correlation in DPD, demonstrating different system behaviors and equilibrium states compared to classical models.

Findings

Velocity distribution remains Gaussian

Acceleration exhibits bi-modal distribution

System reaches equilibrium with fluctuation-dissipation balance

Abstract

We investigate the behavior of dissipative particle dynamics (DPD) with time-correlated random noise. A new stochastic force for DPD is proposed which consists of a random force whose noise has an algebraic correlation proportional to 1/t and is generated by the so called Kangaroo process. We stress the benefits of a time correlated noise in stochastic systems. We show that the system exhibits significantly different properties from classical DPD, driven by Wiener noise. While the probability distribution function of the velocity is Gaussian, the acceleration develops a bi-modal character. Although the fluctuation dissipation theorem may not strictly hold, we demonstrate that the system reaches equilibrium states with fluctuation-dissipation balance. We believe that our explorative research on the DPD model may stimulate the application of modified DPD to unconventional problems beyond molecular modeling.

Figures

Tables

Table 1. Table 1 : Input parameters of the C-DPD and DPD solvent.

Domain size ($V=L_x\times L_y$)	$60\times 60$
Total number of particles ($N_p$)	${42}^2$
Mass (m)	$1$
Temperature ($k_{B}T$)	$1$
Time step ($dt$)	$0.001$
Density ($\rho$)	$N_p/V$
Simulation length in time step ($N_t$)	${10}^6$
Cut-off radius ($r_c$)	$1$
C-DPD
Dissipative coefficient ($\gamma$)	$0.158$
Stochastic coefficient ($\alpha$)	$1$
Maximal non-correlated time ($\epsilon$)	$0.05$
DPD
Dissipative coefficient ($\gamma$)	$0.5$
Stochastic coefficient ($\sigma$)	$2\gamma k_{B}T$

Table 2. Table 2 : C-DPD realization’s summary

Name	$(\gamma ,\alpha )$	initial condition
$k_{B}T=1$	$(0.158,1)$	Random velocity
Cooling	$(10,1)$	Random velocity
Heating	$(0.05,1)$	Random velocity
Body force 1	$(0.158,1)$	$F_x=\sin \left(2\pi y/L_y\right)$ for $t\le 100$
Body force 2	$(0.158,1)$	$F_x=1.5\sin \left(2\pi y/L_y\right)$ for $t\le 100$

Table 3. Table 3 : Decay of the VACF

C-DPD $\gamma$	$c$
$0.05$	$0.90$
$0.158$	$0.70$
$10$	$0.62$
DPD $\gamma$	$c$
$0.5$	$0.95$

Equations

$$\text{diff}{\mathbf{r}}_{i}t={\mathbf{v}}_i\text{ , }\text{diff}{\mathbf{v}}_{i}t=\frac{1}{m}\sum _{j\ne i}{\mathbf{F}}_{ij}.$$

$${\mathbf{F}}_{ij}={\mathbf{F}}_{ij}^C+{\mathbf{F}}_{ij}^D+{\mathbf{F}}_{ij}^R.$$

$${\mathbf{F}}_{ij}^C=\{\begin{array}{ll}{a}_{ij}m(1-r_{ij}){\hat{\mathbf{r}}}_{ij}\text{ if }& r_{ij}<r_c\\ 0\text{ if }{r}_{ij}\ge r_c& \end{array}.$$

$$F_{ij}^D=-\gamma w^D(r_{ij})({\hat{\mathbf{r}}}_{ij}\cdot {\mathbf{v}}_{ij}){\hat{\mathbf{r}}}_{ij}\text{ ,}$$

$$F_{ij}^R=\sigma w^R(r_{ij}){\theta }_{ij}{\hat{\mathbf{r}}}_{ij}\text{ ,}$$

$$⟨{\theta }_{ij}(t)⟩=0\text{ ,}$$

$$⟨{\theta }_{ij}(t){\theta }_{kl}(t^{\prime })⟩=({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})\delta (t-t^{\prime }).$$

$$w^D=(w^R{)}^2\text{ and }{\sigma }^2=2\gamma k_{B}T$$

$$w^D=(w^R(r){)}^2=\{\begin{array}{ll}(1-\frac{r}{r_c}{)}^2,& \text{if }r\le r_c\\ 0,& \text{otherwise}\end{array}$$

$$C_R(t)=\{\begin{array}{ll}\alpha /\epsilon ,& t\le \epsilon \\ \alpha /t,& t>\epsilon \end{array}\text{ ,}$$

$${\mathbf{F}}_i^R(t)=\left(\begin{array}{c}{m}_{x,i}\\ m_{y,i}\end{array}\right)=\left(\begin{array}{c}\sqrt{\alpha }\cos ({\varphi }_i/\sqrt{\epsilon }))\\ \sqrt{\alpha }\sin ({\varphi }_i/\sqrt{\epsilon })\end{array}\right)\text{ ,}$$

$$\Delta t=\epsilon /{\varphi }_i.$$

$$m\text{diff}{\mathbf{v}}_{i}t=\sum _{j\ne i}{\mathbf{F}}_{ij}^D+{\mathbf{F}}_i^R$$

$$\frac{\partial }{\partial t}\rho (\mathbf{r},t)=$$

$$+\int d{\mathbf{r}}^{\prime }\left[W(\mathbf{r}|{\mathbf{r}}^{\prime },t)\rho ({\mathbf{r}}^{\prime },t)-W({\mathbf{r}}^{\prime }|\mathbf{r},t)\rho (\mathbf{r},t)\right]$$

$$W(\mathbf{r}|{\mathbf{r}}^{\prime },t)=\underset{\Delta t\to 0}{\lim }p(\mathbf{r},t+\Delta t|{\mathbf{r}}^{\prime },t)/\Delta t,$$

$$\frac{\partial }{\partial t}\rho (\mathbf{r},t)=$$

$$+\sum _i\text{partialderivative}{\mathbf{p}}_i\left[\sum _{j\ne i}{\mathbf{e}}_{ij}\gamma w^D({\mathbf{e}}_{ij}\cdot {\mathbf{v}}_{ij})\right]\rho (\mathbf{r},t)$$

$$-\nu (\mathbf{r})\rho (\mathbf{r},t)+Q(\mathbf{r})\int \nu ({\mathbf{r}}^{\prime })\rho ({\mathbf{r}}^{\prime },t)d{\mathbf{r}}^{\prime }$$

$$m\text{diff}{\mathbf{v}}_i(t)t=-{\int }_0^{t}K(t-s)\sum _{j\ne i}{\mathbf{v}}_i(s)ds+{\mathbf{F}}_i^R(t).$$

$$C_R(s)=k_{B}TK(t-s).$$

$$\epsilon =Ndt$$

$$F_x=\sigma \sin (2\pi y/L_y).$$

$$\alpha =5.23{\gamma }^{0.89}.$$

$$\mathcal{D}={\int }_0^{\infty }C(t)dt.$$

$$g(t)=a{e}^{-(bt{)}^c}.$$

$$P(a)=C\exp (-a^2/((1+|a\beta {|}^{\gamma }){\sigma }^2)).$$

Full text

I Introduction

An important aspect in the choice of a numerical model of a physical problem are the length and time-scale defining the quantities of interest. At macroscopic scales continuum-based methods are appropriate, while Molecular Dynamics (MD) models can capture the microscopic details. A large variety of methods has been proposed for the intermediate mesoscopic scale. E.g. for complex fluids or soft matter mesoscopic scales play an important role. Large time and length scale ranges characterize their behavior and prevent the suitability of a single model to resolve simultaneously different scales. Groot and Warren classified dissipative particle dynamics (DPD) as a mesoscopic simulation method Groot and Warren (1997). DPD was introduced as a coarse-grained particle-based stochastic model in a Lagrangian reference frame by Hoogerbrugge and Koelman Hoogerbrugge and Koelman (1992). Effectiveness has been demonstrated for a wide range of problems, such as multiphase phenomena Pagonabarraga and Frenkel (2001), interaction of polymers, surfactants and water Groot (2000) and dynamics of membranes Fedosov et al. (2010).

In this paper we investigate DPD as a numerical model abstracted from a specific application. In Azarnykh et al. (2018) similarities and differences between the Langevin model and DPD have been highlighted. It has been found that the current auto-correlation functions of the $N$-particle subregime of DPD and the Langevin equations coincide. The Langevin principle Langevin (1908) consists in splitting the motion into two parts: the slowly varying motion, and the rapidly changing properties resolved by a random variable. This separation is appropriate if the characteristic time-scale of the system is larger than the time-correlation of the rapidly varying variables that are modeled by noise. In this case the system is Markovian. The stochastic procedure proposed by Langevin disregards the non-Markovian property of the system. Non-Markovian extensions of the Langevin equation have been applied successfully to many physical systems, e.g. Brustein et al. (1991); Ceriotti et al. (2009).

Recently, Li et al. Li et al. (2015) proposed to incorporate explicit memory effects into DPD via the Mori-Zwanzig formalism. It was found that when there is a lack of time-scale separation, the non-Markovian DPD shows little improvement on the velocity autocorrelation function compared with the standard DPD model. Motivated by this observation, we aim at manipulating memory effects of DPD implicitly. Unlike Li et al. (2015), where the random noise generator is derived according to the second fluctuation dissipation theorem (FDT) from the memory-kernel, we select a specific colored-noise generator and subsequently analyze the friction term together with the FDT. No explicit memory kernel is used here. The Kangaroo process is used as colored-noise generator Peter Hänggi and Peter Jung (1995) which has a slowly decaying noise-correlation proportional to $1/t$. In Srokowski and Płoszajczak (1998) a non-Markovian Langevin equation with a colored noise generated by the Kangaroo process has been proposed. In this case the velocity autocorrelation function (VACF) can be analytically derived, and in the case of the $1/t$ autocorrelation noise, the VACF has an algebraic tail. Also, behavior which is typical for intermittent structures of Lévy flights has been observed in Srokowski and Płoszajczak (1998). We note that physical systems exhibit a slowly decaying correlation are e.g. atomic diffusion through a periodic lattice Igarashi and Munakata (1988) and ligand migration in biomolecules Hanggi (1986).

The paper is structured as follows: In section II, the model is presented in detail together with the theoretical background. Section III examines the fluctuation dissipation theorem (FDT) in our context. Numerical results are analized in section IV. The latter also provides comparisons between the standard DPD and the DPD with colored-noise (C-DPD). Finally, we conclude with a discussion, summary of our results and a brief outlook of future research and ongoing work.

II The DPD Model

For standard DPD, the motion of each DPD particle is described by

[TABLE]

The total force acting on particle $i$ is composed of three pairwise interactions

[TABLE]

All forces vanish beyond a cutoff radius $r_{c}$ which we choose to be unity. The conservative force acts along the connecting line between two particles and is soft repulsive

[TABLE]

In the following we do not consider the conservative force $F^{C}$. The dissipation and random forces are given by

[TABLE]

where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$, $\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$, $r_{ij}=|\mathbf{r}_{ij}|$ and $\hat{\mathbf{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$. The functions $w^{D}$ and $w^{R}$ are weighting functions. Furthermore, $\theta_{ij}$ is a Gaussian random variable with

[TABLE]

In order to satisfy the fluctuation dissipation theorem, Español and Warren Español and Warren (1995) showed that

[TABLE]

must hold. In Eq. (6) $k_{B}$ is the Boltzmann constant and $T$ the temperature. We choose the standard kernel function

[TABLE]

where $r_{c}$ is the radius cut-off. The velocity Verlet algorithm is used for time integration. The equation of motion given by (1) and (2) is Markovian since the additive noise $\theta$ is not correlated in time.

Let us now consider a noise with the following auto-correlation function Srokowski and Płoszajczak (1998)

[TABLE]

for some $\varepsilon>0$. The variable $\alpha$ is the amplitude of the noise. The stochastic process called ”the Kangaroo process” (KP) provides such a behavior Brissaud and Frisch (1974). Its name originates from its resemblance to Lévy flights. The process is defined as a stepwise random function which is constant during a time interval $\Delta t$ which is a random number. Although the Kangaroo process can be formulated for higher dimensions, here we focus on two dimensional systems. The resulting stochastic force for particle $i$ is given by

[TABLE]

where $\phi_{i}$ is a random number uniformly distributed in the interval $(0,1)$. The force $F^{R}_{i}$ stays constant during an interval $\Delta t$ which is defined as

[TABLE]

If $t>t_{n}+\Delta t$ a new $\phi_{i}$ is generated as well as a new time interval within which the process stays unchanged, and so on. In Fig. (1) the force distribution is shown.

The random force (8) is inserted into the equation of motion

[TABLE]

and gives the C-DPD version of DPD. Assuming the system to be in equilibrium with a heat bath of temperature $T$, it has to satisfy the fluctuation dissipation theorem for this given $T$ Berne and Harp (2007). The fluctuations must be linked to the dissipation in such a way that no loss or gain of energy occurs in the system. The temperature has to fluctuate around a constant. It is for many reasons not straightforward to fulfill this condition. Firstly, the forces on the particle $i$ are of different kind, one is a pairwise interaction while the other is a background random fluctuation. Secondly, the stochastic force does not vanish for $r>r_{c}$. The dissipative force is defined through a kernel and has finite support.

Since the stochastic force is time correlated a discussion on the necessity of a memory kernel as in Li et al. (2015) is in order. One standard procedure would be to derive the Fokker-Planck equation and solve it for its steady state $\rho^{eq}$ Marsh and College ; Español and Warren (1995). The Fokker-Planck equation is a special case of the forward Chapman-Kolmogorov equation (CKE)

[TABLE]

for which the integral term vanishes. For processes with jumps we have to determine $A$, $B$ and

[TABLE]

where $p$ is a transition probability defined by the nature of the chosen random process. According to Kamińska and Srokowski (2003), for KP functions $A$ and $B$ vanish. The CKE for DPD with colored noise generated by the KP can be superimposed by the functions $A$ and $B$ for Eq. (10) without stochastic force and by the jump integral for the KP random term

[TABLE]

where $\nu$ is the jump frequency and $Q$ the size of a jump. In order to address the question whether a FDT can be satisfied we assume that the dissipative term in Eq. (10) is convolved with a memory kernel $K$. Eq. (10) takes the form of the generalized Langevin equation

[TABLE]

From the latter equation one can see that the FDT is satisfied if

[TABLE]

As approximation, the memory kernel $K$ can be neglected beyond a time cut off Li et al. (2015), i.e. the history of the system is taken into account only for $t\leq Ndt$, where $dt$ is the time step of the integration-scheme. So, if $t\leq Ndt$ the shape of the memory kernel follows from the noise auto-correlation function, otherwise it vanishes. However for the Kangaroo process, upon choosing

[TABLE]

the random force produces a constant memory kernel, see Li et al. (2015), Eq. (16), for $0<t\leq\varepsilon$, and time can be collapsed into a recalibrated DPD dissipation term.

III The fluctuation dissipation theorem

Unless specified differently, simulations presented in this paper have been performed with parameters presented in Table 1. Each particle is initialized with a random velocity.

The fluctuation dissipation theorem (FDT) is a crucial element in determining the prevailing thermodynamic regime. When it applies, the system temperature coincides with the actual physical temperature. For this reason it is necessary to check whether a FDT is satisfied or not Español and Warren (1995). In cases where an analytic derivation is difficult, one may verify whether the conditions for the FDT are satisfied by simulations. First, the system should reach a stationary state, which is a necessary condition. Furthermore, the stationary state has to be unique. The system has to fluctuate around a constant temperature, independently of the chosen initial condition.

For $\alpha=1$, different friction coefficients $\gamma$ are tested. The results are compared with a DPD simulation with $\gamma_{DPD}=0.5$. The first observation is that for a specific choice of $\gamma$, the temperature of C-DPD fluctuates around the same equilibrium value as that of DPD, see Fig. (2). When the stochastic force dominates the dissipation of the system, the system heats up and reaches a constant temperature after some transient. When the dissipation force overwhelms the effects of the stochastic term, the system cools down. In Fig. (2), for the yellow set of data points, the temperature fluctuates around $T=0.34$ but remains stationary after $t\approx 40$, which implies that the system has cooled to its equilibrium temperature. For the green data set the temperature grows until it also reaches an equilibrium value. These observations are further discussed in sections IV and V.

In Fig. (3) the probability density function (pdf) of the velocity in the $y$-direction at different times is represented. It appears that apart from small deviations due to the finite sample size (finite number of particles) the pdf is stationary since it does not change with time. Furthermore, in Fig. (4) the variance of the squared velocity is plotted in a lin-log scale which shows that a stationary state is reached. Also the variance of the latter quantity stays in the same range as the one obtained from standard DPD.

In the following we address the issue whether the C-DPD system reaches a unique stationary state, independently of the initial state. Three different simulations are performed, all with the same coefficients for the dissipative and random forces. For each one different initial data are chosen. In Fig. (5) the temperature evolution for these three data sets are plotted. The data ”Body force 1” corresponds to a pre-run with a body force $F_{x}=\sin(2\pi y/L_{y})$, data ”Body force 2” corresponds to body force $F_{x}=3/2\sin(2\pi y/L_{y})$ each applied until $t=1$ and turned off afterwards. See Table 2. The first case, without body force, corresponds to a simulation for which the initial velocity has been chosen randomly but no body force is applied to the flow. For the other two cases, with body force 1 and 2, for the first 1000 iterations an external body force $F_{x}$ is applied and then turned off. We choose $F_{x}$ to be sinusoidal

[TABLE]

We see that independently of the initial state, the system relaxes towards the same stationary state. The pdf of the velocity and the acceleration of these three examples do not differ from each other.

We here performed many realizations with different stochastic coefficients $\alpha=1,2,3,6$. For each value of $\alpha$ one can find a friction coefficient for which the system reaches a stationary state. By iterating this process and a fitting procedure, we calibrate a $\alpha(\gamma)$ function, Eq. (16). Empirically we have established that a FDT may be satisfied for the relation

[TABLE]

IV Discussion of results

Results presented in this section have been produced using mostly the empirical fluctuation-dissipation relation (16). In Fig. (7(b)), a Gaussian (orange curve) has been used as fit to indicate that the pdf of velocity is Gaussian. While for DPD the pdf of the acceleration also follows a Gaussian, for C-DPD it results in a bimodal Gaussian. The realization used to generate Fig. (7) correspond to $k_{B}T=1$ in Table 2. The system has two different most likely states.

This kind of behavior may be explained by the coexistence of two different populations at stationary state. They are categorized by their acceleration. A particle may reside in one population for a finite time (which depends on the parameter $\varepsilon$) before it jumps to the other. One expects trajectories with short segment of free flight before a sudden change of direction, see Fig. (8), where trajectories for the first $20$ time steps are shown. Note that the notion of free flight is used in a descriptive sense. During such free flights, the particle undergoes constant acceleration. While for the orange trajectory no change in direction yet has happened, the red as well as the green trajectories clearly show such events. Also in comparison to the DPD, for the three C-DPD simulations free flights are observed.

The velocity auto-correlation function (VACF) exhibits a power law such as Eq. (18) for the DPD as well as for the C-DPD, Fig. (10). The VACF is connected to the mean free path Zumofen and Klafter (1993). The C-DPD and the DPD VACF are very similar. Two regimes can be distinguished, first an exponential decay and at later times an algebraic dependence. The VACF decays rapidly towards their stationary value $C_{\infty}$. For both models $C_{\infty}$ is almost zero. No significant negative tail can be observed, which indicates essentially classical diffusion.

The diffusion coefficient can be computed by means of the VACF

[TABLE]

The results in Fig. (9) show that the relation between the dissipation coefficient and the friction coefficient $\gamma$ for both versions of DPD are very similar.

Referring to the issue of the cooling and heating behavior, Fig. (2), we observe that for realizations where the dissipation overwhelms the fluctuation the equilibrium temperature of the system is smaller than the input temperature. However, it assumes a stationary value after a transient. When the dissipation is almost negligible, the temperature grows, but the system again reaches a stationary temperature value.

We investigate both cases where the parameter $\gamma$ is either too big or too small to satisfy relation (16). In Fig. (11) the VACF for these cases are plotted for constant $\alpha$. The black lines correspond to the fit

[TABLE]

The calibrated parameters $c$ for the corresponding realization are presented in table 3.

For the heating system ($\gamma=0.05$) the VACF reveals a time interval where the auto-correlation stays negative for finite time which cannot be observed for the cooling system ($\gamma=10$). However, in Fig. (12) we see that the acceleration pdf of the latter differs strongly from the behavior of the heating C-DPD and the standard DPD system. We observe non-Gaussian tails of the pdf for which the correlated noise is responsible. We point out that experimental measurements in La Porta et al. (2001) reported that the Lagrangian acceleration pdf in fully developed turbulence follows

[TABLE]

The fit of the realization for the cooling system with the latter function is in surprisingly good agreement with the experiment. The corresponding parameters $C$, $\alpha$, $\gamma$ and $\sigma$ are given in Fig. (12).

V Summary and discussions

The dissipative particle dynamics (DPD) model with correlated stochastic force reveals phenomena beyond of what has been observed so far for standard DPD. In the present paper we have formulated a colored DPD system for which a fluctuation dissipation relation is found empirically without explicit memory kernel in the dissipative force Li et al. (2015). Results show that a proper calibration of the coefficients of the different forces results in a system fluctuating around a constant temperature. The system reaches a unique stationary state, e.g. also shown by the fact that the second moment of the velocity distribution remains bounded with time.

Transport properties determined by the velocity autocorrelation function alone do not differ from standard DPD. Typical trajectories for Levy flights in the configuration space have been observed with the numerical realizations.

The results of the present paper encourage further investigations of C-DPD with a non-zero conservative force. Also we will consider to investigate the formulation and the effect of a memory kernel which is consistent with the random force.

VI Acknowledgments

We would like to acknowledge support by the Deutsche Forschungsgemeinschaft within the Priority Programme Turbulent Superstructures under Grant no AD 186/30.