.. _DPD_Intro:

Intro to Dissipative Particle Dynamics
======================================

Summary
^^^^^^^

Dissipative Particle Dynamics (DPD) is a mesoscale modelling method with many similarities to molecular dynamics, i.e. calculating and integrating forces over a discrete time step to evolve the positions and velocities of particles. DPD features a special pairwise thermostat to control the system temperature while maintaining momentum conservation to ensure correct hydrodynamics. A broad definition of a DPD particle (or 'bead') at length scales larger than those for atoms allows the use of soft interaction potentials, allowing longer times to be achieved with fewer time steps than ordinarily available for molecular dynamics.

Background
^^^^^^^^^^

DPD was originally conceived as an off-lattice form of gas automata [Hoogerbrugge1992]_ to describe complex fluid flows at larger length scales than those available with molecular dynamics (MD). This algorithm introduces dissipative and random forces between pairs of particles:

.. math :: \mathbf{F}_{ij}^{D} = - \gamma w^{D}\left( r_{ij} \right) \left( {\widehat{\mathbf{r}}}_{ij} \cdot \mathbf{v}_{ij} \right){\widehat{\mathbf{r}}}_{ij}
  :label: dissipative

.. math :: \mathbf{F}_{ij}^{R} = \sigma w^{R}\left( r_{ij} \right)\zeta_{ij}\Delta t^{- \frac{1}{2}\ }\ {\widehat{\mathbf{r}}}_{ij},
  :label: random

which control the particles' kinetic energy while ensuring the system momentum is conserved. In the above equations, :math:`\gamma` and :math:`\sigma` are the dissipative and random force parameters, :math:`w^{D}` and :math:`w^{R}` are functions of distance for dissipative and random forces respectively, :math:`\mathbf{r}_{ij} = \mathbf{r}_{j} - \mathbf{r}_{i}` is the vector between particles :math:`i` and :math:`j`, :math:`{\widehat{\mathbf{r}}}_{ij} = \frac{\mathbf{r}_{ij}}{r_{ij}}` is the unit vector between the same particles, :math:`\mathbf{v}_{ij} = \mathbf{v}_{j} - \mathbf{v}_{i}` is the relative velocity between the two particles, :math:`\zeta_{ij}` is a Gaussian random number with zero mean value and unity variance, and :math:`\Delta t` is the simulation timestep.

The connection between the dissipative and random forces to allow them to act as a Galilean-invariant thermostat was made later by Español and Warren [Español1995]_ who made use of the Fokker-Planck equation for fluctuation-dissipation. To ensure that any equilibrium structure is not affected by the dissipative and random forces, the following relationships need to apply:

.. math :: \sigma^{2} = 2\gamma k_B T
  :label: sigmagamma

.. math :: w^{D}\left( r_{ij} \right) = \left\lbrack w^{R}\left( r_{ij} \right) \right\rbrack^{2}
  :label: screen

where :math:`k_B` is the Boltzmann constant and :math:`T` is the required system temperature. With these conditions, the dissipative and random forces make up the DPD thermostat. (In essence, this can be considered as a pairwise form of the Langevin thermostat.) Any flow field applied to a DPD system will be treated correctly and the correct hydrodynamics may be observed.

Beyond controlling the system temperature, no restriction on how particles in a DPD simulation interact is generally applied - indeed, it is possible to carry out atomistic or coarse-grained MD simulations with the DPD thermostat. For the purposes of modelling at the mesoscale, the particles should ideally be comparatively large and soft (compared with atomistic MD) to ensure larger simulation timesteps can be used. One very common form of conservative interaction between DPD particles (or 'beads') is that proposed by Groot and Warren [Groot1997]_, using a pairwise force that is linear with distance:

.. math :: \mathbf{F}_{ij}^{C} = \left\{ \begin{matrix} A_{ij}\left( 1 - \frac{r_{ij}}{r_{c}} \right){\widehat{\mathbf{r}}}_{ij} & (r_{ij} < r_{c}) \\ 0 & (r_{ij} \geq r_{c}) \\ \end{matrix} \right.\ ,
  :label: grootwarren

where :math:`r_c` is an interaction cutoff distance and :math:`A_{ij}` is a conservative force parameter. This results in a quadratic potential :math:`U_{ij} = \frac{1}{2} A_{ij} r_c \left(1 - \frac{r_{ij}}{r_c} \right)^2` and also, for systems with a single particle species, a quadratic equation of state:

.. math :: p = \rho k_B T + \alpha A_{ij} \rho^2 ~~~~~ (\rho > 2)

where :math:`\alpha \approx 0.101 r_c^4` and :math:`\rho` is the overall particle density with units of :math:`r_c^{-3}`. The above equation of state applies when the particle density is greater than 2.

In spite of the simplicity of this interaction, it has been shown to provide enough detail for many complex systems, particularly when used in conjunction with bonded interactions (harmonic springs etc.) between beads to represent amphiphilic molecules with hydrophobic (water-hating) and hydrophilic (water-loving) parts. The softness of the potential makes it possible to attain equilibrium structures, e.g. mesophases, micelles, vesicles, in a relatively short time thanks to the use of larger timesteps.

Further extensions to the basic Groot-Warren 'DPD' interaction have included density-dependent (many-body DPD) potentials, electrostatic interactions with short-range charge smearing etc.

.. [Hoogerbrugge1992] PJ Hoogerbrugge and JMVA Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, *EPL*, **19**, p. 155-160, 1992, doi: 10.1209/0295-5075/19/3/001.

.. [Español1995] P Español and P Warren, Statistical mechanics of dissipative particle dynamics, *EPL*, **30**, p. 191-196, 1995, doi: 10.1209/0295-5075/30/4/001.

.. [Groot1997] RD Groot and PB Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, *Journal of Chemical Physics*, **107**, p. 4423–4435, 1997, doi: 10.1063/1.474784.