On coarse-graining by the Inverse Monte Carlo method: Dissipative Particle Dynamics simulations

made to a precise tool in soft matter modeling

Alexander P. Lyubartsev[a], Mikko Karttunen[b] , Ilpo Vattulainen[c], and Aatto Laaksonen[a]
[a]
Division of Physical Chemistry, Arrhenius Laboratory, Stockholm University, S–106 91, Stockholm, Sweden
[b]
Biophysics and Statistical Mechanics Group, Laboratory of Computational Engineering,
Helsinki University of Technology, P. O. Box 9400, FIN–02015 HUT Helsinki, Finland
[c]
arXiv:cond-mat/0303583v1 [cond-mat.soft] 27 Mar 2003

Laboratory of Physics and Helsinki Institute of Physics,

Helsinki University of Technology, P. O. Box 1100 FIN–02015 HUT, Helsinki, Finland

We present a promising coarse-graining strategy for linking micro- and mesoscales of soft matter systems.
The approach is based on effective pairwise interaction potentials obtained from detailed atomistic molecular
dynamics (MD) simulations, which are then used in coarse-grained dissipative particle dynamics (DPD) simu-
lations. Here, the effective potentials were obtained by applying the Inverse Monte Carlo method [Lyubartsev
& Laaksonen, Phys. Rev. E. 52, 3730 (1995)] on a chosen subset of degrees of freedom described in terms of
radial distribution functions. In our first application of the method, the effective potentials were used in DPD
simulations of aqueous NaCl solutions. With the same computational effort we were able to simulate systems
of one order of magnitude larger as compared to the MD simulations. The results from the MD and DPD
simulations are found to be in excellent agreement.

Keywords: Computer simulations, atomistic force fields, effective potentials, dissipative particle dynamics,
mesoscale modeling, coarse-graining

I. INTRODUCTION have recently introduced a hybrid method to combine particle

and continuum level descriptions. This framework allows a
derivation of the DPD model and in practice uses a Voronoi
Classical molecular dynamics computer simulations with
tessellation based technique for simulations of coarse-grained
site-site pair potentials can currently be run for molecular sys-
particles at different levels. Another very promising approach
tems consisting of the order of 104 atoms, corresponding typi-
was presented by Malevanets and Kapral [7, 8] who proposed
cally to a system size of 50 – 60 Å. Corresponding time scales
a framework to couple a molecular level description with a
using the same approach are typically of the order of 10 ns.
mesoscale treatment of solvent that conserves hydrodynam-
This may be enough for simulations of isotropic liquids of
ics. This approach seems particularly appealing for studying
simple molecules, while for many complex systems the re-
dilute systems with hydrodynamics.
quirements for the length and time scales are much more ex-
tensive. A lot of impetus for multiscale modeling has come from
polymer research due to both fundamental and practical rea-
The restrictions of the detailed all-atom description are es- sons. Akkermans and Briels [9] applied a projector operator
pecially severe in soft matter systems such as biomembranes formalism based approach to coarse-grain polymer chains us-
and polyelectrolytes, although these are just two examples of ing a microscopic model as a starting point. Their approach
the diversity of soft materials. This is mainly due to a wide is appealing and deserves further studies since it accounts for
range of time and length scales associated with these systems, the fluctuating forces at the coarse-grained level. Bolhuis et
and it is obvious that these problems become particularly pro- al. [10, 11] used the Ornstein-Zernike equation [12] with a
nounced in studies of dynamic processes which usually take hypernetted-chain closure to invert the pair correlation func-
place under hydrodynamic conditions. The primary aim in tion in order to derive effective pair potentials for neutral poly-
coarse-graining is to build a bridge between different time mers in a solution. The coarse-grained polymers were repre-
and length scales. Although there is no unique way to do sented as soft colloidal particles. Their approach is based on
this, the coarse-graining procedure always leads to a reduction a theorem proved by Henderson in 1974 [13] that stipulates
of the number of degrees of freedom and thus alleviates the a unique, point-to-point correspondence between pair poten-
computational burden. The challenge is to establish system- tials and radial distribution functions g(r). The importance of
atic coarse-graining schemes that allow one to develop sim- this theorem becomes clear when one notices that for any g(r)
plified model systems in terms of the information extracted under given conditions there is a unique pair potential that in-
from the underlying microscopic systems. Thus far, various cludes corrections from the many-body interactions, and that
approaches have been used to coarse-grain atomic and molec- by definition the pair potential then yields the original g(r).
ular systems. Below, we will briefly discuss some of them. The problem is thus reduced to inverting the g(r) obtained
As a relatively generic approach, the dissipative particle dy- from a microscopic model – how to obtain a coarse-grained
namics (DPD) [1, 2, 3] has been used to study coarse-grained model from that information lies at the heart of the problem
particle and polymeric systems, see e.g. Ref. [4] for a review. of coarse-graining. At the representational level the soft el-
In our approach we use DPD as a thermostating method which lipsoid model for polymer melts by Murat and Kremer [14]
is discussed in more detail in the following sections. To men- is very close to the one of Bolhuis et al. [10, 11], but in-
tion other relatively generic approaches, Flekkøy et al. [5, 6] stead of using the pair correlation function Murat and Kremer
 2

