To simulate the dynamics of polymers or lipids in the membrane we use the DPD (Dissipative Particle Dynamics)
DPD Dissipative Particle Dynamics
simulation method. The original work from [Koeleman 12
dissipative particle dynamics CHAPTER 1. MODEL and Hoogerbrugge (1993)] was successive upgraded by [Pagonabarraga and Frenkel
(2001)] into a MDPD (Multi body Dissipative Particle Dynamics) scheme, MDPD Multi-body Dissipative Particle Dynamics
which we explain below.
DPD is a method to integrate the equations of motion, like molecular dynamics sim-ulations but its range of validity is the mesoscopic scale and therefore requires Brown-ian noise which includes hidden interactions with the “ghost” microscopic particle. In contrast to aLangevindescription, DPD uses noise and friction that locally conserve angular and linear momentum. Conservation of the hydrodynamics is important in annealing defects [Gonnella et al. (1997)]. For the DPD thermostat, the equation of
motion is the sum of conservative (c), dissipative (d) and random (r) forces Fcconservative Fddissipative Fr random
mv˙ =Fc+Fd+Fr (1.28)
Every component is limited to a distance interval, ∆L, by a weighting functionw(rij) that depends on the relative distance between the particles. This function is 1 for r= 0 and goes to zero at the cut-off distance r= ∆L. Every force is pair-wise and
thus locally conserves momentum. w(rij)fijstrength of the force
acting between theiandj particle,rijdistance between the iandjparticle.
mv˙ =wc(rij)fcij+wd(rij)fdij+wr(rij)frij (1.29) The weighting function makes the forces soft and permits increase of the time-step of the simulated system [Pastorinoet al. (2007)].
The conservative force depends on the particular system and it is obtained from the derivative of the potential. The dissipative and random forces should obey the fluctuation-dissipation theorem. If we determine a dissipative termγwe have to define the strength,ξ. As shown in [Orlandini (2008)], if we choose an uncorrelated Gaussian random noise with zero average
θijg(t)
= 0
θgij(t)θklg(t0)
= (δikδjl+δilδjk)δ(t−t0) (1.30) the strength terms should satisfy the relation
ξ= r
2γkBT
m∆t (1.31)
In practise, following [D¨unweg and Paul (1991)], we can use a uniform random number generator,θuij, which is faster to compute, instead of a Gaussian. In this case we note that if a Gaussian distribution has a variance ofσ, a uniform distribution equal to 1 between [−σ/2, σ/2] has a variance σ/√
12. Hence, our equation of motion takes the form
mv˙ =fcijwc(rij)−γwd(rij)(vijˆrij)ˆr+
r24γkBT
m∆t θijuwr(rij) (1.32) InGrootandWarren’s work [Groot and Warren (1998)] the DPD simulation method is widely investigated and following their results, we setγ= 0.1 and ∆t= 0.01 in the system’s units: kBT = ∆L=m= 1. On the suggestion of the same work we use the velocityVerletintegration scheme instead of theEuler’s. EspañolandWarren (1995) [Espanol and Warren (1995)] have shown that the dissipative weighting function can be chosen arbitrarily and is connected to the random weighting function by the relation wd(r) = (wr(r))2. In [Pastorinoet al.(2007)] different weighting functions in different polymer system are tested. The particular choice of a weighting function concerns the thermostat of the system and the computational efficiency. We define the number of particles thermostated, i.e. the particles included in the sphere within the cut-off multiplied by the weighting function, as
NT P =ρ0
Z r 0
wr(r)g(r)4πr2dr (1.33)
13
dissipative particle dynamics CHAPTER 1. MODEL where theg(r) in the pair correlation function that in our case we approximate as 1.
The use of the following weighting function wd(r) = (wr(r))2=
1− r
∆L 2
if r <∆L (1.34) provides a fast computing efficiency. In agreement with the suggestions of [H¨omberg and M¨uller (2010)] and [Trominov et al. (2002)], we observe in our simulations that if the cut-off distance include an average number of thermostated particle larger than 3-4 the system conserves temperature keeping a time step of ∆t= 0.01.
The integration of the equations of motion is performed via a velocityVerletalgorithm which is composed of two steps
ri(∆t) =ri(0) + ∆tr˙i(0) + ∆t2
and its derivation is briefly discussed in the appendix.
multibody dissipative particle dynamics
As was noted in [Trominovet al.(2002)], the DPD model, as it was initially formulated, does not take into account the local density of the system. If the conservative forces depend only on the mutual distance between the particles, then the system produces an equation of state of up to second order [Groot and Warren (1998)]. This means that we can not simultaneously map the pressure and the compressibility of the system, because the inverse of compressibility is limited by the value of the pressure
k−1T =ρkBT+vAAρ2= 2P−ρkBT <2P P =ρkBT+vAA
2 ρ2 (1.36) which is not the case for more compressible liquids. [Pagonabarraga and Frenkel (2001)] proposed to redefine the conservative forces dependent upon the local density of the system. This suggests defining the free energy as a functional of the density and to proceed with the calculation of the conservative forces by deriving the free energy with respect to the spatial coordinates.
conservative force
The non-bonded interactions between the beads should come from the negative deriva-tive of the potential.
Fαβ(ri, rj) =−∇rU(rα, rβ) (1.37) where the Greek indices refer to the different types of beads,A, B. The forces depend on the type of beads involved and on their mutual distances. It is now important to show how we obtain a pairwise interaction, which is required by the DPD simulation scheme. The Hamiltonian we use for a one component system is,v2:=vAA,v3:=vAAA If we consider our particles to be points, we write the density as a sum of delta functions
ρ(r) =X
i
δ(r−ri) (1.39)
14
homopolymer melt CHAPTER 1. MODEL Since the product of delta functions is not defined, we rewrite the product of density functions as a delta times a weighting functionw(r)
ρ2(r) := X This weighting function should be normalised by
4π Z ∆L
0
dr r2w(r) = 1 (1.42)
which helps us formulate the Hamiltonian as a sum over all the positions of the
parti-cles,rij :=|ri−rj| rij=|ri−rj| The weighting function,w(r), represents the smoothing operation, which implies that we cannot resolve the system below a certain length scale ∆L. We finally obtain the conservative force by deriving the potential with respect to the coordinateri
Fnbi := −∂ri In this way we obtain a pairwise expression for computing the conservative force, as required by the MDPD simulation method.