University of Groningen
Efficient Algorithms for Langevin and DPD Dynamics
Goga, N.; Rzepiela, A. J.; de Vries, A. H.; Marrink, S. J.; Berendsen, H. J. C.
Published in:
Journal of Chemical Theory and Computation
DOI:
10.1021/ct3000876
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from
it. Please check the document version below.
Document Version
Publisher's PDF, also known as Version of record
Publication date:
2012
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Goga, N., Rzepiela, A. J., de Vries, A. H., Marrink, S. J., & Berendsen, H. J. C. (2012). Efficient Algorithms
for Langevin and DPD Dynamics.
Journal of Chemical Theory and Computation
,
8
(10), 3637-3649.
https://doi.org/10.1021/ct3000876
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the
author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
The publication may also be distributed here under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license.
More information can be found on the University of Groningen website: https://www.rug.nl/library/open-access/self-archiving-pure/taverne-
amendment.
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately
and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the
number of authors shown on this cover page is limited to 10 maximum.
Download date: 25-08-2022
Efficient Algorithms for Langevin and DPD Dynamics
N. Goga, A. J. Rzepiela, A. H. de Vries, S. J. Marrink, and H. J. C. Berendsen*
Groningen Biomolecular Sciences and Biotechnology Institute, Zernike Institute for Advanced Materials,
University of Groningen, Nijenborgh 7, 9747 AG Groningen, The Netherlands
ABSTRACT: In this article, we present several algorithms for stochastic dynamics, including Langevin dynamics and different
variants of Dissipative Particle Dynamics (DPD), applicable to systems with or without constraints. The algorithms are based on
the impulsive application of friction and noise, thus avoiding the computational complexity of algorithms that apply continuous
friction and noise. Simulation results on thermostat strength and diffusion properties for ideal gas, coarse-grained (MARTINI)
water, and constrained atomic (SPC/E) water systems are discussed. We show that the measured thermal relaxation rates agree
well with theoretical predictions. The influence of various parameters on the diffusion coefficient is discussed.
1. INTRODUCTION
The purpose of this article is to present simple and effi cient
algorithms for stochastic dynamics, in particular for simple
Langevin dynamics and for various variants of the pairwise
DPD (Dissipative Particle Dynamics) thermostat. The latter are
Galilean invariant; i.e., the motion is the same in a coordinate
system that moves with constant velocity, which is equivalent to
the conservation of total linear momentum. In contrast, the
simple Langevin dynamics will damp all velocities, including a
bulk flow component. Galilean invariance is essential for
simulations that are required to follow Navier−Stokes behavior
in the macroscopic limit but is irrelevant for systems at rest. We
require of such algorithms that they have the same order of
accuracy as the Verlet (or leap-frog or velocity-Verlet)
algorithms that are used in the frictionless case; in the limit
of small friction, they should be equivalent to these algorithms.
They should be valid for any value of the friction constant, not
only in the limit of small friction. They should maintain the
canonical distribution and act as a proper thermostat: i.e., the
temperature derived from the average kinetic energy should
converge to the reference temperature used in the algorithm.
The design of stochastic algorithms for molecular simulation
was an important subject of research in the 1980s for van
Gunsteren et al.
1−4
Van Gunsteren and Berendsen's 1988
paper
4
describes a sophisticated algorithm that fully maintains
the accuracy of the Verlet algorith m by integrating the
stochastic term over the time step. This algorithm is still
standard in both the GROMOS and Gromacs simulation
packages.
Stochastic thermostats have an advantage over the usual
global thermostats that they maintain the correct canonical
distribution and show a robust first-order decay of temperature
deviations toward the reference temperature. Global thermo-
stats of the weak-coupling type
5
show a convenient first-order
decay of temperature de viations but do not maintain a
canonical distribution; in fact, the resulting distribution is in
between a canonical and a microcanonical distribution,
depending on the coupling constant used.
6
Moreover, it is
known that such thermostats may cause an uneven distribution
of kinetic energy among different collective degrees of freedom,
resulting in overheating of collective degrees of freedom that
are only weakly coupled to the other degrees of freedom at
the expense of cooling of these other degrees of freedom (the
flying ice-cube effect
7
). The velocity rescaling thermostat of Bussi
et al.
8
is a global weak-coupling thermostat in which the velocities
are scaled such that the temperature follows a stochastic first-order
process that ensures a canonical distribution for the global
temperature. This solves the problem of an unknown ensemble,
but it is not known if it also solves the flying ice-cube effect. Global
thermostats of the extended-system type, such as the Nose
-Hoover
thermostat,
9
will maintain a canonical distribution in configura-
tional space but show strong oscillatory behavior of the
temperature deviation as a result of the order of the differential
equation related to the system’s extension. The Nose
-Hoover
chain thermostat,
10
using a cascade of thermostats, alleviates but
does not solve this problem. See Berendsen,
11
pp 194−204, for a
comparative discussion of thermostats.
The simple Langevin equations of motion
11
(i.e., Markovian
and no frictional coupling between degrees of freedom) are for
the ith degree of freedom in a Cartesian system of coordinates:
=
x
tvt() ()
ii
(1)
γη
=− +
v
tat vt t() () () ()
ii
i
i
i
(2)
where a
i
(t) is the acceleration of the ith degree of freedom:
=
a
t
Ft
m
r
()
({ ( )})
i
i
i
(3)
and η
i
(t) is a random process with zero mean, no correlation
with any past or present x or v, and with an autocorrelation
function
ηη τ γδτδ
⟨
+⟩=tt
kT
m
() ( ) 2 ()
ij
i
i
ij
Bref
(4)
Here, γ
i
is the friction rate applied to the ith degree of freedom.
[We reserve the name “friction constant” for the proportionality
Special Issue: Wilfred F. van Gunsteren Festschrift
Received: February 1, 2012
Published: June 13, 2012
Article
pubs.acs.org/JCTC
© 2012 American Chemical Society 3637 dx.doi.org/10.1021/ct3000876 | J. Chem. Theory Comput. 2012, 8, 3637−3649
constant between friction force and velocity, which equals the
mass of the particle times the friction rate.]
This equationand its DPD equivalent,
12−14
which applies
friction and noise to velocity differences between particle pairs
is a stochastic differential equation. We note that eq 2 is a
formal equation that is mathematically incorrect because v(t)is
not differentiable; it should be implemented as a difference
equation over a small time step h:
γγξ
Δ
=− +vath vth kThm() () 2 /
ii
i
i
i
iBref
(5)
for small h, where ξ is a normally distributed random number
with unit variance. The last term is obtained by integrating the
noise term in eq 2 over the time step h.
It can be shown that the velocity distribution will converge to
a Maxwellian distribution at the reference temperature T
ref
. The
correct behavior as a thermostat is implicit in the equations and
will be guaranteed only if the equations are solved exactly.
Algorithms to solve such equations in finite time steps are
only exact in the limit of small h because of the position-
dependent force; however, the stochastic terms can be
integrated exactly over any time step. Algorithms
2,4,15
that
integrate the stochastic equations of motion over a time step
become very complex, requiring the sampling of two random
variables from a bivariate distribution.
The same applies to most algorithms used for DPD. While
for Langevin dynamics integration over the frictional and
stochastic terms can be carried out over a time step, for the
pairwise DPD friction and noise, this is not possible. Thus, for
DPD, an additional complicating factor arises: friction and
noise that are needed to update the velocity depend on the
velocity itself. Accurate algorithms require iterative solutions for
the velocity update. Various algorithms based on the velocity-
Verlet scheme have been compared by Nikunen et al.
16
Adifferent approach, pioneered by Peters
17
for the case of
DPD, leads to simpler and still correct algorithms. The
principle is to consider the physical process as a sequence of
a Hamiltonian evolution over one time step, followed by an
impulsive action of friction and noise. The latter modifies the
velocities without advancing the time. Thus, the evolution in
phase space is the approximate application of the Liouville
operator over a time step, followed by a transformation defined
by the impulsive friction and noise. If it can be proven that this
impulsive action also leads to convergence to a canonical (i.e.,
Maxwellian) distribution at the reference temperature, this
physical process is equally val id to achieve our goal of
introducing an effective thermostat. The energy and
momentum transfer implicit in the impulsive friction will
influence transport properties in a controllable way. The
behavior as a thermostat will be robust, while the algorithm
remains very straightforward and simple to implement. The
principal difference wit h the usual stochastic differential
equation is that the time evolution of the system is not
described by a single stochastic differential equation but by a
sequential application of a Hamiltonian evolution over a time
step and an impulsive stochastic action on the velocities. Both
steps should conserve the canonical distribution in phase space.
This approach is reminiscent of the principle of the Andersen
thermostat,
18
which applies impulsive redistributions to particle
velocities with a probability Γh from a Maxwellian distribution.
The same thermostatting method was introduced for DPD by
Lowe.
19
The Andersen and Lowe thermostats do not reduce
velocities with an adjustable friction but reduce the velocity
completely to zero, followed by a noise term equal to a sampling
from the full Maxwellian distribution. The impulsive change
introduced this way is locally (in time and space) quite large and
disrupts the smooth evolution of the trajectory. The impulsive
friction and noise as applied by Peters spreads the changes much
more smoothly over time and space by applying small impulsive
changes to every particle at every step.
In the following, we extend Peters’ DPD-type impulsive
friction and noise (which applies relative velocity changes in the
interparticle direction only) to more general types. These
include the traditional Langevin dynamics, which does not
conserve linear momentum, and pairwise Galilean-invariant
friction and noise acting on velocity differences. The latter can
be applied to all three velocity components, or restricted to
components either parallel or perpendicular to the interparticle
direction. For each case, we investigate the effective thermostat
strength implied by the impulsive application of friction and
noise. We also investigate whether the pairwise application can
be made more efficient by restricting friction and noise to one
or a few neighbors, rather than all neighbors within a given
range. For the ideal gas, we derive the diffusion coefficient,
imposed by friction and noise. While Peters implements friction
and noise in the velocity-Verlet integration scheme, we employ
the popular leap-frog integration scheme. In all cases, we
consider the implementation when the system contains
holonomic constraints, and we consider the consequences for
the computation of pressure.
This article is organized in the following way. In the
following two sections, the impulsive scheme is discussed for
Langevin dynamics using the leap-frog algorithm, for systems
without constraints (section 2) and with constraints (section 3).
Section 4 presents its application to pairwise interactions. In
section 5, an effective average friction rate is defined that serves
as a comparable friction rate measure, with a comparable
influence on diffusion, for any type of impulsive friction
imposed on the system. Section 6 describes the simulation
details and reports the computational efficiency of the different
methods; section 7 gives the results of various tests of the
algorithms on different molecular systems (ideal gas, coarse-
grained water without constraints, and atomic water with
constraints). Both t hermostat and diffusion behavior are
considered. Section 8 discusses the results and summarizes
the conclusions.
2. THE IMPULSIVE LANGEVIN LEAP-FROG
ALGORITHM FOR SYSTEMS WITHOUT
CONSTRAINTS
Consider a system of n particles with 3n degrees of freedom,
and consider every degree of freedom separately. Assume v(t −
(1/2)h), x(t), a nd F(t)=ma are the known velocity,
coordinate, and force components, and a is the acceleration
at time t of that degree of freedom. The impulsive Langevin
extension of the leap-frog algorithm then reads as follows:
For all degrees of freedom, do:
=− +
⎜⎟
⎛
⎝
⎞
⎠
vvt h ah
1
.
1
2
ξΔ=− + −vfvf fkTm
2
.(2)(/)
Bref
+= ++Δ
⎜⎟
⎛
⎝
⎞
⎠
xt h xt v vh
3
.( ) ()
1
2
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct3000876 | J. Chem. Theory Comput. 2012, 8, 3637−36493638
+=+Δ
⎜⎟
⎛
⎝
⎞
⎠
vt h v
v4
.
1
2
Here, step 1 is the usual MD velocity-update of the leap-frog
scheme; step 2 is the impulsive application of friction (reducing
the velocity by a fraction f:0≤ f ≤ 1) and noise (ξ is a random
sample from a normal distribution). Step 3 updates the
coordinates, taking into account that Δv is applied only
between t + (1/2)h and t + h: in fact, step 3 can be considered
as two half steps (see Figure 1):
+=+
⎜⎟
⎛
⎝
⎞
⎠
xt h xt v h
3
a.
1
2
()
1
2
+= + ++Δ
⎜⎟
⎛
⎝
⎞
⎠
xt h xt h v v h
3
b. ( )
1
2
()
1
2
Step 4 assigns the modified velocity to the velocity at the end of
the time step. It is irrelevant whether these steps are carried out
sequentially per degree of freedom, sequentially per 3-D vector
per particle, or each performed on a 3n-D vector over all
degrees of freedom.
The variance of the noise term is chosen such that the
variance of the velocity
⟨
+Δ ⟩= − ⟨ ⟩+ −vv fvf f
kT
m
()(1) (2)
222
Bref
(6)
tends to the stationary value of k
B
T
ref
/m. This is easily seen
from eq 6 by substituting k
B
T
ref
/m for each of the mean
squared velocities.
In this algorithm, it is assumed that all degrees of freedom are
subjected to the same friction and noise at every time step; this
is just a convenient scheme that could be replaced by other
variants, e.g., different f ’s for different particles, or application of
friction and noise to a (randomly) selected subset at every step.
Note that the limiting case f = 1 completely removes the
velocity and replaces it by a sample from a Maxwellian
distribution. Thus, if the impulsive friction and noise is applied
with f = 1 and with a probability Γh per particle per step, the
Andersen thermostat is recovered. A smoothed Andersen
thermostat with the same average velocity reduction factor will
be obtained by applying the impulsive friction and noise with
f = Γh every step to every degree of freedom.
The algorithm is expected to be robust, in the sense that the
impulsive term is exact, independent of the time step used.
There is a lot of freedom of choice in the way the impulsive
term is applied: the damping factor f may be applied to a
random selection of particles and may differ for different
particles. The application to velocity differences as in DPD is
straightforward (see section 4).
It is obvious that the new velocity of any particle does not
have the same direction and magnitude as its old velocity. As
the random changes are uncorrelated, the total momentum is
not conserved and neither is the total energy conserved.
However, the average kinetic energy and hence the temperature
will be stable: they tend toward the values determined by the
reference temperature.
2.1. Is the Velocity Distribution Canonical? In order to
judge the acceptability of the proposed procedure, we ask the
following questions: Assume the velocity distribution before the
impulse is ρ
0
(v): (a) What will be the distribution ρ
1
(v + Δv)
after the impulse? (b) What is the stationary distribution?
After the impulse, the distribution ρ
1
(v + Δv) is the
convolution of the original distribution ρ
0
(v) and the Gaussian
distribution of the random term Δv + fv, which has a variance
of f(2 − f)k
B
T
ref
:
∫
ρ
πρ=−
×−
−−
−
−
−∞
∞
⎡
⎣
⎢
⎤
⎦
⎥
wffkT vv
wfv
ffkT
() [2 (2 ) ] d ()
exp
{(1)}
2(2 )
1
Bref
1/2
0
2
Bref
(7)
where w = v + Δv. This is the answer to question a. The answer
to question b is found by inserting the canonical distribution
for v:
ρ
π=−
−
⎡
⎣
⎢
⎤
⎦
⎥
vkT
v
kT
() [2 ] exp
2
0
Bref
1/2
2
Bref
(8)
into eq 7. Carrying out the integration over v,wefind that
ρ
1
(w) is exactly equal to the same canonical distribution ρ
0
. So,
ρ
0
(v) is the stationary distrib ution. Thus, the impulsive
application of friction and noise not only preserves the variance
(as designed) but preserves the complete canonical distribution
as well.
2.2. How Does the Temperature Behave with Time? A
good thermostat should force a deviation from the reference
temperature back to zero. How does the impulsive Langevin
thermostat behave in this respect?
Consider one dimension. The temperature is given by
=⟨⟩
T
m
k
v
B
2
(9)
The energy change ΔE
1
resulting from a single application of
friction and noise to one degree of freedom is
Δ
=+Δ−Emvv mv
1
2
()
1
2
1
22
(10)
Using eqs 6 and 9, this rewrites to
Δ
=− −EffkTT
1
2
(2 ) ( )
1Bref
(11)
the total energy change per time step h is the sum over all one-
dimensional frictional events that occur per time step:
∑
Δ
=−−EkffTT
1
2
(2 )( )
tot B ref
(12)
The energy change is initially supplied to the kinetic energy of
the system, thus changing the temperature. However, when the
Figure 1. (Top) Traditional leap-frog scheme. (Bottom) Leap-frog
scheme with impulsive phase.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct3000876 | J. Chem. Theory Comput. 2012, 8, 3637−36493639
rate of change is small, the energy change will be distributed
over kinetic and potential energy. The temperature change is
then determined by the total heat capacity C
V
of the system:
Δ
=
Δ
T
E
C
V
tot
(13)
yielding a differential equation for the time dependence of the
temperature:
=
∑
−
−
dT
dt
k
ff
Ch
TT
1
2
(2 )
()
V
Bref
(14)
In the case of a three-dimensional application to N particles, eq
14 has the form
=
−
−
dT
dt
k
c
ff
h
TT
3
2
(2 )
()
V
B
ref
(15)
where
=c
C
N
V
V
(16)
is the specific heat per particle.
Equation 15 shows that any deviation from the reference
temperature will decay to zero according to a first-order kinetic
process
=−
dT
dt
kT T()
th ref
(17)
with rate constant
=
−
k
k
c
ff
h
3
2
(2 )
V
th
B
(18)
Alternatively, the decay can be characterized by a time constant
τ
T
=1/k
th
. Note that for an ideal gas, c
V
= (3/2)k
B
, reducing the
left fraction in eq 18 to 1; for atomic fluids, this fraction is
usually 2 to 3 times smaller.
The thermal rate constant can be expressed in an effective
friction rate γ
eff
,defined by the continuous friction rate that
would reduce the velocity per time step h by a fraction f:
γ
=− −
h
f
1
ln(1 )
eff
def
(19)
yielding
γ
γ=
−−
≈k
k
c
h
h
k
c
3
2
[1 exp( 2 )]
3
2
2
VV
th
B
eff
B
eff
(20)
the latter value being a good approximation for small γh.
Thus, the thermostat is robust: the system temperature
automatically decays to the reference temperature, and the
velocity distribution evolves into the proper canonical
(Maxwellian) distribution.
Slow and Fast Thermostats. We note that this rate equation
is valid when the rate constant of the thermostat is smaller than
the rate of exchange between the kinetic and potential energy
of the system, which gives the system time to equilibrate
between kinetic and potential degrees of fzreedom. Usually,
this condition is fulfilled. If, on the other hand, k
th
is much
larger (“fast” thermostat s), C
V
in eq 13 should be replaced by
E
kin
/T = (1/2)k
B
n
dof
, where n
dof
is the number of degrees of
freedom in the system. The equivalent of eq 14 for fast
thermostats is
=
∑
−
−
dT
dt
ff
nh
TT
(2 )
()
dof
ref
(21)
where the sum is taken over all one-dimensional frictional
events that occur per time step. This means that for a fast
thermostat, c
V
in eq 20 should be replaced by its ideal-gas value
3k
B
/2.
3. IMPULSIVE LANGEVIN ALGORITHMS FOR
SYSTEMS WITH CONSTRAINTS
Here, we give the algorithm for simple Langevin dynamics with
constraints (section 3.3). For clarity, we first give the leap-frog
algorithms for frictionless systems without constraints (section
3.1)already mentioned in section 2and with constraints
(section 3.2).
3.1. No Constraints, No Friction. Consider a conservative
system with n particles (3n degrees of freedom). Assume v(t −
(1/2)h) and x(t) are the known velocity and coordinate vectors
(length 3n) at the beginning of the time step. The leap-frog
algorithm then advances one time step as follows (i = 1, ..., 3n):
1. Compute forces F = −∇V(x(t)) and compute virial using
x(t).
2. Advance velocities ∀ i: v
i
(t + (1/2) h)=v
i
(t − (1/2) h)+
(F
i
/m
i
)h (optional: adjust ∀i, v
i
(t + (1/2)h) according to
global thermostat).
3. Advance positions ∀i: x
i
(t + h)=x
i
(t)+v
i
(t + (1/2)h)h
3.2. Constraints, No Friction. Constraints are satisfied by
calling a routine constr (SHAKE,
20
LINCS,
21
SETTLE
22
) that
makes corrections to an unconstrained configuration x
u
to yield
a configuration x
c
that satisfies the set of m holonomic
constraints σ
k
(x)=0,k = 1, ..., m. The displacements are in the
direction of a reference configuration x
ref
:
=
x
xxconstr( ; )
curef
(22)
The routine constr also computes the constraint forces and
evaluates the contribution to the virial (using the reference
positions) due to the constraint forces. The algorithm of
section 3.1 is now modified to:
1. Compute unconstrained forces F
u
= −▽V(x(t)) and
compute virial using x(t).
2. Advance velocities ∀i: v
i
= v
i
(t − (1/2)h)+(F
i
u
/m
i
)h
(optional: adjust ∀i, v
i
according to global thermostat).
3. Advance positions ∀i: x
i
= x
i
(t)+v
i
h.
4. Apply constraints to new positions: x(t + h) = constr
(x; x(t)) and add constraint contribution to virial.
5. Reconstruct velocities ∀i: v
i
(t + 1/2) = [x
i
(t + h) −
x
i
(t)]/h.
Steps 4 and 5 ensure that the motion takes place on the
hypersurface in configurational space, defined by the con-
straints. Any deviation from this hypersurface due to force
components in the constraint directions (steps 2 and 3) is
removed by the projection onto the constraint hypersurface in
steps 4 and 5. We note that the reconstruction of velocities
from coordinate differences introduces a noticeable, but
avoidable, error when single-precision arithmetic is employed.
23
3.3. Constraints Plus Impulsive Friction and Noise. We
now modify the algorithm of section 2 to include an impulsive
friction and noise:
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct3000876 | J. Chem. Theory Comput. 2012, 8, 3637−36493640