base their approach on separating the free energy into intra- mophore where the time scales of interest are within a few
and interchain parts. In an attempt to model more realistic hundreds of femtoseconds. For studies of systems that are
systems, Reith et al. [15] introduced a coarse-graining pro- characterized by large time and length scales, such as a pro-
cedure that starts from a detailed chemical description of a tein in water solution over tens or hundreds of nanoseconds,
polymer. In order to obtain a coarse-grained approach, a sim- however, atomistic simulation techniques are not very feasi-
plex algorithm based optimization is performed to obtain a ble. Our aim in this section is to discuss ways how this prob-
coarse-grained force field. The structural features produced lem can be resolved by coarse-graining microscopic systems.
by this approach are in good agreement with experimentally In particular, we concentrate on an approach based on solv-
observed ones. For a more complete overview of the current ing the “inverse problem” which yields effective interaction
coarse-graining approaches, see e.g. Ref. [16] and references potentials which in turn can be coupled to mesoscopic sim-
therein. ulation techniques such as the dissipative particle dynamics
Our approach is also based on inverting the radial distri- method discussed below. Full details of the methods can be
bution function but instead of the integral equation based ap- found in Refs. [17] and [19, 20], respectively.
proach of Bolhuis et al. [10, 11] we use the inverse Monte
Carlo procedure of Lyubartsev and Laaksonen [17] which will
be discussed in the following section. This systematic coarse- A. The Inverse Problem in statistical mechanics
graining method is then combined with a momentum conserv-
ing dissipative particle dynamics (DPD) thermostat. Our ap- In the standard statistical mechanical description of soft
proach consists of three conceptually simple steps. First, we matter systems one begins from a formulation of a model,
use a detailed atomistic description and perform MD simu- which is usually given in terms of a Hamiltonian or interac-
lations. From these simulations we compute radial distribu- tion potentials between atoms and molecules. When a Hamil-
tion functions g(r) between different atoms, molecules, or tonian is specified, one can use different numerical methods
molecular groups. Second, we apply the Inverse Monte Carlo to compute canonical averages such as energies, distribution
[17] procedure to invert the radial distribution functions. This functions, and response functions which can be compared
yields effective interaction potentials V eff (r) between the se- with experiments. As one has an explicit expression for the
lected interaction sites. Finally, we employ these effective, canonical averages, their calculation is in principle a fairly
coarse-grained potentials in the DPD algorithm to study the straightforward task, although in many cases it may be com-
large-scale properties of systems in question. Note that since putationally a very demanding one.
DPD satisfies momentum conservation, the present approach
Occasionally one is interested in the solution of the inverse
allows studies of a given system with full hydrodynamics.
problem, i.e., how to deduce information about molecular in-
Preliminary results of this work were recently reported in teractions from the canonical averages which may be known
Ref. [18], in which we presented only the most essential ideas from experiments. More specifically, the question is how to
of this approach. Our aim in this article is to provide one reconstruct the Hamiltonian based on results for the canon-
with a comprehensive discussion of issues related to coarse- ical averages. Many experimental properties can be chosen
graining and mesoscopic simulations in terms of the DPD as a starting point to this end. Among the most important
method. In particular, we discuss in detail how atomistic ones for the liquid state are the radial distribution functions
MD simulations can be coupled to DPD by the Inverse Monte g(r), which may be known from the structure factor measure-
Carlo method. To this end, we use NaCl as our model system. ments using X-ray or neutron scattering [12]. In 1974, Hen-
Even in this simple system, and with a very modest degree of derson proved a theorem [13] stipulating a unique, point-to-
coarse-graining, we obtain a computational speed-up of one point correspondence between pair potentials and radial dis-
order of magnitude. Our procedure is well-defined and truly tribution functions. In essence it states that for a given system
allows easy and controllable tuning of the level of coarse- under given conditions of temperature and density, two pair
graining. Finally, we show through coarse-grained DPD sim- potentials which give rise to the same radial distribution func-
ulations that this approach is physically sound and provides tion cannot differ more than by an additive constant. As this
results in excellent agreement with the MD simulations and constant can be defined from the condition that the interaction
experiments. potential tends to zero at an infinite distance, the potential is
This paper is organized as follows. First, in Sect. II, we uniquely defined. In practice the solution of an inverse prob-
describe the methods: The Inverse Monte Carlo procedure and lem is not a straightforward calculation, however, since g(r)
the DPD algorithm. In Sect. III we present the results and does not provide one with a formal expression for the poten-
finally, in Sect. IV, we discuss the results and some general tial. Therefore, some special techniques are required for the
aspects related to the general applicability of the method. solution of the inverse problem.
In 1988, McGreevy and Pusztai suggested a Reverse Monte
Carlo (RMC) method in which the starting point for a simula-
II. METHODS AND MODELS tion was the radial distribution function [21]. In this method,
Monte Carlo simulations are carried out without any prior
As discussed in the Introduction, detailed atomistic simu- knowledge of the interaction potential to fit g(r) that serves
lation techniques are suitable for studies of microscopic fea- as an input. The set of configurations obtained in RMC may
tures of soft matter systems, e.g. the excited states of a chro- be used for further structural analysis, for example for calcula-
 3

