## Feature #1885

### DPD Thermostat

**Description**

Dear Gromacs Developers,

I propose to develop a new patch that will include the DPD thermostat. Our MD Group from Groningen is willing to take care of the coding and further maintenance. The theory is developed together with Herman Berendsen. I selected bellow the algorithms that we want to implement and explain as compact as possible (selection made from two articles that we published). Please have a look bellow. It will take no more than 5-10 minutes.

I already started a patch in Gromacs development. Let me know of any comment that you have.

Best regards and Happy New Year

Nicu

1. Dissipative Particle Dynamics
***********************************

The essence of Didssipative Particle Dynamics (DPD)-like friction and noise is that it is applied to velocity differences (or relative velocities) between pairs of particles. If done properly, the pairwise application ensures conservation of total linear momentum. In Gromacs the friction and noise is applied isotropically (ISO) to the velocity difference vector, irrespective of its direction. After determining the change in the relative velocity due to the impulsive friction and noise, this change is distributed over the two particles such that the c.o.m. is conserved. ISO is an option of the mdrun program.

The particle pair selection can be done in several ways. In Gromacs implementation we will chose a pair particle at random from the particles of a simulated molecular system. The DPD algorithm looks like follows.

For each selected pair (i, j) with velocities v_i, v_j, do the following:

- Choose the velocity reduction factor f.

`- Determine the velocity noise factor g , defined as`

g = \sqrt{f(2-f) k_B*T_ref/ \mu,

where \mu = m_i * m_j/ (m_i+ m_j) is the reduced mass of the two particles.

- Construct the relative velocity vector 'v' of the selected pair:

v = v_i - v_j.

- Choose 3 random numbers xi= (xi_1,xi_2,xi_3) from a standard normal distribution (mean=1, sd=1).

- Construct the vector

v = - f * v + g * \xi

Distribute the relative velocity change over the two particles:

v_i = v_i + \mu/m_i * v

v_j = v_j - \mu/m_j * v

In this way the velocities of the particles are updated while the total momentum m_i * v_i+ m_j * v_j is conserved.

Note that the particle velocities are updated after each impulsive event, not at the end of each step. This is necessary, as a single particle may be involved in more than one pair events; first adding the velocity changes and updating the velocities with the sum of the velocity changes gives erroneous results.

1.2 Dissipative Particle Dynamics With Constraints
************************************************************

The algorithm from the previous section is modified for the constraints case in the following way.

1. Compute unconstrained forces $\ve{F}^u = -\ve{\nabla} V(\ve{x}(t))$ and compute virial using $\ve{x}(t)$

2. Advance velocities $\forall{i}: v_i = v_i(t - \frac12 \Delta t) + (F_i^u/m_i) \Delta t$

3. Advance positions $\forall{i}: x_i = x_i(t) + v_i \Delta t$

4. Apply constraints to new positions: $(\ve{x}^c, \ve{v}^c) = $\texttt{constr}$(\ve{x};\ve{x}(t), \ve{v})$ and add constraint contribution to virial

5. Obtain a Maxwellian velocity distribution $v^m$, constrain it and obtain $ v^{m,c}$

6. Add to velocity the Maxwelian contribution $\forall{i}: v_i^c = v_i^c + (v_i^m - v_i^{m,c})$

7. Advance the positions with the Maxwelian contribution $\forall{i}: x_i = x_i + (v_i^m - v_i^{m,c}) \Delta t$

8. Apply DPD friction and noise isotropically (as described in the previous section) and obtain $v_i'$

9. Correct positions for velocity change: $\forall{i}: x_i' = x_i^c \frac12 (v_i' - v_i^c) \Delta t $

10. Apply constraints to corrected positions: $(\ve{x}(t\Delta t), \ve{v}_i(t+\frac12 \Delta t) = $\texttt{constr}$(\ve{x}';\ve{x}(t), \ve{v}'_i)$; do \textbf{not} add constraint contribution to virial

Note that the constraint algorithm is applied twice: the first time in step 4 to correct the effect of the unconstrained forces (this application also provides the conservative constraint force and computes the constraint contribution to the virial) and the second time in step 10 to annihilate the effects of friction and noise in the constraint directions. The latter concerns non-conservative forces acting on the velocities, which do not contribute to the virial. Steps 5 to 7 are needed for compensating the loss of degrees of freedom through the application of constraints (first time in step 4). This compensation is needed for the correct functioning of DPD-ISO algorithm in step 8 (for a larger discussion on this topic we refer to~\cite{Goga2014}).

Bibliography
***************

[Goga2012] N. Goga, A. J. Rzepiela, A. H. de Vries, S. J. Marrink and H. J. C. Berendsen, Efficient Algorithms for {Langevin} and {DPD} Dynamics, Journal of Chemical Theory and Computation, 2012, vol. 8, pp 3637--3649

[Goga2014] E. A. J. F. Peters, N. Goga, H. J. C. Berendsen, Stochastic Dynamics with Correct Sampling for Constrained Systems, Journal of Chemical Theory and Computation, 2014, vol 10,pp 4208--4220

**Related issues**

### History

#### #1 Updated by Nicolae goga about 5 years ago

Friends,

Happy New Year to everyone!

I see no reply yet to this open issue. I modified the patch to have no conflicts with accepted patches (suggestion of Teemu Murtola - thanks). Please have just a quick look and comment please. Code I will improve if the basis is Ok. Thank you in advance.

Best regards

Nicu

#### #2 Updated by Berk Hess almost 5 years ago

Looking at the code, for each particle you seem to randomly take another particle in the system (without DD) or in the DD domain and apply distance independent DPD friction and noise. This can't give proper DPD, since that the friction and noise should be spatially localized. Or am I missing something?

#### #3 Updated by Nicolae goga almost 5 years ago

Dear Berk,

Thanks. This is a variant of DPD, ISO type (from the article) that has whatever is needed:

- maintain temperature at the desired one

- conserves the total momentum.

Is perfectly possible to do it also like this.

Best regards

Nicu

Berk Hess wrote:

Looking at the code, for each particle you seem to randomly take another particle in the system (without DD) or in the DD domain and apply distance independent DPD friction and noise. This can't give proper DPD, since that the friction and noise should be spatially localized. Or am I missing something?

#### #4 Updated by Berk Hess almost 5 years ago

But is there a use for such a DPD thermostat? Conserving the overall momentum is useless, since the overall momentum is an independent degree of freedom.

#### #5 Updated by Nicolae goga almost 5 years ago

Dear Berk,

Thanks. Sorry for a late reply. It is not about the total momentum of the system. But about the total momentum of any two particles 'i' and 'j' that participate in thermostating ( mi*vi + mj*vj) that is conserved.

Best regards

Nicu

#### #6 Updated by Mark Abraham about 2 years ago

**Related to***Feature #1627: DPD integrator*added