™

ri = vi (11.1)

™

vi = F i = F ij , (11.2)

j=i

where

F ij = F C + F D + F R . (11.3)

ij ij ij

The conservative force on particle i due to j is repulsive with a range 1 given

by:

r ij

F C = aij (1 ’ rij ) rij < 1, (11.4)

ij

rij

= 0 rij ≥ 1, (11.5)

with r ij = r i ’ r j . This corresponds to a quadratic repulsive potential with

one parameter aij . Note that the distance is scaled such that the maximum

interaction range equals 1. The dissipative force is given by

r ij

F ij = ’γwD (rij )(v ij · r ij ) 2. (11.6)

rij

It acts in the direction of r ij and is proportional to the component of the

velocity di¬erence in the interparticle direction, being repulsive when par-

ticles move towards each other and attractive when they move away. Thus

308 Dissipative particle dynamics

it damps the relative motion of the two particles. The parameter γ mea-

sures the strength of the damping; wD (rij ) is a weight function vanishing

for rij > 1.

The random force also acts in the interparticle direction:

r ij

F R = σwR (rij ) θij , (11.7)

ij

rij

where σ is the strength, wR (rij ) a weight function vanishing for rij > 1, and

θij a random function with average zero and with no memory: θ(0)θ(t) =

δ(t), uncorrelated with the random function on any other particle pair. The

distribution function of θ can be chosen to be normal, but that is not

a requirement. Espa˜ol and Warren (1995) showed that the ¬‚uctuation“

n

dissipation theorem, ensuring that the energy changes from dissipation and

random force cancel, requires that

wD (rij ) = [wR (rij )]2 , (11.8)

and that the noise intensity must be related to the friction coe¬cient γ and

the temperature T , just as is the case for Langevin dynamics:

σ 2 = 2γkB T. (11.9)

The form of one of the weight function is arbitrary; Groot and Warren (1997)

chose for wR the same functional form as for the conservative force:

r ij

wR (rij )C = (1 ’ rij ) rij < 1, (11.10)

rij

= 0 rij ≥ 1. (11.11)

11.2 Prescribing ¬‚uid parameters

The unit of length has been set by the choice of the maximum interaction

range. The number density n (per cubic length unit) can be chosen, and, to-

gether with a choice of a, the strength of the conservative force, the pressure

is ¬xed. The pressure is found from the usual virial equation:

N 1

r ij · F C .

p= kB T + (11.12)

ij

V 3V

i,j>i

In the virial part also the total force may be used, but in equilibrium the

contributions of friction and random force cancel. One may wish to match

the isothermal compressibility with a desired value for a given liquid. This

11.3 Numerical solutions 309

is best expressed as the dimensionless value

1 ‚p

κ’1 = , (11.13)

kB T ‚n

which has a value of 16 for water at 300 K and 30 for 1-propanol at 300 K.

From a series of simulations Groot and Warren (1997) found that

an

, ± = 0.101 ± 0.001.

p = nkB T + ± (11.14)

kB T

This determines an/kB T , which is equal to 75 for water.

11.3 Numerical solutions

Since the force F(t) depends on the equal-time velocity v(t), the normal

Verlet-type algorithms cannot be used, because they have the equal-time

velocities v(t) available after a step that has already used the force F(t).

This applies also to the velocity-Verlet version. If earlier velocities are used,

the order of the algorithm degrades and the performance becomes unaccept-

able. Possibilities are to predict the velocity for the force calculation and

correct it afterwards, or solve the velocity iteratively, requiring more than

one force evaluation per step. Lowe (1999) has devised an alternative algo-

rithm which adds a thermostat much like Andersen™s velocity rescaling. The

equations of motion are integrated with the velocity-Verlet scheme, but in

addition randomly selected pairs of particles exchange their relative velocity

for a sample drawn from a Maxwellian distribution, in such a way that mo-

mentum (and angular momentum) is conserved. This solves problems with

temperature drift that otherwise occur.

11.4 Applications

Applications have been published in many di¬erent ¬elds, such as polymer

rheology (Schlijper et al., 1995), rheology of colloidal suspensions (Koelman

and Hoogerbrugge, 1993; Boek et al., 1997), ¬‚ow of DNA molecules in micro-

channels (Fan et al., 2003). The method can be applied to mixtures and

to microphase separation (Groot and Warren, 1997). Figure 11.1 shows the

impact and subsequent coalescence of two drops of liquid moving towards

each other in a liquid environment in which the drops don™t mix. The sim-

ulation comprises 3.75 million DPD particles and was carried out by Florin

O. Iancu, University of Delft, the Netherlands (Iancu, 2005). The collision

2

is characterized by the dimensionless Weber number W e = ρDUr /“, where

ρ is the density, D the diameter of the drops, Ur the relative velocity of

310 Dissipative particle dynamics

the drops just before impact and “ the interfacial surface tension. At low

Weber numbers the drops bounce o¬ each other without mixing, and at

higher Weber numbers they coalesce after impact. The DPD parameters in

3

this case were a density of 10 particles per rc (rc being the cut-o¬ range

of the repulsive potential) and a repulsion parameter a (see (11.5)) of 14.5

mutually between drop particles or between environment particles, but 41.5

between drop particles and environment particles. This choice leads to an

2

interfacial tension “ = 28 (units of kB T /rc ). When the collision is o¬-center,

the drops elongate before they coalesce and may even split up afterwards

with formation of small satellite droplets.

11.4 Applications 311

Figure 11.1 Impact and coalescence of two liquid drops moving towards each other

in an inmiscible liquid environment at Weber™s number of about 4, simulated with

3.75 million DPD particles (courtesy of Dr F.O. Iancu, University of Delft, the

Netherlands)

Part II

Physical and Theoretical Concepts

12