# softmatter:Dissipative Particle Dynamics Simulation (DPD)

Dissipative particle dynamics (DPD) is a mesoscale simulation technique invented by Hoogerbrugge and Koelman to study suspension flows. The method was therefore designed to satisfies the Navier-Stokes equations [1, 2]. The first advantage of DPD is the inclusion of hydrodynamic interactions while avoiding computationally expensive calculations. The second advantage is that DPD is able to use time steps up to an order of magnitude larger than those typically used in Brownian Dynamics Simulation (BD) or Molecular Dynamics Simulation (MD) to integrate the discretized equations of motion, e.g. ΔtDPD = 0.06 vs. ΔtMD = 0.005 vs. ΔtBD = 0.01. The ability to integrate the equations of motion using such large time steps results from particles interacting via a soft interaction potential, as described below. DPD rigorously samples the canonical ensemble; the system is run as NVT where the Number of particles, Volume of the box, and Temperature are all fixed (for a more detailed explanation see Thermodynamic Ensembles).

## Method

As in BD and MD, the time evolution of a DPD particle is governed by Newton’s equations of motion such that the positions, velocities and forces are related by the following equations.

$\frac{d\mathbf{r}_{i}}{dt} = \mathbf{v}_{i}$

$\frac{d\mathbf{v}_{i}}{d'} = \mathbf{f}_{i}$

ri, vi, and fi are the position and velocity of particle i and the force acting on particle i respectively. The forces acting on pairs of particles can be separated into three parts:

$\mathbf{f}_{i} = \sum_{j \neq i}(\mathbf{F}_{ij}^C + \mathbf{F}_{ij}^R + \mathbf{F}_{ij}^D)$


a conservative force, random force, and dissipative force, respectively. The conservative force acts between centers of particles and is usually assumed to be a linear function of the distance between the particles.

An example of a linear conservative force is:

$\mathbf{F}_{ij}^C = \begin{cases} a_{ij}(1-r_{ij})\mathbf{\hat{r}}_{ij} & r_{ij} < 1 \\ 0 & r_{ij} \geq 1 \end{cases}$


aij is the repulsion parameter between species i and j, rij is the distance between the centers of particles i and j, and $\mathbf{\hat{r}}_{ij}$ is a unit vector pointing along the line of centers from particle i toward particle j. Note that the equations of motion are stable for particle separations of zero. This is possible because of the linear conservative force. Particles can pass through each other and therefore DPD does not describe the excluded volume of real particles. The random and dissipative forces are described by:

$\mathbf{F_{ij}^R} = \sigma w^R(rij)\theta_{ij}\mathbf{\hat{r}}_{ij}$

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


wR and wD are weighting functions dependent on the distance between particles i and j and $\mathbf{v}_{ij}$ is the relative velocity between particles i and j. θij is a randomly fluctuating variable that satisfies Guassian statistics and γ and σ are the amplitudes of the dissipative and random forces respectively. The amplitudes of these forces are related to the thermal energy kBT such that σ2 = 2γKBT.

## Weighting Functions

The random and dissipative forces are not as straightforward as they may seem. To address these subtleties we consider the equilibrium distribution of the particles. Espanol and Warren derived the following Fokker-Planck equation for DPD [3].

$\frac{\partial \rho(\mathbf{r}_{i},\mathbf{p}_{i},t)}{\partial t} = \mathcal{L}^C \rho(\mathbf{r}_{i},\mathbf{p}_{i},t)+\mathcal{L}^D \rho(\mathbf{r}_{i},\mathbf{p}_{i},t)$

where $\rho(\mathbf{r}_{i},\mathbf{p}_{i},t)$ is the probability of finding the system in a state where the particles have positions $\mathbf{r}_{i}$ and momentum $\mathbf{p}_{i}$ at a given time t. $\mathcal{L}^C$ is the Liouville operator of the conservative force and $\mathcal{L}^D$ is the Liouville operator for the dissipative and random forces. For a system interacting with just conservative forces, the equilibrium probability distribution $\rho(\mathbf{r}_{i},\mathbf{p}_{i})$ is not a function of time and therefore the equation

$\frac{\partial \rho(\mathbf{r}_{i},\mathbf{p}_{i})}{\partial t} = \mathcal{L}^C \rho(\mathbf{r}_{i},\mathbf{p}_{i},t =0)$

is satisfied.

For the canonical ensemble

$\rho^{eq}(\mathbf{r}_{i},\mathbf{p}_{i}) = exp\left( -\frac{1}{k_{B}T}\sum_{i} \left[ \frac{\mathbf{p}_{i}^2}{2m}-U(\mathbf{r}_{i}) \right]\right)$

satisfies the equilibrium distribution. Since, the distribution should not change when the dissipative and random forces are included $\mathcal{L}^D \rho(\mathbf{r}_{i},\mathbf{p}_{i},t)$ must be equal to zero at equilibrium. Therefore, the random and dissipative forces must be related in a manner that satisfies this equilibrium condition. This condition can be satisfied by setting the weighting functions such that wD(rij) = (wR(rij))2 [3].

To simplify the equations of motion the weighting functions usually take on the following form:

$w^D(r_{ij})=[w^R(r_{ij})]^2= \begin{cases} (1-r_{ij})^2 & r_{ij} < 1 \\ 0 & r_{ij} \geq 1 \end{cases}$.


However, any form of the weighting functions can be used providing they satisfy the above constraint.

## Integrating the equations of motion

The equations of motion are integrated with a modified Velocity Verlet Algorithm described by the following equations [4]:

$\mathbf{r}_{i}(t+\Delta t) = \mathbf{r}_{i} +\mathbf{v}_{i}\Delta{t}+\frac{\mathbf{f}_{i}(t)(\Delta t)^2}{2},$

$\tilde{\mathbf{v}}_{i}(t+\Delta t) = \mathbf{v}_{i}(t) + \lambda \mathbf{f}_{i}(t)\Delta t,$

$\mathbf{f}_{i}(t+\Delta t)=\mathbf{f}_{i}(\mathbf{r}_{i}(t+\Delta t), \tilde{\mathbf{v}}_{i}(t+ \Delta t)),$

$\mathbf{v}_{i}(t+ \Delta t)= \mathbf{v}_{i}(t)+\frac{\mathbf{f}_{i}(t)}{2}+\mathbf{f}_{i}(t+ \Delta t))\Delta t$