tion of three-dimensional spatial or orientational correlations. where Sα is the number of pairs with interparticle distances
However, as the interaction potential is not recovered, the in- inside the α-slice. Evidently, Sα is an estimator of the radial
verse problem is not solved completely, and it is not possible distribution function: hSα i = 4πr2 g(rα )N 2 /(2V ). The av-
to calculate thermodynamical or dynamical properties using erage values of Sα are some functions of the potential Vα and
this approach. For the same reason this approach in not suit- can be written as an expansion
able for the development of coarse-grained models.
X ∂hSα i
In another approach, Reatto et al. [22] used the Hypernet- ∆hSα i = ∆Vγ + O(∆V 2 ). (5)
ted Chain (HNC) approximation to solve the inverse problem. γ
∂Vγ
This and similar approaches, based on some closure of the
Ornstein-Zernike equation, have been used in a number of The derivatives ∂hSα i/∂Vγ can be expressed using exact sta-
works during the last decade to compute interaction potentials tistical mechanics relationships [17]
from radial distribution functions [23, 24]. It should be noted, R P
however, that computations within the HNC theory are fea- ∂hSα i ∂ dq Sα (q) exp(−β λ Kλ Sλ (q))
= R P
sible only for relatively simple models. Moreover, although ∂Vγ ∂Vγ dq exp(−β λ Kλ Sλ (q))
yielding sometimes very accurate results [25], the HNC the- (6)
hSα Sγ i − hSα ihSγ i
ory is not an exact mathematical solution of the statistical- = − .
mechanical problem, and its accuracy should be investigated kB T
carefully in every specific case. Other works devoted to the Equations (5) and (6) allow us to find the interaction poten-
inverse problem are listed in Refs. [26, 27, 28]. tial Vα iteratively from the radial distribution functions hSα i.
A practical way to solve the inverse problem has been sug- (0)
Let Vα be a trial potential for which a natural choice is the
gested by Lyubartsev and Laaksonen [17]. This is the method
potential of mean force (PMF)
we will concentrate on in the following, and we will use it to
compute effective potentials from radial distribution functions Vα(0) = −kB T ln g(rα ), (7)
obtained from detailed MD simulations. The effective poten-
tials then allow us to to study the large-scale properties of a By carrying out standard Monte Carlo simulations, one can
given model system through coarse-grained DPD simulations. evaluate the averages hSα i and their deviations from the ref-
erence values Sα∗ defined from the given radial distribution
function as ∆hSα i(0) = hSα i(0) − Sα∗ . By solving the sys-
B. The Inverse Monte Carlo Method – A Simple Case tem of linear equations [Eq. (5)] with coefficients defined by
Eq. (6), and omitting terms O(∆V 2 ), we obtain corrections
To illustrate the Inverse Monte Carlo method, we consider (0)
to the potential ∆Vα . The procedure is then repeated with
the case of a single-component system consisting of identi- (1) (0) (0)
the new potential Vα = Vα + ∆Vα until convergence is
cal particles with pairwise interactions. The general case of
reached. The whole procedure resembles a solution of a mul-
a multicomponent system is a straightforward extension of it.
tidimensional non-linear equation using the Newton-Raphson
The Hamiltonian for this system is given as
method.
X If the initial approximation of the potential is poor, some
H= V (rij ), (1) regularization of the iteration procedure is needed. In this
i,j
case we multiply the required change of the radial distribu-
where V (rij ) is the pair potential and rij is the distance be- tion function by a small factor that is typically between 0 and
tween the interaction sites i and j. Let us assume that we 1. By doing so, the term O(∆V 2 ) in Eq. (5) can be made
know the radial distribution function g(rij ). Our aim is now small enough to guarantee convergence of the whole proce-
to find the corresponding interaction potential V (rij ). dure, although the number of iterations will increase.
Let us introduce the following grid approximation to the The above algorithm has also the advantage that it provides
Hamiltonian, us with a method to evaluate the uncertainty of the inverse
procedure. The radial distribution function has normally some
Ve (r) = V (rα ) ≡ Vα (2) uncertainty. An analysis of the eigenvalues and eigenvectors
of the matrix in Eq. (6) allows one to make conclusions of
for which changes in g(r) correspond to which changes in the
potential. For example, eigenvectors with eigenvalues close
1 1
rα − < r < rα + and rα = (α−0.5) rcut /M, to zero correspond to changes in the potential which have al-
2M 2M most negligible effect on the radial distribution function. The
(3)
where α = 1, . . . , M , and M is the number of grid points presence of these small eigenvalues makes the inverse prob-
within the interval [0, rcut ], and rcut is a chosen cut-off dis- lem not well-defined, however. In some cases, such as for
tance. Then, the Hamiltonian in Eq. (1) can be rewritten as liquid water, that may pose serious problems in the inversion
procedure [29].
It is instructive to note a relation between the present ap-
X
H= Vα Sα , (4) proach and the renormalization group Monte Carlo method
α [30, 31] used to study phase transitions in lattice models (e.g.,
 4

polymer models and Ising models for ferromagnets) as well quate. It was shown by Forrest and Suter in 1995 that coarse-
as in the quantum field theory. The renormalization procedure graining of molecular representation leads to this type of soft-
introduces a scale change in the system. During this procedure ening [33]. However, in many cases this formulation is of too
one consecutively obtains a more and more coarse-grained de- generic nature. Although it is possible to use mean field theo-
scription of the system. Equations (5) and (6) were in fact ries such as the Flory-Huggins theory for polymers to find the
used by Swendsen et al. [30, 31] to describe how the pa- amplitudes for the conservative forces between different types
rameters of a Hamiltonian change during the renormalization of particles [32], there are cases where this approach does not
procedure. It now appears that the applications of this method retain enough information about the actual character of the
are more general than the original lattice systems near a phase different atoms, molecules, or interactions (see discussion be-
transition, and cover even soft matter problems allowing us to low).
“renormalize” Hamiltonians of molecular systems in such a The weight functions for the dissipative and random forces,
way that only the degrees of freedom of primary interest are ω D and ω R , cannot be chosen arbitrarily. This is easy to
kept in the coarse-grained system. understand by the intuitive argument that the thermal heat
generated by the random force must be balanced by dissipa-
tion. Español and Warren [3] studied this situation analyti-
C. Dissipative Particle Dynamics (DPD) cally and derived a fluctuation-dissipation relation connecting
the weight functions and amplitudes of the dissipative and ran-
Dissipative particle dynamics was introduced in 1992 by dom forces via
Hoogerbrugge and Koelman for simulations of hydrodynamic ω D (rij ) = [ω R (rij )]2 and σ 2 = 2 γ kB T. (11)
phenomena in complex fluids [1, 2]. The original formulation
did not obey detailed balance, though, and in 1995 Español These conditions guarantee convergence to the canonical en-
and Warren [3] formulated a new DPD algorithm which they semble as required. It is important to notice that ω D and ω R
showed to be fully consistent with statistical mechanics. This are completely decoupled from the conservative part. In other
algorithm is nowadays known as DPD. In the following we words, we can consider the DPD algorithm as a momentum
present the basic DPD formalism and discuss some of its fea- conserving thermostat for an arbitrary conservative potential.
tures that are relevant to the present work. This is in fact the way how we apply the DPD algorithm in
In the basic formulation of DPD, the interactions are pair- our coarse-grained simulations: The conservative potential is
wise additive and the force exerted by particle j on particle defined through the Inverse Monte Carlo procedure and then
i is given as a sum of conservative, dissipative, and random used in the DPD algorithm.
forces through F~ij = F~ijC + F~ijD + F~ijR , respectively. These Finally, to study the dynamics of the system one needs to
forces are typically given as evolve the system in time. In DPD this is simply done by inte-
grating the Newton’s equations of motion. In our simulations
(c)
F~ijC = Fij (rij ) ~eij , (8) we used a velocity-Verlet based algorithm adapted for DPD
(sometimes called DPD–VV) [19, 20].
F~ijD = −γ ω D (rij )(~vij · ~eij ) ~eij , (9)
F~ijR = σ ω R (rij ) ξij ~eij , (10) D. Combining MD and DPD
where ~rij = ~ri − ~rj is the position vector connecting two par-
ticles, rij = |~rij | is the interparticle distance, ~eij = ~rij /rij With the above definitions we now have all the tools for
is the corresponding unit vector, and ~vij = ~vi − ~vj is the coarse-graining. We remove fast internal and orientational de-
relative velocity of particles i and j. The terms ξij are sym- grees of freedom from all water molecules, and represent wa-
metric random variables with zero mean and unit variance, ter molecules as one-site particles interacting with other wa-
and are independent for different pairs of particles and dif- ter particles and ions by spherically symmetric potentials. In
ferent times. The conservative force is essentially given by comparison to the original all-atom model, we have approxi-
(c)
Fij (rij ). Constants γ and σ give the amplitudes of the dis- mately one third of the degrees of freedom left. This can be
interpreted as intermediate level of coarse-graining.
sipative and random forces, respectively, and ω D and ω R are Without any loss of generality we can choose the weight
the corresponding weight functions. functions to be of the standard form ω R (rij ) = 1 − rij /rc
Let us now return to the different forces. The DPD for- as discussed above. This choice has been made invariably
mulation does not specify the form of the conservative force. in DPD simulations. Since the Inverse Monte Carlo proce-
The most common choice in standard DPD simulations is dure provides us with the effective potentials, thus yielding us
(c) (c)
Fij (rij ) = αij (1 − rij )/rc for r ≤ rc and Fij (rij ) = 0
F~ijC (r), the only task left is to ensure that the MD and DPD
otherwise. In other words, the conservative force is derived
systems correspond to each other. In this study, that was done
from a soft pair potential U = αij (1 − rij /rc )2 /2, where rc
by adjusting the amplitude of the dissipative force γ which de-
is the cutoff distance. Interactions between different types of
termines the dissipation rate in the DPD system. Here, γ was
coarse-grained particles can then be described by varying the
determined by requiring the short-time decay of the single-
amplitude of the conservative interaction αij [32].
particle velocity autocorrelation function
The above is the preferred formulation for the conserva-
tive force when generic, simple, and soft interactions are ade- φ(t) ≡ h~vi (t) · ~vi (0)i (12)
 5

