Splitting for dissipative particle dynamics
Shardlow, Tony
2003
MIMS EPrint: 2006.261
Manchester Institute for Mathematical Sciences
School of Mathematics
The University of Manchester
Reports available from: http://eprints.maths.manchester.ac.uk/
And by contacting: The MIMS Secretary
School of Mathematics
The University of Manchester
Manchester, M13 9PL, UK
ISSN 1749-9097
SPLITTING FOR DISSIPATIVE PARTICLE DYNAMICS
∗
TONY SHARDLOW
†
SIAM J. S
CI. COMPUT.
c
2003 Society for Industrial and Applied Mathematics
Vol. 24, No. 4, pp. 1267–1282
Abstract. We study numerical methods for dissipative particle dynamics, a system of stochastic
differential equations for simulating particles interacting pairwise according to a soft potential at
constant temperature where the total momentum is conserved. We introduce splitting methods and
examine the behavior of these methods experimentally. The performance of the methods, particularly
temperature control, is compared to the modified velocity Verlet method used in many previous
papers.
Key words. stochastic differential equations, numerical methods, molecular dynamics
AMS subject classifications. 60H10, 82C80
PII. S1064827501392879
1. Introduction. Dissipative particle dynamics (DPD) was first introduced by
Hoogerbrugge and Koelman in [9] to simulate complex hydrodynamic behavior. The
technique is based on simulation of particles that evolve according to a system of SDEs
(see (1) below). DPD particles should not be thought of as individual molecules; the
particles represent mesoscopic groups of fluid molecules that interact at short range
and according to a soft potential. This coarse graining provides a model, thought to be
realistic in some regimes, that is convenient computationally, allowing for large time
steps in simulations and for the exploration of long time scales. Recent applications of
DPD include polymer simulation, spinodal decomposition, and suspension rheology
[21, 12, 7, 4, 1]. Theoretical investigations of DPD have focused on understanding the
physics (see [5, 24, 8, 18]). The issue we tackle is the numerical simulation of a DPD
fluid.
DPD is described by the following system of SDEs. Consider N particles with
positions q
i
and momenta p
i
for i =1,...,N evolving in dimension d:
dq
i
= p
i
dt,
dp
i
= −
j=i
a
ij
V
(q
ij
)
ˆ
q
ij
dt − γ
j=i
w
D
(q
ij
)(
ˆ
q
ij
· p
ij
)
ˆ
q
ij
dt + σ
j=i
w
R
(q
ij
)
ˆ
q
ij
dβ
ij
(t).
(1)
We take the simplest case and prescribe periodic boundary conditions on the posi-
tions q
i
in the domain [0,L]
d
(see [18] for a discussion of other boundary conditions).
The relative positions and momenta are denoted by q
ij
and p
ij
and the unit direction
from q
j
to q
i
by
ˆ
q
ij
and the length of q
ij
by q
ij
. We define
ˆ
q
ij
= 0 for q
ij
= 0. The
pair potential
V (r)=
1
2
(1 − r/r
c
)
2
,r<r
c
;
0 otherwise,
where r
c
is the radius of interaction. The parameters a
ij
are positive and symmetric
(a
ij
= a
ji
≥ 0) and describe the strength of repulsion between particles i and j. For
∗
Received by the editors July 26, 2001; accepted for publication (in revised form) July 30, 2002;
published electronically February 6, 2003.
http://www.siam.org/journals/sisc/24-4/39287.html
†
Department of Mathematics, Oxford Road, University of Manchester, Manchester M13 9PL, UK
(shardlow@maths.man.ac.uk).
1267
1268 TONY SHARDLOW
i<j, β
ij
are IID Brownian motions and for i>j, β
ij
= −β
ji
. We will use E to
denote averages with respect to realizations of the Brownian motion.
The functions w
D
and w
R
describe the dissipative and random forces,
w
D
(r)=w
R
(r)
2
=
(1 − r/r
c
)
2
,r<r
c
;
0 otherwise.
The strength of the dissipation is parameterized by γ and of the noise by σ.
The mathematical theory of solutions of SDEs is discussed in [10], for example.
A rigorous study of existence, uniqueness, and regularity is beyond the scope of the
present paper. The difficulty in (1) is that neither drift nor diffusion coefficients are
continuous at q
ij
= 0, and this could lead to nonuniqueness of solutions (even of weak
solutions) for initial data with q
ij
= p
ij
= 0, some i = j. This can be understood by
considering
dq = p dt, dp =
ˆ
q dβ(t), q(0) = Q, p(0) = P,
for a standard Brownian motion β(t) and
ˆ
q := q/q (where ·denotes the standard
Euclidean norm). For initial data Q = P = 0, no matter how
ˆ
q is defined at q =0,
the following
q(t)=
ˆ
e
t
0
β(s) ds, p(t)=
ˆ
e β(t)
is a solution on a random time interval [0,τ], where τ>0 and
ˆ
e ∈ R
d
is a unit
vector. For dimension d>1, this gives rise to multiple weak solutions. Even if we
assume the initial data is nonzero, there is a positive probability of reaching the zero
state. An extra condition is needed at q = p = 0 to specify the solution uniquely.
To the author’s knowledge, a rigorous analysis of this issue for DPD has not been
undertaken. We assume the existence of a Markov process ({q
i
}, {p
i
}) on the space
of configurations, which converges to the Gibbs canonical distribution (defined below)
and whose generator A satisfies
A =
i
p
i
·∇
q
i
−
i=j
a
ij
V
(q
ij
)
ˆ
q
ij
·∇
p
i
− γw
D
(q
ij
)(
ˆ
q
ij
· p
ij
)
ˆ
q
ij
·∇
p
i
+
1
2
σ
2
i<j
w
R
(q
ij
)
2
ˆ
q
T
ij
(∇
p
i
−∇
p
j
)
2
ˆ
q
ij
(2)
and has domain D(A), a subset of the space of continuous functions on R
dN
× R
dN
that have period L in the spatial components.
The fluctuation dissipation theorem applies to (1): Suppose that σ
2
=2γk
B
T ,
where k
B
is Boltzmann’s constant and T is the equilibrium temperature. If
H({q
i
}, {p
i
}):=
1
2
i
p
i
2
+
1
2
i=j
a
ij
V (q
ij
),
then the Gibbs canonical distribution
µ({q
i
}, {p
i
}):=
1
Z
exp
−H({q
i
}, {p
i
− p})
k
B
T
,(3)
SPLITTING FOR DISSIPATIVE PARTICLE DYNAMICS 1269
where Z is chosen so that µ is the density of a probability measure and p is the
average momentum (over the N particles). As shown in [5], µ is an invariant measure
of the SDEs (1) and
1
d
1
N
E
µ
N
i=1
p
i
2
−
p
2
d
= k
B
T.(4)
(E
µ
denotes expectation with respect to µ.)
The main emphasis of this paper is numerical methods for solving (1). Many
papers have used a modified version of the velocity Verlet method (defined in sec-
tion 2.1) to do this [8, 5, 6, 24]. The velocity Verlet scheme itself is a method for
solving the Hamiltonian equations for H(q, p)=
1
2
p
2
+ V (q) and consists of iterating
the following for a time step ∆t:
p
n+
1
2
= p
n
−
1
2
∆tV
(q
n
),q
n+1
= q
n
+∆tp
n+
1
2
,p
n+1
= p
n+
1
2
−
1
2
∆tV
(q
n+1
).
This scheme is second order, symplectic, and will conserve linear and angular mo-
mentum if this holds for the Hamiltonian system. The force function, which is the
most expensive part of the iteration, is computed only once per time step. For further
discussion, see, for example, [19]. The modified Verlet method accounts for the noise
and dissipation terms in (1). In DPD linear momentum is conserved; angular momen-
tum is conserved when q
i
evolve on the spatial domain R
d
rather than the periodic
domain [0,L]
d
. The modified Verlet scheme inherits these properties. The method
gives satisfactory results in many applications, but the time step must be small for
the method to be stable. This becomes a more severe problem when the dissipation
parameter γ or the density of particles is large. Marsh and Yeomans [13] discuss this
issue in some detail and derive critical temperatures and densities for the method to
be stable.
A large density or dissipation in (1) gives rise to stiffness, which is usually dealt
with by applying an implicit numerical method. This approach has been tried by
Pagonabarraga, Hagen, and Frenkel [15], who investigate an implicit method but find
it expensive to run. Another implicit method was attempted during the writing of
the present paper, based on replacing the relative momenta p
n
ij
at time step n with
a semi-implicit momenta p
n
i
− p
n+1
j
. This leads to an implicit numerical method
where only block diagonal matrices need be inverted. This method is cheap to run
but destroys the underlying conservation of linear and angular momentum, and poor
behavior is observed [22].
Splitting is a widely used and important technique for solving differential equa-
tions; one important application to the Langevin equations is described in [20]. The
main contribution of this paper is an implicit method based on splitting the vector
field into a sum of conservative terms and pairwise fluctuation-dissipation terms. We
solve the conservative system defined by the Hamilitonian H (and the boundary con-
ditions) using the velocity Verlet method. Following this, we solve an SDE for the
fluctuation and dissipation between particles i and j. We choose an implicit method
that conserves the invariance of the momenta and that may be solved efficiently. This
is possible because of the simple structure of the pair equations. Looping over all
possible pairs, we can find an approximation for the equation as a whole. There are
a number of ways of combining the solvers for the split systems.
To introduce the splitting method, we restrict attention to weak convergence.
That is, for a functional φ (such as temperature) on the set of configurations {q
i
}, {p
i
},
1270 TONY SHARDLOW
we want to generate approximations {q
n
i
}, {p
n
i
} such that the averages (over the set of
Brownian paths) of φ({q
n
i
}, {p
n
i
}) converge to the average of φ({q
i
(n∆t)}, {p
i
(n∆t)}).
Weak convergence is of most interest in molecular dynamics and gives convergence
of approximations to time correlations and, under ergodic assumptions, to averages
with respect to the canonical distribution.
For weak convergence, we approximate the action of the semigroup, e
At
, with
generator A on the functional φ. This allows us to introduce splitting methods by
writing the generator A = A
1
+ A
2
. We will consider first order splitting,
e
An∆t
=(e
A
1
∆t
e
A
2
∆t
)
n
+ O(n∆t
2
) ,
and the following second order splitting,
e
An∆t
=(e
A
1
∆t/2
e
A
2
∆t
e
A
1
∆t/2
)
n
+ O(n∆t
3
) .
The first order splitting was introduced by Trotter; the second order splitting was
introduced by Strang [23] for numerical approximation. Numerical analysts exploit
these relations by approximating e
A
1
∆t
, e
A
2
∆t
separately and using the above rela-
tions to generate a consistent approximation to e
A∆t
. The splitting formula may be
applied recursively to deal with a generator A split into multiple components. A
rigorous statement of the convergence properties depends on the regularity of the
test function φ and the generators A
1,2
and is not described herewith. We do define
split methods for (1) in section 2.2 and demonstrate that the methods are valuable
for DPD. It is not trivial to establish regularity of the generator in (1), because the
diffusion is not uniformly elliptic nor are the coefficients smooth. Further discussion
of splitting methods for SDEs includes [17, 3, 14].
DPD is used to simulate fluids in thermal equilibrium, where the solution of (1)
is evolving in the canonical distribution (3). The temperature with respect to the
canonical distribution is known and this provides a convenient way to evaluate the
numerical methods. Assuming ergodicity of DPD, we expect time averages of the
temperature to converge to the average given in (4). To evaluate the methods, we
compute the time average
1
d
1
N − 1
1
T
1
− T
0
T
1
/∆t
n=T
0
/∆t
N
i=1
p
n
i
2
∆t,(5)
where T
0
is chosen to allow the system to reach equilibrium and T
1
is chosen large to
reduce the variance in the result. Notice that the average is computed by dividing by
N − 1 (the number of degrees of freedom) rather than N (the number of particles).
We will take initial conditions with zero total linear momentum, in which case the
computed temperature (5) will be compared to k
B
T .
After defining the methods in section 2, results are presented in section 3 for pa-
rameter values that test the modified Verlet scheme with large density and dissipation.
The results indicate the advantage of using the first order splitting method: 1% tem-
perature control is achieved with ∆t =0.04 in section 3.1 for the first order splitting
scheme whilst the modified Verlet method needs ∆t =0.02. Though the modified
Verlet method is less computationally expensive per time step, it is faster to use the
splitting technique to get 1% temperature control (the first order splitting method
with ∆t =0.04 is 40% faster than the modified Verlet method with ∆t =0.02). The
improvements are more dramatic when the dissipation (see section 3.2) or density is
increased (see section 3.3).