where $\tilde{\mathbf{v}}_{i}$ represents a prediction of the new velocity which is corrected once the force has been calculated. More indepth details of this method can be found in reference [3] . Guidelines for choosing the time step and corresponding amplitudes of the random and dissipative forces are provided in reference [5].

## Issues

Several issues arising from the soft conservative force should be briefly discussed. First, as mentioned earlier, DPD does not properly capture the excluded volume interactions of particles and therefore the geometric affects of the particle are not captured in the local packing. Second, because the particles can pass through each other the Schmidt number has a value close to unity. The Schmidt number is a measure of the momentum diffusion compared to the mass diffusion and values of unity correspond to gas-like phases where the particles can move freely. For dense phases such as liquids, the Schmidt number is usually several orders of magnitude higher due to the relatively slow mass diffusion resulting from the restricted mobility of the particles. Therefore, caution should be taken whenever dynamics are considered to be important.

## References

• [1] J. M. V. A. Koelman and P. J. Hoogerbrugge, "Dynamic Simulations of Hard-Sphere Suspensions Under Steady Shear," Europhysics Letters, 21, 363-368, (1993)
• [2] P. J. Hoogerbrugge and J. Koelman, "Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics," Europhysics Letters, 19, 155-160, (1992)
• [3] P. Espanol and P. Warren, "Statistical-Mechanics of Dissipative Particle Dynamics," Europhysics Letters, 30, 191-196, (1995)
• [4] M. P. Allen and D. J. Tildesley, " Computer Simulation of Liquids," Clarendon Press: Oxford (1987)
• [5] R. D. Groot and P. B. Warren, "Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation," Journal of Chemical Physics, 107, 4423-4435, (1997)