scispace - formally typeset
Search or ask a question
Journal ArticleDOI

Efficient Algorithms for Langevin and DPD Dynamics.

13 Jun 2012-Journal of Chemical Theory and Computation (AMER CHEMICAL SOC)-Vol. 8, Iss: 10, pp 3637-3649
TL;DR: Several algorithms for stochastic dynamics, including Langevin dynamics and different variants of Dissipative Particle Dynamics (DPD), applicable to systems with or without constraints are presented, showing that the measured thermal relaxation rates agree well with theoretical predictions.
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.

Content maybe subject to copyright    Report

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

Ecient 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 dierent
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 diusion 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 inuence of various parameters on the diusion coecient is discussed.
1. INTRODUCTION
The purpose of this article is to present simple and e 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 ow component. Galilean invariance is essential for
simulations that are required to follow NavierStokes 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.
14
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 rst-order decay of temperature
deviations toward the reference temperature. Global thermo-
stats of the weak-coupling type
5
show a convenient rst-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 dierent 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 rst-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 ying ice-cube eect. Global
thermostats of the extended-system type, such as the Nose
-Hoover
thermostat,
9
will maintain a canonical distribution in congura-
tional space but show strong oscillatory behavior of the
temperature deviation as a result of the order of the dierential
equation related to the systems extension. The Nose
-Hoover
chain thermostat,
10
using a cascade of thermostats, alleviates but
does not solve this problem. See Berendsen,
11
pp 194204, 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, 36373649

constant between friction force and velocity, which equals the
mass of the particle times the friction rate.]
This equationand its DPD equivalent,
1214
which applies
friction and noise to velocity differences between particle pairs
is a stochastic dierential equation. We note that eq 2 is a
formal equation that is mathematically incorrect because v(t)is
not dierentiable; it should be implemented as a dierence
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 nite 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
Adierent 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 modies 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 dened
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 eective thermostat. The energy and
momentum transfer implicit in the impulsive friction will
inuence 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 dierence wit h the usual stochastic dierential
equation is that the time evolution of the system is not
described by a single stochastic dierential 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 dierences. 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 eective thermostat
strength implied by the impulsive application of friction and
noise. We also investigate whether the pairwise application can
be made more ecient 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 diusion coecient,
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 dened that serves
as a comparable friction rate measure, with a comparable
inuence on diusion, for any type of impulsive friction
imposed on the system. Section 6 describes the simulation
details and reports the computational eciency of the dierent
methods; section 7 gives the results of various tests of the
algorithms on dierent molecular systems (ideal gas, coarse-
grained water without constraints, and atomic water with
constraints). Both t hermostat and diusion 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)(/)
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, 363736493638

+=+Δ
⎜⎟
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 modied 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., dierent f s for dierent 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 dier for dierent
particles. The application to velocity dierences 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,wend 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, 363736493639

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 dierential 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 specic heat per particle.
Equation 15 shows that any deviation from the reference
temperature will decay to zero according to a rst-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 uids, this fraction is
usually 2 to 3 times smaller.
The thermal rate constant can be expressed in an effective
friction rate γ
e
,dened 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 fullled. 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 rst 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 satised by
calling a routine constr (SHAKE,
20
LINCS,
21
SETTLE
22
) that
makes corrections to an unconstrained conguration x
u
to yield
a conguration x
c
that satises the set of m holonomic
constraints σ
k
(x)=0,k = 1, ..., m. The displacements are in the
direction of a reference conguration 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 modied 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 congurational space, dened 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 dierences 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, 363736493640

Citations
More filters
Journal ArticleDOI
TL;DR: The Martini model, a coarse-grained force field for biomolecular simulations, has found a broad range of applications since its release a decade ago and is described as a building block principle model that combines speed and versatility while maintaining chemical specificity.
Abstract: The Martini model, a coarse-grained force field for biomolecular simulations, has found a broad range of applications since its release a decade ago. Based on a building block principle, the model combines speed and versatility while maintaining chemical specificity. Here we review the current state of the model. We describe recent highlights as well as shortcomings, and our ideas on the further development of the model.

1,022 citations