FIG. 1: The decay of the single-particle velocity autocorrelation function φ(t) at early times. Results shown here are for water, Na+ and Cl−
ions. As the data illustrates, the early-time decay of φ(t) is essentially identical between MD and DPD simulations for γ = 0.72 used in the
DPD simulations.

to be approximately the same in both MD and DPD systems may be constructed by the Inverse Monte Carlo method from
as shown in Fig. 1. It shows how the early-time decay of φ(t) ion-ion radial distribution functions generated in high-quality
can be matched, leading to a value of γ = 0.72 for the DPD all-atomic molecular dynamics simulations [17, 25, 34]. This
model. The long-time behavior of the velocity autocorrelation is the approach we used in obtaining the effective potentials.
function between MD and DPD is different, however. This is The MD simulations were performed in the NVT ensem-
an obvious result since some microscopic degrees of freedom ble using the flexible SPC water model [35] and Smith–Dang
have been coarse grained and thus effects due to hydrogen parameters for Na+ and Cl− ions [36], i.e., the Lennard-
bonds, for example, are implicitly included in the effective in- Jones parameters for the sodium ions were σ = 2.35 Å and
teractions. This leads to an enhanced diffusion rate at inter- ε = 0.544 kJ/mol, and for chloride σ = 4.4 Å and ε =
mediate times, and is demonstrated in Fig. 1 as a positive tail 0.419 kJ/mol. In the results reported here we set the temper-
for the DPD particles. ature to 300 K. The electrostatic interactions were taken into
This approach is meaningful since (i ) it makes sure that account by Ewald summation. Other simulation details are
the early-time dynamics is described properly, while (ii ) it described in Refs. [17, 34]. From the MD simulations, the ra-
does not fix the tracer diffusion coefficient. One should no- dial distribution functions between different pairs of particles
tice that the tracer diffusion coefficient is an integral over the were calculated and fed as an input in the Inverse Monte Carlo
velocity correlation function over long times, and in hydro- procedure. The results are shown in Figs. 2 and 3.
dynamic systems the long-time tail is important. With this Figure 2 shows effective potentials between water
choice, we can assume the diffusion coefficient to be an inde- molecules, presented as one-center particles, and between
pendent quantity so that its behavior, found by DPD, can be other water molecules and ions. For comparison, water-water
compared with both MD simulations and experiments. Fixing effective potential calculated from MD simulations for pure
γ thus determines the transport properties of the system. water is also displayed. It is interesting to note that the pres-
ence of ions has practically no effect on the water-water effec-
tive potential.
III. RESULTS Figure 3 displays the effective potentials between the ions.
They are compared to the effective potentials calculated with-
A. Constructing the potentials out water [17, 34]. The effective potentials are very similar
in the two cases both of them having an oscillating charac-
One of the typical simplifications used in molecular simu- ter with one or two oscillations. The potential approaches the
lations is the replacement of solvent molecules by continuum primitive model potential with dielectric constant close to 80.
media. For example, in the primitive electrolyte model, ions At distances above 10 Å the effective potential coincides al-
in water are substituted by charged spheres moving in dielec- most perfectly with the Coulombic potential.
tric media with the dielectric constant set to about 80. This is Within the inverse Monte Carlo procedure, calculation of
a serious simplification at small distances (a few ångströms) effective potentials without water is much easier than in the
where it is impossible to define a dielectric constant. Besides presence of water. On the other hand, the presence of particles
this, in the primitive electrolyte model the ion radius is an ad- representing water is necessary in DPD simulations. Our test
justable parameter without any clear physical meaning. A bet- studies have shown that while the use of effective potentials
ter model for effective ion-ion interactions in an aqueous solu- calculated without water may give qualitatively satisfactory
tion must take into account the solvation structure around the results, inclusion of water in the Inverse Monte Carlo simula-
ions. Practically, effective solvent-mediated ion-ion potentials tions essentially leads to an improvement of results, especially
 6

B. Coarse-grained simulations

In our first paper using this approach [18], we showed that

the radial distribution functions computed from the coarse-
grained DPD simulations for the pairs between the Na+ and
Cl− ions agreed with those obtained from the all-atom MD
simulations. In addition to the static properties, the diffusion
coefficients were shown to be in good agreement between the
MD and DPD results at various NaCl concentrations between
0.5 M and 4.1 M. What makes those results very interesting is
the fact that, for all the DPD simulations, we used the effec-
tive potentials obtained at one single NaCl concentration only.
In that case, the effective potentials were obtained from MD
simulations at a salt concentration of 0.87 M and then used
in DPD simulations for various salt concentrations between
0.5 M and 4.1 M. When compared with MD simulations, the
agreement was found to be excellent. Here, we present results
FIG. 2: Effective potentials for coarse-grained water in NaCl solu- for static correlations, coordination numbers, velocity auto-
tion. Thin dotted line is water-water effective potential in pure water. correlation functions, and viscosity.
“O” refers to an oxygen atom in the water molecule.

molarity: 0.55 M 0.87 M 2.2 M 4.1 M

MD: Na+ : 5.8 5.5 5.42 5.3
Cl− : 6.9 6.8 6.9 7.1
DPD: Na+ : 5.67 5.66 5.55 5.40
Cl− : 6.91 6.96 6.97 7.00