Journal ArticleDOI
TL;DR: QuantumATK as discussed by the authors is an integrated set of atomic-scale modelling tools developed since 2003 by professional software engineers in collaboration with academic researchers, which enable electronic-structure calculations using density functional theory or tight-binding model Hamiltonians, and also offers bonded or reactive empirical force fields in many different parametrizations.
Abstract: QuantumATK is an integrated set of atomic-scale modelling tools developed since 2003 by professional software engineers in collaboration with academic researchers. While different aspects and individual modules of the platform have been previously presented,a#13; the purpose of this paper is to give a general overview of the platform. The QuantumATK simulation engines enable electronic-structure calculations using density functional theory or tight-binding model Hamiltonians, and also offers bonded or reactive empirical force fields in many different parametrizations. Density functional theory is implemented using either a plane-wave basis or expansion of electronic states in a linear combination of atomic orbitals. The platform includes a long list of advanced modules, including Green's-function methods for electron transport simulations and surface calculations, first-principles electron-phonon and electron-photon couplings,a#13; simulation of atomic-scale heat transport, ion dynamics, spintronics, optical properties of materials, static polarization, and more.a#13; Seamless integration of the different simulation engines into a common platform allows for easy combination of different simulation methods into complex workflows. Besides giving a general overview and presenting a number of implementation detailsa#13; not previously published, we also present four different application examples. These are calculations of the phonon-limited mobility of Cu, Ag and Au, electron transport in a gated 2D device, multi-model simulation of lithium ion drift through a battery cathode in an external electric field, and electronic-structure calculations of the composition-dependent band gap of SiGe alloys.a#13;

658 citations

Journal ArticleDOI
TL;DR: This study investigates how six well-established thermostat algorithms applied with different coupling strengths and to different degrees of freedom affect the dynamics of various molecular systems.
Abstract: Temperature control algorithms in molecular dynamics (MD) simulations are necessary to study isothermal systems. However, these thermostatting algorithms alter the velocities of the particles and thus modify the dynamics of the system with respect to the microcanonical ensemble, which could potentially lead to thermostat-dependent dynamical artifacts. In this study, we investigate how six well-established thermostat algorithms applied with different coupling strengths and to different degrees of freedom affect the dynamics of various molecular systems. We consider dynamic processes occurring on different times scales by measuring translational and rotational self-diffusion as well as the shear viscosity of water, diffusion of a small molecule solvated in water, and diffusion and the dynamic structure factor of a polymer chain in water. All of these properties are significantly dampened by thermostat algorithms which randomize particle velocities, such as the Andersen thermostat and Langevin dynamics, when...

283 citations

Journal ArticleDOI
TL;DR: An implicit-solvent version of the popular CG Martini model, nicknamed "Dry" Martini, is introduced, to account for the omitted solvent degrees of freedom, and the nonbonded interaction matrix underlying the Martini force field was reparametrized.
Abstract: Coarse-grained (CG) models allow simulation of larger systems for longer times by decreasing the number of degrees of freedom compared with all-atom models. Here we introduce an implicit-solvent version of the popular CG Martini model, nicknamed "Dry" Martini. To account for the omitted solvent degrees of freedom, the nonbonded interaction matrix underlying the Martini force field was reparametrized. The Dry Martini force field reproduces relatively well a variety of lipid membrane properties such as area per lipid, bilayer thickness, bending modulus, and coexistence of liquid-ordered and disordered domains. Furthermore, we show that the new model can be applied to study membrane fusion and tether formation, with results similar to those of the standard Martini model. Membrane proteins can also be included, but less quantitative results are obtained. The absence of water in Dry Martini leads to a significant speedup for large systems, opening the way to the study of complex multicomponent membranes containing millions of lipids.

240 citations

Journal ArticleDOI
TL;DR: Free energy calculations based on molecular dynamics and thermodynamic cycles accurately reproduce experimental affinities of diverse bromodomain inhibitors.
Abstract: Accurate prediction of binding affinities has been a central goal of computational chemistry for decades, yet remains elusive. Despite good progress, the required accuracy for use in a drug-discovery context has not been consistently achieved for drug-like molecules. Here, we perform absolute free energy calculations based on a thermodynamic cycle for a set of diverse inhibitors binding to bromodomain-containing protein 4 (BRD4) and demonstrate that a mean absolute error of 0.6 kcal mol−1 can be achieved. We also show a similar level of accuracy (1.0 kcal mol−1) can be achieved in pseudo prospective approach. Bromodomains are epigenetic mark readers that recognize acetylation motifs and regulate gene transcription, and are currently being investigated as therapeutic targets for cancer and inflammation. The unprecedented accuracy offers the exciting prospect that the binding free energy of drug-like compounds can be predicted for pharmacologically relevant targets.