TABLE I: Coordination numbers obtained from MD and DPD sim-

ulations at various salt concentrations. The error bars for both cases
are ±0.05. The MD data is from Ref. [39].

It is important to notice that due to the self-consistency

conditions for the effective potentials, the DPD simulations
will always produce the correct pair correlation behavior at
the point in phase space where the effective interactions were
determined in the first place. As an example, this is demon-
strated in Fig. 4 for all different kinds of pairs of particles stud-
FIG. 3: Ion-ion effective potentials in aqueous solution. Solid lines
describe potentials in presence of one-site water molecules (this
ied in this work. It is very relevant to ask, however, how the
work), while dashed lines are for effective potential without water effective interactions change when system variables such as
from Ref. [34]. the salt concentration or the temperature of the system change.
We have explored this matter through combined MD and DPD
simulations and found that for a wide range of salt concentra-
tions in the present system it is not necessary to recalculate
in the case of dynamical properties. the effective potentials. This is demonstrated by the data in
It is worth pointing out that the effective potential dif- Table I which compares the coordination numbers obtained
fers from the potential of mean force defined as ψP MF = from DPD and MD simulations for g(r). In the DPD sim-
−kB T ln g(r), which corresponds to the Kirkwood approx- ulations we have again used the effective potentials obtained
imation for the n-particle correlation function. The potential at 0.87 M. We find that both MD and DPD follow the same
of mean force is screened by the other ions in a system and qualitative behavior. For sodium the coordination number de-
decays as (1/r)e(−r/rD ) , i.e., not as the 1/r Coulombic po- creases with increasing salt concentration whereas for chlo-
tential. The potential of mean force is therefore not a very ride we find an opposite trend. Furthermore, the quantitative
suitable choice to present ion-ion interactions within a con- agreement between the MD and DPD results is remarkably
tinuum solvent model at finite ion concentrations even in a good.
coarse-grained solvent system. The effective potentials in- Our finding that the effective interactions depend only
clude a contribution from Coulombic forces, however. weakly on thermodynamic conditions is also supported by
 7

FIG. 4: Radial distribution functions g(r) for different pairs of particles in the aqueous NaCl system studied by both DPD (shown on the left)
and MD (on the right). The studies were made at a salt concentration of 0.87 M. “O” refers to the ogygen atom in the water molecule. Note
the agreement between MD and DPD results.

the present results for shear viscosity (see below) and ear- However, since the shear viscosity coefficient is a collective
lier studies of the same system for tracer diffusion [18]. In quantity, meaning that all particles in a system give rise to a
Ref. [18], we found that DPD simulations with effective po- single sample, an accurate determination of the shear viscos-
tentials yielded tracer diffusion coefficients in good agreement ity coefficient through MD simulations is both difficult and
with MD simulations. Further support for this result is given very time-consuming. Thus it was not done in the present
by previous work [25, 34], where the dependence of effective case, and as a matter of fact, to the best of our knowledge,
potentials on the salt concentration was studied in detail. It there are no published reports of MD simulations of shear vis-
was found that effective potentials depend very little on the cosity in NaCl solutions. Thus, we compare our DPD data
salt concentration: An increase in salt concentration leads to to experimental results [37] instead. The results shown in
a slow decrease in the dielectric permittivity. Fig. 5 indicate that the qualitative behavior of the shear viscos-
Our results thus suggest that the effective interactions ity coefficient obtained through DPD simulations is in good
change only slightly when the system parameters are varied, agreement with the experiments. This result is truly remark-
as is here the case for salt concentration. However, while this able as it shows that both the equilibrium (static) and dynamic
result holds true in the present system for an aqueous solution behavior of the system are properly described by the coarse-
of monovalent ions, the generality of this finding remains an grained approach. This provides us with strong support that
open question. Further studies of divalent systems and more the present approach of coarse graining an MD system and ap-
complex molecules, among others, are therefore called for. plying the obtained effective interactions in DPD simulations
Next, we determined the shear viscosity using a Green- is a promising approach indeed.
Kubo relation. As with the other quantities, the shear viscos- We would like to note that we have very recently performed
ity should be compared to results from the MD simulations. simulations using LiCl and CaCl2 [38] to study further aspects
 8

2.5 this work, the computational gain is still remarkable, around

a factor of 20. This speed-up is due to the fact that no explicit
Ewald summation was required for the coarse-grained system
2.0 (since the effective potentials include information about the
[mPa s]

electrostatics, and the screening length was shorter than the

applied cut-off) and that SPC water used in molecular dynam-
1.5 ics simulations is replaced by a single coarse-grained particle
with the electrostatic interactions included in the effective po-
1.0 tentials.

0.5 This approach allows one to coarse-grain complex macro-

DPD molecular systems to a level where a detailed molecular de-
Experiments scription is replaced with a simplified model of coarse-grained
0.0 particles. For example, a lipid-acyl chain “. . . – CH2 – CH2 –
CH2 – CH2 – CH2 – CH2 . . .” can be described as “. . .– B –
0 1 2 3 4 5 B – . . .” where each of the beads “B” represent 2 – 4 methyl
groups CH2 . Since the beads “B” describe the interaction site
Salt concentration [M] for a whole molecular group, it is clear that the interactions
between different beads will be relatively soft in the spirit of
the original DPD formulation, yet retaining their characteristic
FIG. 5: The shear viscosity coefficient η from coarse-grained features well beyond the mean-field description of the original
DPD simulations and experiments. The experimental data is from
DPD formulation. For macromolecular systems, this coarse-
Ref. [37]. The data by DPD has been scaled to allow comparison.
graining procedure thus has the potential to lead to a major
reduction in computational effort as well as to softer interac-
tion potentials. The computational gain depends on the type
of this issue, such as the effects of salt concentration, valence,
of interactions: While here in the case of NaCl we did not
and temperature. Despite the very different natures of these
need to use Ewald summation, the situation may be different
solutions (as compared to NaCl), the conclusions with regard
in the case of divalent ions such as calcium. In addition, one
to the coarse-graining method and its general applicability re-
should note that coarse-graining of a macromolecular system
main the same.
requires a lot of insight, and there is no unique way of select-
ing the units to be coarse-grained. Work is in progress to ad-
dress these issues and the coarse-graining of macromolecular
IV. CONCLUSIONS
systems in general.
In this article we have described in detail a method that
combines the Inverse Monte Carlo method and the momen- The approach presented in this article has close connections
tum conserving dissipative particle dynamics thermostat for to other methods as well. The Inverse Monte Carlo procedure
coarse-grained simulations. This approach allows systematic has the same foundations as the integral equation methods of
coarse-graining of molecular systems and conserves hydro- Bolhuis et al. [10, 11]. Second, the coarse-graining procedure
dynamics. Here, the level of coarse-graining has been quite used by Reith, Mayer and Müller-Plathe [15] has a close re-
modest and the potentials used did not have particularly soft semblance to the Inverse Monte Carlo procedure. The latter
cores. This is due to the fact that in this system of water and group has also applied their method to coarse-grain polymeric
NaCl, we coarse-grained the molecular representation of wa- systems in the spirit described in the previous paragraph. The
ter molecules only. Thus, the corresponding effective interac- novelty here is the application of the effective potentials in
tions retained features which are typical for systems described coarse-grained simulations by using the dissipative particle
on an atomic level. However, it is important to notice two dynamics thermostat. The results for the present system in-
facts. First of all, further coarse-graining is not possible for dicate that dynamic properties could be studied using this ap-
this kind of a system if one wishes to preserve the identity proach. Whether or not this holds true in more complex sys-
of individual ions and the essential characteristics of water tems, and what are the limits of this approach are topics of
molecules. One could, of course, continue the coarse-graining further studies.
procedure to obtain softer potentials by clustering several wa-
ter molecules together and treating the clusters as new “effec-
tive” solvent particles. This would lead to further computa- In conclusion, the agreement between the DPD simula-
tional gains but the ions and molecules would then lose essen- tions, the MD results and experiments strongly indicates that
tial parts of their identity. As a result of this the interactions the method presented in this work is both physically sound
would become softer and resemble the quadratic potentials and provides a controllable method for performing coarse-
used for conservative interactions in the original DPD formu- graining and mesoscale simulations at least to a moderate de-
lation. Second, despite the modest level of coarse-graining in gree.
 9

Acknowledgments (M.K.). Further support has been obtained from the Academy
of Finland though its Centre of Excellence Program (I.V.) and
This work has been supported by the Swedish Science from the Finnish Academy of Science and Letters (I.V.).
Council (A.P.L. and A.L.) and the Academy of Finland