238 citations

References
More filters
Journal ArticleDOI
TL;DR: In this paper, a method is described to realize coupling to an external bath with constant temperature or pressure with adjustable time constants for the coupling, which can be easily extendable to other variables and to gradients, and can be applied also to polyatomic molecules involving internal constraints.
Abstract: In molecular dynamics (MD) simulations the need often arises to maintain such parameters as temperature or pressure rather than energy and volume, or to impose gradients for studying transport properties in nonequilibrium MD A method is described to realize coupling to an external bath with constant temperature or pressure with adjustable time constants for the coupling The method is easily extendable to other variables and to gradients, and can be applied also to polyatomic molecules involving internal constraints The influence of coupling time constants on dynamical variables is evaluated A leap‐frog algorithm is presented for the general case involving constraints with coupling to both a constant temperature and a constant pressure bath

25,256 citations

Journal ArticleDOI
TL;DR: In this paper, a numerical algorithm integrating the 3N Cartesian equations of motion of a system of N points subject to holonomic constraints is formulated, and the relations of constraint remain perfectly fulfilled at each step of the trajectory despite the approximate character of numerical integration.

18,394 citations

Journal ArticleDOI
TL;DR: The dynamical steady-state probability density is found in an extended phase space with variables x, p/sub x/, V, epsilon-dot, and zeta, where the x are reduced distances and the two variables epsilus-dot andZeta act as thermodynamic friction coefficients.
Abstract: Nos\'e has modified Newtonian dynamics so as to reproduce both the canonical and the isothermal-isobaric probability densities in the phase space of an N-body system. He did this by scaling time (with s) and distance (with ${V}^{1/D}$ in D dimensions) through Lagrangian equations of motion. The dynamical equations describe the evolution of these two scaling variables and their two conjugate momenta ${p}_{s}$ and ${p}_{v}$. Here we develop a slightly different set of equations, free of time scaling. We find the dynamical steady-state probability density in an extended phase space with variables x, ${p}_{x}$, V, \ensuremath{\epsilon}\ifmmode \dot{}\else \.{}\fi{}, and \ensuremath{\zeta}, where the x are reduced distances and the two variables \ensuremath{\epsilon}\ifmmode \dot{}\else \.{}\fi{} and \ensuremath{\zeta} act as thermodynamic friction coefficients. We find that these friction coefficients have Gaussian distributions. From the distributions the extent of small-system non-Newtonian behavior can be estimated. We illustrate the dynamical equations by considering their application to the simplest possible case, a one-dimensional classical harmonic oscillator.

17,939 citations

Journal ArticleDOI
TL;DR: A new implementation of the molecular simulation toolkit GROMACS is presented which now both achieves extremely high performance on single processors from algorithmic optimizations and hand-coded routines and simultaneously scales very well on parallel machines.
Abstract: Molecular simulation is an extremely useful, but computationally very expensive tool for studies of chemical and biomolecular systems Here, we present a new implementation of our molecular simulation toolkit GROMACS which now both achieves extremely high performance on single processors from algorithmic optimizations and hand-coded routines and simultaneously scales very well on parallel machines The code encompasses a minimal-communication domain decomposition algorithm, full dynamic load balancing, a state-of-the-art parallel constraint solver, and efficient virtual site algorithms that allow removal of hydrogen atom degrees of freedom to enable integration time steps up to 5 fs for atomistic simulations also in parallel To improve the scaling properties of the common particle mesh Ewald electrostatics algorithms, we have in addition used a Multiple-Program, Multiple-Data approach, with separate node domains responsible for direct and reciprocal space interactions Not only does this combination of a

14,032 citations

Journal ArticleDOI
TL;DR: Although the derivation of the algorithm is presented in terms of matrices, no matrix matrix multiplications are needed and only the nonzero matrix elements have to be stored, making the method useful for very large molecules.
Abstract: In this article, we present a new LINear Constraint Solver (LINCS) for molecular simulations with bond constraints. The algorithm is inherently stable, as the constraints themselves are reset instead of derivatives of the constraints, thereby eliminating drift. Although the derivation of the algorithm is presented in terms of matrices, no matrix matrix multiplications are needed and only the nonzero matrix elements have to be stored, making the method useful for very large molecules. At the same accuracy, the LINCS algorithm is three to four times faster than the SHAKE algorithm. Parallelization of the algorithm is straightforward. (C) 1997 John Wiley & Sons, Inc.

12,699 citations