[1] Hoogerbrugge, P. J.; Koelman, J. M. V. A. Europhys. Lett. 1992, Rev. E 2000, 62, R7611–R7614.
19, 155–160. [20] Vattulainen, I.; Karttunen, M.; Besold, G.; Polson, J. M. J.
[2] Koelman, J. M. V. A.; Hoogerbrugge, P. J. Europhys. Lett. 1993, Chem. Phys. 2002, 116, 3967–3979.
21, 363–368. [21] McGreevy, R. L.; Pusztai, L. Mol. Sim. 1988, 1, 359–367.
[3] Español, P.; Warren, P. Europhys. Lett 1995, 30, 191–196. [22] Reatto, L.; Levesque, D.; Weis, J. J. Phys. Rev. A 1986, 33,
[4] Warren, P. B. Curr. Opin. Coll. Interf. Sci. 1998, 3, 620–624. 3451–3465.
[5] Flekkøy, E. G.; Coveney, P. V. Phys. Rev. Lett. 1999, 85, 2522. [23] Rosenfeld, Y.; Kahl, G. J. Phys.: Cond. Mat. 1997, 9, L89–L98.
[6] Flekkøy, E. G.; Wagner, G.; Feder, J. Europhys. Lett. 2000, 52, [24] Gonzalez-Mozuelos, P.; Carbajal-Tinoco, M. D. J. Chem. Phys.
271–276. 1998, 24, 11074–11084.
[7] Malevanets, A; Kapral, R. J. Chem. Phys. 1999, 110, 8605– [25] Lyubartsev, A. P.; Marcelja, S. Phys. Rev. E 2002, 65, 041202.
8613. [26] Ostheimer, M.; Bertagnolli, H. Mol. Sim. 1989, 3, 227–233.
[8] Malevanets, A; Kapral, R. J. Chem. Phys. 2000, 112, 7260– [27] Soper, A. K. Chem. Phys. 1996, 202, 295–306.
7269. [28] Toth, G.; Baranyai, A. J. Mol. Liquids 2000, 85, 3–9.
[9] Akkermans, R. L. C; Briels, W. J. J. Chem. Phys. 2000, 113, [29] Lyubartsev, A. P.; Laaksonen, A. Chem. Phys. Lett. 2000, 325,
6409–6422. 15–21.
[10] Louis, A. A.; Bolhuis, P. G.; Hansen, J. P.; Meijer, E. P. Phys. [30] Swendsen, R. H. Phys. Rev. Lett. 1979, 42, 859–861.
Rev. Lett. 2000, 85, 2522. [31] Pawley, G. S.; Swendsen, R. H.; Wallace, D. J.; Wilson, K. G.
[11] Bolhuis, P. G.; Louis, A. A.; Hansen, J. P.; Meijer, E. P. J. Chem. Phys. Rev. B 1984, 29, 4030–4040.
Phys. 2001, 114, 4296–4311. [32] Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423–
[12] Hansen, J. P.; McDonald, I. R. Theory of simple liquids; Aca- 4435.
demic press: London, 1986. [33] Forrest, B. M.; Suter, U. W. J. Chem. Phys. 1995, 102, 7256–
[13] Henderson, R. L. Phys. Lett. 1974, 49A, 197–198. 7266.
[14] Murat, M; Kremer, K. J. Chem. Phys. 1998, 108, 4340–4348. [34] Lyubartsev, A. P.; Laaksonen, A. Phys. Rev. E 1997, 55, 5689–
[15] Reith, D; Meyer, H.; Müller-Plathe, F. Macromolecules 2001, 5696.
34, 2235–2245. [35] Toukan, K.; Rahman, A. Phys. Rev. B 1985, B31, 2643–2648.
[16] Baschnagel, J.; Binder, K.; Doruker, P.; Gusev, A.; Hahn, O.; [36] Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757–
Kremer, K.; Mattice, W. L.; Müller-Plathe, F.; Murat, M.; Paul, 3766.
W.; Santos, S.; Suter, U. W.; Tries, V. Adv. Polym. Sci 2000, [37] Lide, D. R., Ed. CRC Handbook of Chemistry and Physics;
152, 41–156. CRC Press: Boca Raton, 82nd ed., 2001.
[17] Lyubartsev, A. P.; Laaksonen, A. Phys. Rev. E 1995, 52, 3730– [38] Terämä, E.; Vattulainen, I.; Patra, M.; Karttunen, M.; Lyubart-
3737. sev, A. P.; Laaksonen, A. in preparation 2002.
[18] Karttunen, M.; Laaksonen, A.; Lyubartsev, A. P.; Vattulainen, [39] Lyubartsev, A. P.; Laaksonen, A. J. Chem. Phys. 1996, 100,
I. submitted 2002. 16410–16418.
[19] Besold, G.; Vattulainen, I.; Karttunen, M.; Polson, J. M. Phys.