Eﬀects of Temperature Control Algorithms on Transport Properties
and Kinetics in Molecular Dynamics Simulations
Joseph E. Basconi and Michael R. Shirts*
Department of Chemical Engineering, University of Virginia, Charlottesville, Virginia 22904-4741, United States
*S Supporting Information
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 diﬀerent coupling strengths
and to diﬀerent degrees of freedom aﬀect the dynamics of
various molecular systems. We consider dynamic processes
occurring on diﬀerent times scales by measuring translational and rotational self-diﬀusion as well as the shear viscosity of water,
diﬀusion of a small molecule solvated in water, and diﬀusion and the dynamic structure factor of a polymer chain in water. All of
these properties are signiﬁcantly dampened by thermostat algorithms which randomize particle velocities, such as the Andersen
thermostat and Langevin dynamics, when strong coupling is used. For the solvated small molecule and polymer, these
dampening eﬀects are reduced somewhat if the thermostats are applied to the solvent alone, such that the solute’s temperature is
maintained only through thermal contact with solvent particles. Algorithms which operate by scaling the velocities, such as the
Berendsen thermostat, the stochastic velocity rescaling approach of Bussi and co-workers, and the Nosé-Hoover thermostat, yield
transport properties that are statistically indistinguishable from those of the microcanonical ensemble, provided they are applied
globally, i.e. coupled to the system’s kinetic energy. When coupled to local kinetic energies, a velocity scaling thermostat can have
dampening eﬀects comparable to a velocity randomizing method, as we observe when a massive Nose-Hoover coupling scheme
is used to simulate water. Correct dynamical properties, at least those studied in this paper, are obtained with the Berendsen
thermostat applied globally, despite the fact that it yields the wrong kinetic energy distribution.
1. INTRODUCTION
Temperature control algorithms are an important component
of many molecular dynamics simulations. Using a method to
enforce constant temperature is necessary to compare
simulation results with laboratory experiments conducted at
constant temperature and either constant pressure or
volume.1−3
Precise control of the temperature is necessary
when studying temperature-dependent properties and phenom-
ena,4,5
as a constant energy simulation may not have either the
correct average kinetic energy or correct distribution of kinetic
energies. Controlling the temperature can prevent the system
from heating signiﬁcantly over long time scales (though good
integration algorithms should already minimize this drift) and
can improve the eﬃciency of conformational sampling in some
cases.6
A molecular dynamics thermostat couples a ﬁctitious heat
bath to the system or some portion of the system, such that the
time-averaged instantaneous kinetic energy of the coupled
degrees of freedom corresponds to a target temperature. This
coupling is commonly achieved by altering Newton’s equations
of motion or the particle velocities themselves during the
simulation, so that the kinetic energy is steered toward its
target. For a thermostat to be entirely consistent with the
canonical ensemble, it should generate total energies according
to the Boltzmann distribution for that system and generate
kinetic energies consistent with the Maxwell−Boltzmann
distribution. Additionally, the dynamics should be ergodic,
meaning all states with the same energy are visited with equal
probability in the long time limit. However, this thermodynamic
speciﬁcation of the correct distribution of energies for
the NVT ensemble does not deﬁne what the dynamics of the
system should be. There are many diﬀerent equations of
motion to choose from for a classical molecular system that
result in the same correct thermodynamics.
If an accurate representation of transport properties or other
dynamical quantities is required, then the thermostat should,
besides just maintaining the correct kinetic energy distribution,
operate in such a way that minimally disturbs the Newtonian
dynamics. Since the equations of motion are altered when using
a thermostat, some disturbance of the dynamics must occur,
though the extent of such a disturbance is in general not
known.
Received: February 13, 2013
Published: May 17, 2013
Article
pubs.acs.org/JCTC
© 2013 American Chemical Society 2887 dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−2899
However, some common algorithms are known to cause
signiﬁcant dampening of the energy ﬂuctuations or dynamics of
the system. As the ﬁrst problem of obtaining the correct energy
ﬂuctuations is generally well understood,7
this study deals
primarily with the latter problem. Speciﬁcally, we investigate
the eﬀects of diﬀerent temperature control schemes on basic
dynamical properties relevant to experiments by studying three
test systems: a pure solvent, a solvated small molecule, and a
solvated polymer chain. These properties illustrate the eﬀects
on dynamic processes that occur over diﬀerent time scales. The
tested temperature control schemes include diﬀerent thermostats,
coupling strengths, and degrees of freedom to which the
thermostat is applied.
Section II of this paper brieﬂy reviews the diﬀerent
thermostat algorithms we consider as well as related research
on how temperature control inﬂuences dynamics. In Section III
we describe the details of the simulations and the transport
properties calculated. In Sections IV and V we discuss the
eﬀects of the various temperature coupling schemes, particularly
with regard to the dynamics, and consider implications of
these ﬁndings.
2. BACKGROUND
2.1. Thermostat Algorithms. A variety of methods have
been developed to perform isothermal MD simulations. Here
we brieﬂy describe the thermostat algorithms tested in this
study and particular aspects of their implementation in the
GROMACS simulation package used throughout the work.8
Conceptually, the use of any thermostat is analogous to
coupling a ﬁctitious heat bath to the system in order to
maintain the average temperature T at some target T0. This
heat bath is applied to a given particle i, either by directly
altering the particle’s velocity or by modifying Newton’s
equation of motion:
=m
d
dt
r
F r( )i
i
i i
2
2 (1)
In the absence of temperature control, the system evolves
classically according to eq 1 and has the microcanonical (NVE)
distribution of energies. This ensemble provides the “true”
dynamics, i.e. classical Newtonian dynamics for a system
described by a given force ﬁeld, at the accuracy level
determined by the force calculations and integration algorithm
used. In this study we therefore compare dynamical properties
obtained from Newtonian (constant energy) dynamics with
those from NVT simulations with diﬀerent types of temperature
control.
The coupling strength τT determines how strongly the heat
bath is applied to the system or how quickly deviations in T are
“pulled back” to T0. While the physical meaning of τT varies
depending on the algorithm, in general T is more tightly
controlled as τT → 0 and less tightly controlled as τT increases,
with the system approaching the NVE ensemble in the limit τT
→ ∞. If an algorithm implements the canonical distribution
correctly, τT will not aﬀect the distribution of kinetic energies
sampled, only the rate of change of the instantaneous kinetic
energy.
Among the six thermostat algorithms that we test are
methods which randomize particle velocities in order to control
the temperature. In the Andersen thermostat, at each time step
in the simulation a subset of atoms are randomly selected and
reassigned new velocities chosen from the Maxwell−Boltzmann
distribution corresponding to the target temperature.9
Each
particle experiences stochastic collisions which occur on
average every τT and evolves according to Newton’s equation
of motion in between collisions. Δt/τT (where Δt is the time
step) is the probability that any given particle is selected for
velocity reassignment at a given time step or, eﬀectively, the
size of the subset of particles to be randomized at a given step.
Tighter temperature control is obtained when the chance for a
particle’s velocity to be randomized is large, corresponding to
low values of τT. We also test a “massive Andersen” thermostat
in which the velocities of all particles are randomized at
intervals of τT. Both Andersen schemes generate energy
distributions consistent with the canonical ensemble7,10
and
require no direct modiﬁcation of the integration equations
themselves.
Another stochastic approach to temperature control uses the
Langevin equation of motion to describe the evolution of the
system8
ξ= − + + ̂m
d
dt
m
d
dt
r
r r
F r( )i
i
i i
i
i i i
2
2 (2)
where ξi is a friction constant, typically though not necessarily
the same for all particles, and r̂i is a stochastic noise term
proportional to ξi and the target T0. Here τT is inversely
proportional to the friction constant. Thus, as τT → 0 the
temperature is more tightly controlled due to both greater
friction (removing more energy) and stochastic noise (adding
more energy) in the system. With a carefully chosen coupling
strength, Langevin dynamics may be used to approximate
solvent interactions in an implicit solvent environment, as the
stochastic forces on the molecules simulated in vacuum
represent the friction and random thermal noise of collisions
with an explicitly represented solvent.6
Langevin dynamics has
been shown to be ergodic for a wide range of systems.11
In contrast to approaches that randomize the velocities, the
temperature may be controlled by scaling the velocities at
regular intervals. Here the equation of motion is
γ= − +m
d
dt
m
d
dt
r r
F r( )i
i
i
i
i i
2
2 (3)
which is of a similar form as eq 2, only with no stochastic noise
term, and a term γ that is not restricted to positive values (as is
the friction constant in eq 2) and is usually the same for all
degrees of freedom. γ as well as the velocity scaling factor, λ, are
both functions of τT
−1
; therefore, the velocities are scaled more
signiﬁcantly to achieve tighter temperature control when τT →
0. The magnitude of the velocity scaling is generally low,
however, because τT typically describes the time constant for
relaxation of the ensemble-averaged temperature (or kinetic
energy) of the entire system to the target. With a suﬃciently
large number of degrees of freedom added together, the
temperature varies slowly, and therefore the required velocity
scalings are small.
The original “weak coupling” scaling algorithm devised by
Berendsen uses an exponential decay of the temperature12
τ= −−dT
dt
T T t[ ( )]1
0
(4)
where the time constant τ = 2CVτT/(Nf kB), with CV
representing the constant volume heat capacity, kB representing
Boltzmann’s constant, and Nf representing the number of
degrees of freedom of the system. However, this approach
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992888
yields an energy distribution with lower variance than the
distribution of the true canonical ensemble, because it
disproportionally samples kinetic energies closer to T0 than
would be observed in the Maxwell−Boltzmann distribution.7
Bussi and co-workers developed an extension to the method
which addresses this distribution problem.13
We refer to this
method as a “stochastic rescaling” thermostat to distinguish it
from the other velocity scaling algorithms, as it not only
rescales the velocities similarly to Berendsen but also includes a
stochastic term in the equation describing the decay of the
kinetic energy K
τ τ
= − +dK K K
dt K
N
dW
( ) 2
T f T
0
0
(5)
where K0 is the kinetic energy corresponding to the target
temperature T0 such that K0 = (NfkBT0)/2, the instantaneous
temperature T is related to the instantaneous K by the same
relationship, Nf is again the number of degrees of freedom, and
dW is a Wiener process noise term. Including this stochastic
noise yields the correct energy ﬂuctuations for the canonical
ensemble.
The widely used Nosé-Hoover and Nosé-Hoover-chains
thermostats utilize an extended system to relax the temperature
to T0 in an oscillatory manner.14−16
The Nosé-Hoover method
is also a scaling thermostat, as the velocities are scaled by a
variable pη = Qdη/dt describing the momentum of a heat bath η
coupled to particle i and the instantaneous temperature of the
system, in the following coupled equations of motion for r and
η
η
= −
=
−
ηd
dt m
p
Q
d
dt
d
dt
T T
Q
r F r r( )i i i
i
i
2
2
2
2
0
(6)
where Q is the equivalent of the mass of the heat bath. The
equations give the correct kinetic energy distribution for any
arbitrary choice of Q, though lighter masses lead to faster
oscillatory behavior. In the GROMACS implementation, Q =
(τT
2
T0Nf)/(4π2
), where Nf is the number of degrees of freedom
to which the thermostat is applied, and the coupling strength τT
describes the period of oscillations of kinetic energy transferred
between the system and the heat bath. With Nosé-Hoover
integration the conserved quantity can be obtained from an
extended Lagrangian and is a function of η, pη, and Q.
The Nosé-Hoover thermostat has a drawback in that it is
provably nonergodic and therefore does not sample the
canonical ensemble for certain systems,10
though it appears
to be statistically indistinguishable from the canonical
distribution for molecular systems when particle number is in
the hundreds or larger.7
Martyna et al. proposed an extension,
known as the Nosé-Hoover-chains method, in which the heat
bath η1 in contact with the system is itself in contact with a
chain of NC additional heat baths.16
The equations of motion
are therefore coupled to a chain of NC additional bath variables
increasing the computation time and complexity but increasing
the ergodicity of the dynamics by increasing the available phase
space of the dynamics.
Thermostats (like Berendsen, Nosé-Hoover, and Bussi’s
algorithm) which scale the velocities based on a single bath
variable shared between all degrees of freedom can leak kinetic
energy from high frequency modes into lower frequency
modes, especially the translational or rotational degrees of
freedom which are decoupled from the rest of the system. This
is another possible way that the ergodic assumption can break
down when dynamics are modiﬁed by thermostats. This
problem most egregiously manifests itself in the “ﬂying ice
cube” problem,17
where all the kinetic energy available to the
system moves into translation and rotational motion. The
general solution is to subtract out the kinetic energy of these
degrees of freedom. There is some evidence that kinetic energy
may accumulate in other low frequency vibrational modes,17,18
but the general scope of the problem appears not to have been
studied signiﬁcantly since then, remaining an open question.
Bussi et al. did ﬁnd that using their algorithm, the phonon
distribution of ice remained the same between NVE and NVT
simulations,13
indicating that at least in standard usage, this
problem may not be signiﬁcant for stochastic velocity rescaling.
2.2. Thermostats and Dynamics. As described above,
diﬀerent thermostat algorithms have various inﬂuences on the
classical equations of motion. These inﬂuences lead to diﬀerent
eﬀects on the dynamics of a simulation, which we aim to
quantify in this study. The dampening of dynamical processes
such as diﬀusion can be attributed to disturbances of the natural
time correlations of a particle’s velocity, which can be identiﬁed
by calculating the velocity autocorrelation function (ACF).
While any algorithm which alters the classical equations of
motion will have some eﬀect on velocity time correlations and
thus dynamical properties, these eﬀects can be negligible
depending on the algorithm and the coupling strength at which
the thermostat is applied. Observed dampening in an ensembleor
time-averaged transport property is indicative of consistent
rapid decorrelation of velocities, throughout the system or over
the course of the simulation.
Thermostats which operate by randomizing velocities, such
as the Andersen and Langevin dynamics methods, signiﬁcantly
disturb velocity time correlations if the coupling is suﬃciently
strong. The decorrelation of velocities is an obvious
consequence of velocity randomization, whether by randomly
reassigning velocities in the case of Andersen or by imposing
stochastic forces on the particles in the case of Langevin
dynamics. While in the limit of weak coupling (large τT), the
Andersen thermostat is expected to have a small eﬀect on
system dynamics, it has still been advised not to use such a
stochastic method when studying dynamical properties.10
We hypothesize that typical implementations of velocity
scaling thermostats generally preserve the correct velocity time
correlations, which should be reﬂected in the dynamical
properties we measure. This assumes the scaling algorithm is
implemented “globally”, i.e., by coupling the scaling factor or
heat bath variable for each particle to the system’s overall
kinetic energy (or temperature). As previously mentioned, the
overall kinetic energy typically exhibits relatively slow
ﬂuctuations at equilibrium, and therefore only small changes
to the velocities should be necessary to relax this quantity to its
target. However, coupling a more rapidly ﬂuctuating energy
(e.g., an individual particle’s energy) to the scaling factor or
heat bath variable of a given particle could lead to larger
changes in the velocity, which could in turn aﬀect the velocity
time correlations.
This sort of coupling is used in “local” implementations of a
thermostat, which apply independent thermostats to individual
molecules or atoms that respond to the kinetic energies of the
individual particles, not the total kinetic energy. They can even
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992889
be applied to individual degrees of freedom,19
although single
atoms are the smallest units that can be acted upon individually
in GROMACS. The independent coupling of all atoms or
degrees of freedom individually is often referred to as a massive
thermostat scheme. This use of “massive” is diﬀerent than for
the previously mentioned massive Andersen algorithm, which
stochastically randomizes all velocities simultaneously, independently
of the kinetic energies of both individual particles
and the overall system. Similarly, the conventional Andersen
thermostat is not local by this deﬁnition because alterations of
particle velocities are again independent of the kinetic energy.
Langevin dynamics is also not local in the same way as the
scaling-based methods; although the frictional term is proportional
to the particle velocity, the random term is not, and the
friction term can never increase the velocity. Massive NoséHoover
coupling is sometimes used in practice, particularly for
ab initio molecular dynamics, as it can be an eﬀective approach
for equilibrating a simulation and improving ergodicity.19,20
Finally, we note that using an extremely strong coupling for a
global velocity scaling thermostat (corresponding to shorter
time scales than examined in this paper) could possibly cause
rapid velocity decorrelation as well, if imposing an overly fast
relaxation of the system’s kinetic energy to the target led to
suﬃciently large changes in the velocities.
We have thus identiﬁed two diﬀerent factors which could
rapidly decorrelate velocities and lead to dampened dynamics:
the direct randomization of velocities and the scaling of
velocities if implemented in a local coupling scheme or with
very strong coupling, such that velocities change rapidly when
the kinetic energy deviates from its target. These factors will be
further discussed in light of simulation results obtained using
various thermostat schemes.
2.3. Thermostat Schemes which Treat Separate
Groups of Atoms Diﬀerently. In addition to the eﬀects of
the thermostat algorithm, we are also interested in how a
solvated molecule’s dynamics may be aﬀected if temperature
coupling is applied only to the solvent. This scheme may be
viewed equivalently as using separate thermostats for the solute
and solvent, with inﬁnite and ﬁnite coupling strengths,
respectively. Lingenheil and co-workers have termed this the
“noninvasive” coupling scheme and have also considered a
minimally invasive scheme in which τT applied to the solute is
much larger than the time constant describing the relaxation of
the separately coupled solute’s T to its steady state, as measured
from simulation.21
Note that these schemes could be
implemented either “globally” or “locally”, in the sense that
the terms are used above; the separate temperature control
applied to the solvent or any subsystem may involve a single
thermostat for each particle in the subsystem or alternatively
independent thermostats for each particle.
The non- and minimally invasive schemes as well as a
commonly used approach in which separate thermostats are
applied with a small τT for both subsystems22−24
are well suited
for systems with inhomogeneous rates of heating. In these
cases, applying a single thermostat to the entire system can give
rise to nonuniform temperature distributions, otherwise known
as the “hot solvent-cold solute” problem, due to the fact that
kinetic energy can ﬂow between the diﬀerent degrees of
freedom through the heat bath variables and create a steady
state distribution inhomogeneous in temperature. Using a nonor
minimally invasive coupling scheme can circumvent this
problem and avoid artifacts in the solute dynamics associated
with the inhomogeneous heating, as long as the heat transfer
rate due to molecular collisions between solvent and solute is
non-negligible.
In principle, the most physically realistic way to control a
system’s temperature is to do so exclusively through thermal
contact with the surroundings, using the unaltered Newton’s
equation of motion. If the system is suﬃciently large, the
solvent will act as a perfect heat bath for the solute, and the
solute must therefore have a canonical distribution of
energy.21,25
This type of temperature control would be diﬃcult
to achieve with current MD simulations, however. A nonperiodic
box whose edges were in physical contact with walls of
a ﬁxed temperature would represent a truly physical thermal
bath, but in most cases the number of solvent particles required
to prevent the edges from inﬂuencing the solute would be
prohibitively large. Alternatively, a thermostat might be applied
only to molecules at the edges of a periodic box, but again very
large system sizes would be required to prevent solutes from
being directly aﬀected by the thermostatted regions, particularly
if multiple solutes are included. This scenario also would
introduce nonphysical, anisotropic heat source/sink terms into
the system.
The noninvasive coupling scheme in which only the solvent
is thermostatted can be seen as the next best approach if the
integration is suﬃciently accurate that the solute’s energy does
not drift signiﬁcantly in the absence of direct temperature
control. A solute with unaltered dynamics is coupled to a heat
bath (in this case, the bulk solvent) through molecular
collisions, not through a extradimensional variable (velocity
scaling methods), through collisions with phantom particles
(the Andersen variants), or through both (Langevin dynamics).
By thermalizing only the solvent, we are closer to the “natural”
dynamics of a solute in contact with a solvent, which is itself in
contact with the heat reservoir. If heat is produced as a result of
molecular motion, it may lead to spatially inhomogeneous
temperature transients through the system; but of course this is
closer to “natural” heating behavior of the system. In this study
we test the extent to which removing temperature coupling of
the solute can also reduce artifacts in its dynamics due to the
thermostat algorithm.
2.4. Previous Findings. This paper aims to draw consistent
comparisons between the eﬀects of diﬀerent temperature
coupling schemes on common dynamical properties. Bussi
and co-workers have considered the dynamics of homogeneous
systems using their stochastic velocity rescaling thermostat used
for temperature control. They found that the diﬀusion of
TIP4P water using this thermostat with a range of coupling
strengths agreed with the result of an NVE simulation. The
thermostat also yielded the correct features of the vibrational
spectrum of hydrogen atoms in ice.13
In simulations of a
Lennard-Jones ﬂuid it was also found that diﬀusion was
independent of τT using stochastic velocity rescaling but
dampened with Langevin dynamics at low τT.11
We expect to
see similar trends for pure ﬂuid self-diﬀusion and solute
diﬀusion in our study and hypothesize that other dynamical
properties will be aﬀected similarly.
Another study assessed the “eﬃciency” of diﬀerent thermostats
by comparing the error in the dynamics resulting from the
thermostat to the characteristic time for equilibration of the
kinetic energy.26
In simulations of a homogeneous system,
stochastic velocity rescaling and a modiﬁed Nosé-Hoover
approach that included stochastic inﬂuences were both found
to be more eﬃcient thermostats than Langevin dynamics, and
they perturbed the dynamics to a lesser extent. The dynamics of
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992890
a pure SPC/E water system were considered in another study,
which measured the power spectra of ﬂuctuating energies.18
It
was found that weak (τT > 1 ps) Berendsen coupling yielded
power spectra of potential energy ﬂuctuations which matched
those from the NVE ensemble, while strong Berendsen
coupling led to signiﬁcant deviations at low frequencies.
Another study considered protein dynamics through the
power spectra of the velocity ACFs of α − C atoms and the
ACFs of dihedral ﬂuctuations, using diﬀerent temperature
control methods: Langevin dynamics, a global Nosé-Hoover
thermostat, and global as well as local Nosé-Hoover-chains
algorithms.19
The various Nosé-Hoover and Nosé-Hooverchains
methods yielded ACFs with the same qualitative features
as the NVE results (though peak intensities varied), but these
dynamical features were not preserved using Langevin
dynamics.
Dissipative particle dynamics (DPD) is another common
simulation approach.27
As with Langevin dynamics, it includes
friction and noise terms that when properly balanced provide a
means for temperature control. However, with DPD these
forces are applied to particle pairs rather than to individual
particles, thus preserving the total linear momentum of each
pair. Goga and co-workers recently applied the thermostat
component of DPD to various molecular systems using an
implementation in which the system is evolved classically for
one step, followed by an impulsive application of friction and
noise.28
They compared how diﬀusion was aﬀected by Langevin
dynamics versus this DPD thermostat approach, with the
pairwise forces applied isotropically, or perpendicular or parallel
to the direction of the velocity diﬀerence of a given particle pair.
In each of these three cases, diﬀusion using DPD was strongly
dependent on the eﬀective friction rate (related to our coupling
strength τT), despite the modiﬁcation of the equations of
motion to conserve pairwise linear momentum. Lowe also
showed that the DPD approach enhances the shear viscosity,
which is useful when applying DPD with smoothed potentials
with artiﬁcially low viscosity to make such simulations more
physical, but in our context, is more problematicDPD
modiﬁes and slows the underlying Newtonian dynamics.29
Existing results for DPD are not entirely unambiguous; in one
study of polymer brushes under ﬂuid ﬂow, the correct dynamic
ﬂuctuations of the polymers were obtained with DPD
temperature control but not with Langevin dynamics,
reportedly because of the momentum conservation of the
former.30
However, it is clear that DPD thermostats do not
inherently preserve transport properties.
A thermostat’s impact on the thermodynamic ensemble can
also be important in order to eﬃciently sample complex
systems and obtain the correct distribution of conformations. It
has been rigorously shown that the Berendsen thermostat does
not generate the true canonical ensemble, with too narrow a
distribution of kinetic, potential, and total energies.7
An
incorrect ensemble of energies can lead to systems being
trapped in local energy minima, as seen in one study which
found the time required for a protein to fold increased with
strong Berendsen coupling.31
In replica exchange molecular
dynamics simulations, the Berendsen thermostat has been
shown to alter the distribution of folded and unfolded
conformations of diﬀerent peptides at low temperatures.32,33
Kinetic processes which depend on the thermodynamics of the
system can be similarly aﬀected. Lingenheil and co-workers
simulated a peptide with separate strongly coupled thermostats
for the solute and solvent as well as with their minimally
invasive scheme in which the solute coupling is relaxed.21
Using
the Berendsen thermostat for the solute, the peptide backbone
exhibited dampened dynamics (as seen in a lower ﬂip rate of
the torsional angles) with strong coupling, relative to the
minimally invasive scheme. No dampening was observed with
the Nosé-Hoover thermostat applied in the same manner,
suggesting that the correct energy ﬂuctuations (which are
suppressed by the Berendsen method) are necessary to allow
the energy barrier crossings associated with certain kinetic
processes.
These studies illustrate that the appropriate temperature
control scheme largely depends on the processes or
phenomena of interest. In this study, we test the commonly
used thermostats and coupling strategies described above, in
order to provide a consistent comparison of how various
thermostatting schemes aﬀect transport properties and the
thermodynamic ensemble, which can in turn aﬀect kinetics.
Understanding the extent to which thermostats may alter these
aspects of a simulation should help improve the analysis of MD
studies at constant temperature.
3. METHODS
3.1. Simulation Details. We simulate three test systems, a
pure solvent, a solvated small molecule, and a solvated
uncharged polymer chain with a variety of temperature control
schemes. Water was used as the solvent in each system. The
small molecule system includes a single united-atom (UA)
methane, while the polymer system includes a single coarsegrained
chain comprised of N = 100 united atom monomers.
Table 1 provides additional details on the systems, which are
studied in a cubic box with periodic boundary conditions, and
the number and length of simulations. The GROMACS
molecular dynamics package (a prerelease version of
GROMACS 4.6) is used for these studies.8
We use the
standard 6-12 Lennard-Jones potential to describe all van der
Waals interactions and the TIP3P water model for the solvent,
with water model constraints maintained using the SHAKE
algorithm.34,35
The Lennard-Jones radius and well depth for the
UA-methane are obtained from the OPLS force ﬁeld.36
For the
polymer chain, these parameters are obtained from the TraPPE
united atom force ﬁeld for hydrocarbons.37
For the middle and
terminal monomers we use the parameters for united-atom
CH2 and CH3 in an alkane, respectively. While the TraPPe
force ﬁeld treats bonds with ﬁxed lengths and bending and
torsional potentials, we use only a harmonic potential to
describe the bonds, with an equilibrium distance and spring
constant corresponding to aliphatic sp3
hybrid carbons.38
This
simpliﬁed model leads to a very ﬂexible chain which does not
accurately describe a real alkane but serves as a simple model
for measuring the eﬀects of thermostats on collective polymer
dynamics.
For all systems and integrators the simulation time step is 2
fs. The Particle-mesh Ewald (PME) method is used to calculate
Table 1. Simulation Details for Test Systems
system Nsolvent
box length
(Å)
number of
sim.
length (ns) of
sim.
pure solvent 895 3.0258 1 10
solvated UA-
methane
893 2.9996 1 200
solvated polymer
chain
10192 6.7500 15 1
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992891
Coulombic forces, along with a switch potential with rswitch =
0.89 nm and rCoulomb = 0.9 nm. van der Waals forces are also
treated with a switch potential with rvdw−switch = 0.8 nm and rvdw
= 0.9 nm, and the neighbor list cutoﬀ is rlist = 1.1 nm. We use
the velocity Verlet integration scheme39
for all simulations
except when Langevin dynamics provides the temperature
control, where the leapfrog integration scheme40
is used.
We simulate each system with six diﬀerent thermostats with
target temperature T0 = 300 K as well as the NVE ensemble (or
coupling with τT → ∞). For each thermostat we test the
coupling strengths τT = 0.1, 1.0, and 10.0 ps, values which were
kept constant for the sake of consistency, despite the somewhat
diﬀerent physical meaning of τT for the various algorithms. In
the small molecule and polymer systems, for each thermostat
and τT we consider the two coupling schemes described
previously, in which the thermostat either is applied to the
entire system or applied only to the solvent particles, as
proposed in the “noninvasive” scheme of Lingenheil and co-
workers.
The tested thermostats include three “velocity randomizing”
algorithms: the Andersen thermostat, a massive Andersen
thermostat which randomizes the velocities of all particles, and
Langevin dynamics. We also test three “velocity scaling”
algorithms: the Berendsen thermostat, Bussi and co-workers’
stochastic rescaling thermostat, and the Nosé-Hoover thermostat.
For the simulations of the pure solvent, we also evaluated
the Nosé-Hoover chains approach with NC = 50 chains.
However, the results are statistically indistinguishable from
those obtained using only a single bath variable, and therefore
we report only the Nosé-Hoover results here.
For every system, the initial conﬁgurations were obtained by
running an energy minimization followed by MD simulation for
1 ns with constant temperature and pressure (stochastic
velocity rescaling thermostat with τT = 0.1 ps and T0 = 300 K;
Berendsen barostat with τP = 0.5 ps and Pref = 1 bar). For the
NVE simulations, initial kinetic energies were chosen so that
the steady state average kinetic energy was within 1% of the
target NVT kinetic energy. For the pure solvent and small
molecule systems, all simulations are started from the same
initial conﬁguration, followed by a 100 ps equilibration period,
long enough for the chosen thermostat to determine the
dynamics.
As shown in Table 1, a single extended simulation is used to
calculate the pure solvent and small molecule transport
properties. Consecutive sections of the trajectories are analyzed
because the properties are not signiﬁcantly correlated from one
section to another. Because the dynamics of the polymer chain
evolve more slowly, however, we analyze 15 independent 1 ns
simulations to avoid biased estimates of means and
uncertainties due to correlated samples. Prior to collecting
data for each independent polymer simulation, we ran for 200
ps with stochastic velocity rescaling (τT = 0.1 ps) to obtain a
polymer conﬁguration uncorrelated with the previous structure,
as measured by the root mean squared deviation (RMSD) of
the conﬁguration relative to an initial structure. This was
followed by 100 ps of equilibration with the desired
temperature coupling scheme, to allow its eﬀects on the
dynamics to become dominant.
3.2. Dynamical Properties. We measure a number of
diﬀerent transport properties to investigate how dynamics
occurring over diﬀerent time scales are aﬀected by temperature
control. The translational diﬀusion constant of a given molecule
i is calculated from its average mean square displacement
(MSD) using the following Einstein relation:41
= ⟨| − | ⟩tD tr r2
1
3
( ) (0)i i
2
(7)
We evaluate the MSD of the molecule’s center of mass,
averaged over time for all systems as well as over the ensemble
of molecules for the pure water system. For the pure water and
small molecule, D is measured from the slope of the MSD
between t = 4 and 20 ps, a region which avoids the ballistic
regime at short times and noise due to insuﬃcient sampling at
long times.42
We average the diﬀusivities measured from
consecutive 100 ps sections of the 10 ns trajectory to obtain our
overall estimate of D and estimate the standard error by
bootstrapping the standard deviation of the mean (using 1000
bootstrap samples).
For polymer diﬀusion, to avoid correlated samples we
calculate only one MSD for each independent 1 ns trajectory
and average the diﬀusivities from these 15 sections to obtain D.
Error estimates for the polymer diﬀusivities are comparable
with those of the other systems despite the smaller number of
samples, because analyzing longer trajectories reduces the
uncertainty in each MSD curve.
For the pure solvent, we also calculate the rotational
correlation time, τα, describing the time required for a water
molecule to rotate about a given axis α. Here we consider
rotation about the axis normal to the plane spanned by the two
O−H bonds. The correlation time is estimated by integrating
the following rotational correlation function, calculated from
consecutive 100 ps sections of the trajectory
= ⟨ · ⟩α α α
C t P te e( ) ( ( ) (0))1 1 (8)
where P1 is a ﬁrst order Legendre polynomial, and e is the unit
vector along the α axis. Because the tail of this autocorrelation
function converges slowly, we estimate the area from t = 0 to 5
ps by direct integration of eq 8, and from t = 5 to 20 ps by
integrating the following exponential ﬁt to the data:43
=α τ− α
C t( ) exp t
1
/
(9)
As before, we estimate τα by taking the mean of the correlation
times from the 100 ps subsections of the trajectory and its
associated error by bootstrapping the standard deviation of the
mean (1000 bootstrap samples).
The shear viscosity of the pure ﬂuid is calculated using the
Green-Kubo relation44
∫η = ⟨ + ⟩αβ αβ
∞V
k T
P t t P t dt( ) ( )
B
t
0
0 0 0
(10)
where T is the average instantaneous temperature, and Pαβ are
the oﬀ-diagonal elements of the pressure tensor, and the
integrand is averaged over multiple reference times t0. The ACF
of Pαβ is calculated from t = 0 to 1.8 ps. To integrate this slowly
converging function, we estimate the area under the curve
directly from t = 0 to 0.4 ps, and the area under an exponential
ﬁt to the data from t = 0.4 to 1.8 ps. We evaluate eq 10 for each
of the three oﬀ-diagonal elements of the pressure tensor, for the
four consecutive 2500 ps sections of the trajectory, and estimate
η and the associated error by averaging these 12 data points and
bootstrapping the standard deviation of the mean (1000
bootstrap samples). Error estimates did not vary signiﬁcantly
with the length of the trajectory used to calculate the ACF;
however, a slight minimum in the standard error was obtained
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992892
using 2500 ps sections. It should be noted that while
nonequilibrium methods are often used to compute the shear
viscosity with good accuracy,44
the Green-Kubo relation is
more appropriate for studying how the viscosity at equilibrium
is aﬀected by a thermostat.
For the polymer chain, we also calculate the dynamic
structure factor S(k,t), which in physical systems can be
measured by dynamic light scattering experiments45
∑=
| + − |
| + − |
S k t
N
k t t t
k t t t
r r
r r
( , )
1 sin[ ( ) ( ) ]
( ) ( )ij
i j
i j
t
0 0
0 0
0 (11)
where N is the monomer-length of the chain, k is the
magnitude of the scattering vector (with units of nm−1
here),
and ri(t) is the position of monomer i. The magnitude of the
scattering vector as well as the molecule’s radius of gyration Rg
determines the types of motion that may be detected, with
translational motion dominating when kRg ≪ 1 and internal
modes of motion becoming important for kRg > 1. Given our
objective of studying the internal dynamics of the polymer and
the Rg ≈ 0.7 nm observed for our model polymer, we use a
scattering vector k = 2.0σ−1
, which corresponds to the desired
kRg > 1 regime. Additionally, inspection of the log−log plot of
the static structure factor S(k,t = 0) vs k revealed that k = 2.0σ−1
was in the linear region corresponding to the dynamic scaling
regime.
We calculate S(k,t) from t = 0 to 1000 ps, which for our
model polymer is suﬃcient time for the function to decay to 0
under any type of temperature control. This decay occurs over
much longer time scales in other dilute polymer systems.46
The
rapid decorrelation observed here may be attributed to the high
ﬂexibility of our model polymer, which as mentioned does not
include bending or torsional terms.
We aim to quantify the rate of decorrelation to compare the
eﬀects of the various thermostats. It has been found
previously,46
and we also observed here that S(k,t) in the kRg
> 1 regime cannot be accurately described by a single
exponential function but that a stretched exponential is often
appropriate. We therefore ﬁt the data with a stretched
exponential function of the form y = A exp(−(t/τB)C
) and
measure a characteristic decorrelation rate as τB
−1
. For a given
thermostat, we determine the mean S(k,t) by averaging the
functions of the 15 independent simulations and the overall
decorrelation rate by ﬁtting the resulting function. The
associated errors were estimated by bootstrapping the standard
deviation of the mean function and its associated decorrelation
rate, again with 1000 bootstrap samples.
4. RESULTS AND DISCUSSION
4.1. Bulk Water Dynamics. Throughout this section we
refer to the Andersen thermostat, the massive Andersen variant,
and Langevin dynamics as “velocity randomizing” algorithms
and classify the Berendsen thermostat, Bussi and co-workers’
stochastic velocity rescaling thermostat, and the Nosé-Hoover
thermostat as “velocity scaling” algorithms. We ﬁrst consider
how the dynamics of the pure water system are aﬀected by
these six thermostats as well as by other schemes which
illustrate diﬀerent factors that can dampen dynamics in general.
With this system we investigate dynamic processes that occur
on relatively short time scales, by considering the selfdiﬀusivity,
rotational correlation time, and shear viscosity of
water.
Figure 1 shows these transport properties obtained from
NVT simulations using various thermostats and coupling
strengths τT. It is clearly evident that the velocity randomizing
thermostats signiﬁcantly dampen the dynamics of the ﬂuid as
compared to the NVE results, particularly with strong coupling.
As previously mentioned, the NVE simulation is used as a
benchmark for comparing the various thermostats, because in
the NVE ensemble a system evolves according to the “true”
Newtonian equations of motion. With the massive Andersen
algorithm coupled with τT = 0.1 ps, the self-diﬀusion of water is
slower by approximately a factor of 4, and the characteristic
time for rotation and shear viscosity are increased by factors of
2. With the traditional Andersen thermostat and Langevin
dynamics, the dampening eﬀects are less severe but still
signiﬁcant. These eﬀects diminish as the coupling strength of
the velocity randomizing thermostat weakens, and at τT = 10 ps
the properties agree, within error, with those of the NVE
simulation. It is interesting to note that with both the
conventional and massive implementation of the Andersen
algorithm, each particle’s velocity is randomized on average at
the same frequency, yet the massive Andersen method leads to
greater dampening of the dynamics for TIP3P water and for the
Figure 1. Self-diﬀusivity (a), rotational correlation time (b), and shear viscosity (c) of pure TIP3P water system with various types of temperature
control. Velocity randomizing thermostats (upper ﬁgures) with strong coupling slow the diﬀusion and increase the correlation time and viscosity,
while the velocity scaling algorithms do not aﬀect the kinetic properties.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992893
other systems as well. This may be because τT represents the
average time between randomizations, and thus with the
conventional method there are a large number of particles that
continue uninterrupted trajectories for longer than τT.
On the other hand, the velocity scaling algorithms--
Berendsen, stochastic rescaling, and Nosé-Hoover--yield
properties that agree with the NVE values, regardless of τT,
down to 0.1 ps. This is true even with stochastic velocity
rescaling, as has been shown previously for other systems.13
The thermostat’s stochastic component does not inﬂuence
particle velocities directly, and their natural time correlations
are generally preserved. Our results using the velocity scaling
thermostats and NVE are consistent with previous studies of
pure water systems.43,47,48
As previously mentioned, dampened dynamical properties
are evidence of disruptions to the time correlations of velocities
throughout the system and/or over an extended time. The
above results for the Andersen methods and Langevin dynamics
agree with many previous observations that randomization of
velocities leads to dampened dynamics.10,11
Velocity rescaling
algorithms also may aﬀect velocity time correlations, and
therefore dynamical properties, if they lead to large changes in
the velocities at each step. As mentioned, this is a possibility
using a massive temperature coupling scheme in which each
particle is independently coupled to its own thermostat, which
is in turn coupled to the rapidly ﬂuctuating kinetic energy of
that same particle.
To determine if such a scheme could appreciably aﬀect
velocities, such that ensemble-averaged dynamics are dampened,
we ran additional simulations of the pure water system
subjected to massive Nosé-Hoover coupling. An independent
Nosé-Hoover thermostat (with a single Nosé-Hoover chain)
was applied to each particle using the same τT and T0 = 300 K
for each thermostat. The Nosé-Hoover masses were chosen to
be proportional to the number of degrees of freedom that they
couple to (in this case Nf = 3), as suggested by Martyna et al.16
Harmonic bonds with a time step of 0.5 fs were required to
implement this scheme for TIP3P water in GROMACS, in
contrast to the explicit constraints and 2 fs time step used for
the other simulations, because holonoic constraints are not
preserved when diﬀerent atoms in the constraint system are
scaled by diﬀerent amounts. This treatment of bonds changes
slightly the magnitude of the transport properties, and therefore
we normalize the massive Nosé-Hoover results by the
properties obtained from NVE simulations with harmonic
bonds and compare to the NVE-normalized results obtained
using the thermostats with constraints. Deviations from unity
reﬂect the extent to which a given thermostat dampens a
property relative to the NVE dynamics.
Figure 2 shows the normalized transport properties using the
massive Nosé-Hoover scheme as well as Langevin dynamics, a
thermostat which we found signiﬁcantly aﬀects the dynamics.
As shown, the massive Nosé-Hoover with strong coupling does
indeed dampen the dynamics, slowing diﬀusion by nearly a
factor of 2 when τT = 0.1 ps. This shows that dampening due to
a massive implementation of a velocity scaling algorithm can be
comparable to that of a velocity randomizing thermostat. The
diﬀerence between the massive and global Nosé-Hoover
schemes is signiﬁcant, though it is not clear exactly to what
extent this is because individual kinetic energies ﬂuctuate more
than the global kinetic energy or because the thermostat masses
Qi are smaller in the massive scheme. It is possible that similar
eﬀects would be observed with massively coupled Berendsen
and stochastic rescaling thermostats, though implementation
issues in GROMACS prevent directly testing this hypothesis at
this time.
For massive Nosé-Hoover and Langevin dynamics, we
conducted additional simulations with τT = 100 ps to determine
how the dynamics are aﬀected with very weak temperature
coupling. In this limit, both schemes yield dynamical properties
that converge to within 2% of the NVE results. At this τT the
thermostats also maintain the average instantaneous temperature
at the target T0 = 300 K (within statistical error) as
eﬀectively as they do with stronger coupling. Therefore, even
algorithms known to dampen dynamics with strong coupling
can be used in constant temperature simulations in which
relatively short time scale dynamics are important, if one uses a
suﬃciently weak τT that only neglibibly impacts the velocity
time correlations. However, it is possible that dynamical
processes occurring in more complex systems and/or on longer
time scales would be more sensitive to velocity randomization
or massive velocity rescaling, even when weakly coupled. We
note that with the Nosé-Hoover approach, the value of τT
depends on exactly how the heat bath variable Q is deﬁned.
Diﬀerent multiplicative factors of the product τT
2
T0 are
sometimes used,8,19
which aﬀects the ﬁctitious mass coupled
to a given particle and therefore the magnitude of changes
made to the particle’s velocities. We can, however, say that as
the coupling becomes stronger, the dynamics will begin to
diverge from the NVE dynamics.
The consistent dampening eﬀects of the Andersen, massive
Andersen, and Langevin dynamics algorithms clearly show that
Figure 2. Self-diﬀusivity (a), rotational correlation time (b), and shear viscosity (c) of pure TIP3P water using massive Nosé-Hoover coupling and
Langevin dynamics temperature control with various coupling strengths: τT = 0.1 ps (no lines), τT = 1 ps (slanted lines, thin space), τT = 10 ps
(slanted lines, thinner space), and τT = 100 ps (slanted lines, thinnest space). Both thermostat algorithms signiﬁcantly dampen the dynamics when
applied with strong coupling.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992894
the reassignment of particle velocities disturbs their time
correlations. “Randomness” per se is not required, as shown in
Figure 3. Here we plot normalized properties from simulations
conducted in the NVE ensemble, in which the velocities of half
the molecules in the system are swapped with the velocities of
the other half of the molecules (the same ones each time) at
regular intervals speciﬁed by τT. This “NVE velocity-swap”
scheme does not control the temperature at all but serves as a
deterministic velocity randomization scheme, as the velocity
reassigned to a given particle comes from another particle
within the system. Figure 3 shows that this scheme leads to a
dampening of dynamics comparable to that of its stochastic
analog, the massive Andersen thermostat. It is the reassignment
of the particles with velocities uncorrelated to their previous
velocities that leads to the slowed dynamical properties, not
“randomization” per se.
These simulations of a pure ﬂuid system illustrate various
properties of a thermostat which can lead to artifacts in the
dynamics. Direct alteration of velocities, whether by stochastic
or deterministic reassignment or by the use of stochastic forces,
signiﬁcantly dampens dynamical processes occurring on short
time scales. Comparable eﬀects are observed with velocity
scaling algorithms if implemented locally, such that individual
particles’ velocities change rapidly as a function of time. The
eﬀects of both of these factors diminish in the limit of weak
coupling, which corresponds to infrequent randomizations or
slow relaxation of the kinetic energy to its target, for velocity
randomizing or rescaling algorithms, respectively.
4.2. Solute Dynamics. We restrict the following study of
solute dynamics to the six standard thermostats listed above,
with the velocity rescaling algorithms applied globally. While
the thermostats are classiﬁed as either velocity randomizing or
rescaling algorithms, we emphasize that these classes do not
necessarily lead to dampended and correct dynamics,
respectively, but that the eﬀects of a given algorithm will
depend on the way in which it is applied.
Figure 4 shows the translational diﬀusivities for the solvated
UA-methane and polymer chain using various types of
Figure 3. Self-diﬀusivity (a), rotational correlation time (b), and shear viscosity (c) of pure TIP3P water obtained from NVE simulations in which
velocities are swapped at intervals of τT and from NVT simulations with massive Andersen temperature control. The coupling strengths are as
follows: τT = 0.1 ps (no lines), τT = 1 ps (slanted lines, thin space), τT = 10 ps (slanted lines, thinner space). The deterministic and the stochastic
reassignment algorithms signiﬁcantly dampen the dynamics in a similar way when applied with strong coupling.
Figure 4. Diﬀusivities of solvated united-atom methane (a) and solvated polymer chain (b) with thermostats applied to either the entire system
(Sys) or to only the solvent (Solv). Velocity randomizing thermostats (upper ﬁgures) with strong coupling (τT = 0.1 ps) increase the self-diﬀusivity
of water and the diﬀusivities of both solutes tested, with a larger eﬀect observed for the polymer. The noninvasive scheme of coupling only the
solvent slightly reduces the dampening of diﬀusion.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992895
temperature control. With these systems we also compare the
“system-wide” scheme of applying a single thermostat
uniformly to the entire system versus the noninvasive approach
of applying the thermostat only to the solvent. Similar to our
results for the self-diﬀusion of water, we see that, in general,
solute diﬀusion is slowed by the velocity randomizing
thermostats with strong coupling as compared to the NVE
results but is unaﬀected by the velocity scaling methods
regardless of the coupling strength. The strategy of coupling
only the solvent particles slightly mitigates the dampening of
diﬀusion observed with velocity randomizing algorithms. For
example, using Langevin dynamics and τT = 0.1 ps, the
noninvasive scheme increases the diﬀusivities of UA-methane
and the polymer by 25% and 90% as compared to the systemwide
scheme, though even with this increase the diﬀusivities are
well below the NVE results. This suggests that velocity
randomization of the solvent particles dampens not only the
solvent dynamics but also the dynamics of the solute in thermal
contact with the thermostatted solvent. Therefore, regardless of
which degrees of freedom are coupled, when such an algorithm
is used it is important to choose an appropriately weak τT if the
correct solute dynamics are to be preserved.
We also observe that the dampening of UA-methane’s
diﬀusion is comparable to the dampening of the self-diﬀusion of
water, while diﬀusion of the larger polymer molecule is
dampened more signiﬁcantly. Using Langevin dynamics applied
to the entire system with τT = 0.1 ps, the polymer’s diﬀusivity is
lower than the NVE result by approximately 90%, while the
diﬀusivities of both pure TIP3P water and UA-methane are
lower than their respective NVE results by approximately 60%.
Using system-wide Langevin dynamics and weak coupling of τT
= 10 ps, the polymer’s diﬀusivity is still approximately 15%
lower than the NVE value, though it approaches to the NVE
result (within error) when coupling is applied only to the
solvent. This suggests there is a relationship between a
molecule’s size and the extent to which velocity randomization
aﬀects its translational diﬀusion. Therefore, when using a
velocity randomizing algorithm to simulate a large solute at
constant temperature, very weak coupling and possibly relaxing
or removing temperature coupling to the solvent may be
necessary to obtain the correct translational diﬀusion. It is
possible that the sensitivity of other dynamical properties to
velocity randomization also depends on the solute size.
The polymer dynamic structure factor S(k,t) describes a
dynamic process which occurs on longer time scales, the
ﬂuctuations of the chain’s structure. Figure 5a shows
representative S(k,t) functions for NVE and NVT simulations,
the latter with various thermostats applied to the entire system
with τT = 0.1 ps. The scattering vector k = 2.0σ−1
is appropriate
for detecting internal motions of the chain, as kRg > 1. As
shown, the rate of decay of S(k,t) is much slower with velocity
randomizing thermostats than with velocity scaling thermostats
or NVE. The former methods slow down the natural
conﬁgurational ﬂuctuations in the chain, leading to structures
that remain correlated over times of up to 1 ns. As previously
discussed, these time scales are much less than those observed
in other polymer simulations and in experiments, due to the
high ﬂexibility of our model chain. As a check we calculated
another dynamic property of the polymer, the autocorrelation
function of the end-to-end distance of the chain, and found that
it decays on time scales of the same order of magnitude as
S(k,t). Despite the generally rapid ﬂuctuations of our model
chain, the impact of velocity randomizing thermostats on the
collective internal dynamics of a polymer is evident with this
simpliﬁed model. We note that there is good agreement in our
results obtained using the Berendsen and Nosé-Hoover
thermostats, whereas Lingenheil and co-workers observe
dampened peptide dihedral transitions under Berendsen
control.21
One possible interpretation is that the conformational
ﬂuctuations of our model polymer chain involved barriers
that are low enough to be crossed with the kinetic energy
ﬂuctuations observed using Berendsen temperature control,
while the dihedral barriers observed by Lingenheil et al. were
not, but resolving that question is beyond the scope of the
current study.
Upon ﬁtting S(k,t) with the stretched exponential y = A
exp(−(t/τB)C
) we ﬁnd that the parameters A and C vary little
(<10%) with respect to the type of thermostat, but there are
signiﬁcant deviations in the time constant τB. Using τB
−1
to
quantify the rate at which the chain’s structure decorrelates, we
Figure 5. Dynamic structure factor S(k,t) of solvated polymer chain with various thermostats applied to the entire system with τT = 0.1 ps (a).
Velocity randomizing algorithms lead to polymer structures that remain correlated over longer times. The rate of decorrelation, as measured by the
time constant of a stretched exponential ﬁt to S(k,t), is shown for all temperature control schemes tested (b). The magnitude of scattering vector is k
= 2.0σ−1
for all S(k,t) reported.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992896
compare the dynamics with various temperature coupling
schemes in Figure 5b. Again we observe a dampening eﬀect of
slower decorrelation rates due to velocity randomizing
thermostats, which is mitigated somewhat by using weaker
coupling and coupling only the solvent. Decorrelation rates
obtained using the global velocity scaling algorithms agree with
the NVE result, within error. The corresponding characteristic
decorrelation time τB varies from 10 ps with velocity scaling
temperature control and NVE, to approximately 80 and 140 ps
with strongly coupled Langevin dynamics and the massive
Andersen thermostat, respectively.
It is possible that the dynamics of diﬀerent modes of motion
could exhibit varying sensitivities to a thermostat which disrupts
velocity time correlations, though the diﬀerences we observe for
our model polymer chain may be of limited practical
signiﬁcance. Whereas the polymer’s translational diﬀusion
using Langevin dynamics and τT = 1 ps is approximately 50%
slower than with NVE, its structure factor decorrelation rate,
which reﬂects higher order modes of motion, is only 15% lower
than the NVE result. In principle, this suggests that when using
a thermostat that can potentially decorrelate velocities, the
appropriate coupling strength depends on the time scale of the
dynamic process of interest. However, diﬀerences between the
τT values needed to obtain accurate results for one dynamic
process versus another might be too small to justify careful
study of these parameters for various systems, in which case the
conservative choice of a weak coupling strength would ensure
that all dynamics are relatively well-preserved.
4.3. Temperatures and Energy Distributions. Because
the thermodynamics of a system are important in and of
themselves and can inﬂuence kinetic processes as well, we
measure the eﬀects of the various thermostats on the average
and distribution of the polymer chain’s instantaneous temperature.
The average temperature of the polymer is maintained at
or near the target T0 = 300 K by all the thermostats tested, for
both the approaches of coupling the entire system and coupling
only the solvent. The Andersen thermostats with strong
coupling and the Nosé-Hoover thermostat with any τT provide
the most precise and accurate temperature control, with
deviations from T0 of less than 0.1 ± 0.01 K. Precise control
is also obtained with the Berendsen thermostat, stochastic
velocity rescaling, and Langevin dynamics when strongly
coupled; however, deviations of up to 1 ± 0.1 K are observed
using Berendsen and stochastic velocity rescaling with τT = 10
ps. These small but statistically signiﬁcant deviations may be a
result of artifacts observed by Lingenheil et al.21
and may point
to issues using thermostats with a single bath. With coupling
applied only to the solvent, larger deviations of up to 2 ± 1 K
are observed with some algorithms; however, as seen in Figures
4b and 5b, these variations are not statistically signiﬁcant
enough to appreciably aﬀect the polymer dynamics tested here.
Average temperatures for the pure solvent and solvated small
molecule systems are similarly well-controlled.
As previously mentioned, an incorrect energy distribution
can impact certain dynamic processes. We therefore compare
the distributions of the polymer’s instantaneous temperature
from our simulations to the predicted distribution for the true
canonical ensemble. To calculate this prediction, we ﬁrst note
that for a system with kinetic energy (3N)/(2kBT), the
associated variance is σKE
2
= kBT2
CV. If one considers only the
kinetic energy, the ideal gas heat capacity may be used for CV,
which yields a kinetic energy variance of σKE
2
= (3N)/(2(kBT)2
).
We ﬁnd the temperature corresponding to this σKE by noting
that, on average, KE = (3/2)NkBT. Finally, because N of our
polymer chain is suﬃciently large, we approximate the
Maxwell−Boltzmann distribution as a Gaussian with mean T
= 300 K and standard deviation σT = (2σKE)/(3NkB).
In Figure 6 we compare this analytic distribution to those
obtained from the simulations with system-wide and solventonly
temperature coupling, respectively. We also include the
NVE distribution. Figure 6a shows that when the Berendsen
thermostat is applied to the entire system with strong coupling,
the polymer chain samples an overly narrow distribution of
temperatures (and therefore kinetic energies). In more complex
systems, this could be reﬂected in dampened kinetics for
processes that depend on large energy ﬂuctuations.21
On the
other hand, Langevin dynamics and stochastic velocity rescaling
yield distributions that reﬂect the true canonical ensemble, and
the NVE simulation deviates only slightly from the canonical
prediction.
Figure 6b shows the polymer temperature distributions when
only the solvent is thermostatted. Slight deviations in the peak
locations reﬂect the previously discussed variations in the mean.
In general, the results match the Maxwell−Boltzmann
Figure 6. Distribution of the polymer chain’s instantaneous temperature with (a) temperature coupling applied to the whole system and (b)
coupling applied to only the solvent. A canonical kinetic energy distribution for the polymer can be obtained even with NVE if the solvent bath is
suﬃciently large and with Berendsen if only the solvent is thermostatted, as the solvent acts as a thermal bath for which only the average temperature
is important.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992897
distribution for all NVT simulations as well as for NVE.
Agreement with the canonical distribution even in the absence
of direct thermostatting supports previous ﬁndings that a
solvent may eﬀectively control a solute’s temperature by virtue
of its thermal contact with the molecule, acting as a heat bath.21
This is true even for the Berendsen thermostat, which does not
give the proper distribution when applied to the entire system.
However, in the standard derivation of the canonical
distribution, all that is required of the heat bath is that it be
much larger than the system and that its temperature is ﬁxed,
conditions that are met with Berendsen coupled to the solvent
and also with no coupling (NVE) if the solvent bath is
suﬃciently large. In some cases an explicit thermostat, rather
than just a large solvent bath at NVE, might be necessary to
avoid signiﬁcant drift in the solute temperature over long times,
though if the integration is suﬃciently accurate, the solute
heating time scale should be much longer than the time scale of
heat transport between the solute and solvent.
5. CONCLUSIONS
A variety of methods have been developed to run isothermal
molecular dynamics simulations. These algorithms modify
either the equations of motion or the particle velocities directly
and therefore can have unintended eﬀects on the system’s
dynamics and thermodynamics. We have compared how six
diﬀerent thermostats and various coupling schemes inﬂuence
these aspects of a simulation, particularly the dynamics as
measured by fundamental transport properties.
The dampening of dynamical properties reﬂects systemic
decorrelation of particle velocities due to factors external to the
molecular description itself, in this case a thermostat algorithm.
As one would expect, velocity time correlations are disrupted by
thermostats that randomize velocities. These include the
Andersen thermostat, which stochastically reassigns velocities
at regular intervals, and Langevin dynamics, which includes
additional terms for stochastic noise and friction in the
equations of motion. When these algorithms are applied with
strong coupling, we ﬁnd the dynamics are signiﬁcantly
dampened relative to the NVE results, which provide a
measure of the true dynamics obtained with the unaltered
equations of motion. The dampening eﬀects are observed in
dynamics measured over various time scales, including shear
viscosities (100
− 01
ps time scale), translational and rotational
diﬀusion (100
−101
ps), and the dynamic structure factor of the
polymer (100
−103
ps for our model). We ﬁnd that a massive
Andersen approach generally leads to greater dampening of
dynamics than the conventional Andersen thermostat, despite
reassigning velocities at the same rate, and that a deterministic
reassignment strategy leads to similar eﬀects as those observed
with Andersen, which uses stochastic reassignment.
For solvated molecules, the dampening of dynamics is only
somewhat mitigated by the noninvasive strategy of applying the
thermostat to the solvent only. This suggests that a solute’s
diﬀusion and higher-order dynamics can be slowed just by
virtue of its thermal contact with a solvent subjected to velocity
randomization. In general, the eﬀects of a velocity randomizing
thermostat diminish as the coupling strength weakens,
approaching the Newtonian dynamics limit if coupling times
are on the order of 10 ps for these systems. We also observe
that the extent to which these methods slow a solute’s diﬀusion
increases with the size of the molecule and that the extent to
which they dampen a given transport property depends on the
time scale of the dynamical process. Translational diﬀusion of
our model polymer was slowed more signiﬁcantly by velocity
randomization than was the decorrelation rate of the polymer’s
dynamic structure factor.
Algorithms which operate by scaling particle velocities
include the Berendsen thermostat and Bussi and co-workers’
stochastic velocity rescaling thermostat, which adjust velocities
so that the temperature relaxes exponentially to its target at a
speciﬁed rate, and the Nosé-Hoover thermostat, which adjusts
velocities by coupling the equations of motion to extended
dynamical variables and gives an oscillatory relaxation to the
target. When applied globally to the system using any coupling
strength, these thermostats yield transport properties that are
statistically indistinguishable from the NVE results. However, a
massive implementation of the Nosé-Hoover thermostat, in
which each particle is coupled independently to its own
thermostat, leads to dampening on the same order as that
observed using Langevin dynamics. This illustrates that even a
velocity scaling algorithm can signiﬁcantly disturb velocity time
correlations if the scalings are suﬃciently large, for instance if
particles are coupled to their own rapidly ﬂuctuating kinetic
energies with light masses.
It appears that the well-documented problem of the
Berendsen thermostat providing an incorrect distribution of
instantaneous temperature may be largely mitigated by applying
the thermostat to only the solvent, as long as one is only
concerned with the dynamics of the solvated system. In more
complex systems with kinetic processes that depend on energy
ﬂuctuations, this approach could possibly preserve the correct
dynamics. However the same results might be obtained with
stochastic velocity rescaling or the Nosé-Hoover thermostat,
without the need to remove the coupling of the solute.
These ﬁndings suggest that velocity scaling methods are most
appropriate for simulations that require both constant temperature
and realistic dynamics. When applied globally, as is
typically done, these thermostats yield correct transport
properties which describe processes that occur on relatively
short time scales. However, these algorithms do not intrinsically
preserve the correct dynamics; if applied to local kinetic
energies they can lead to dampening eﬀects comparable to
those caused by velocity randomizing thermostats, unless the
coupling is suﬃciently weak.
We ﬁnd that the globally applied stochastic rescaling and
Nosé-Hoover thermostats are particularly eﬀective in preserving
the correct transport properties and thermodynamic
distributions. Velocity randomizing algorithms should be used
with caution when accurate representations of dynamical
processes are important; even if the dynamics are not
important, they will slow down the rate of sampling for the
system. However, their dampening eﬀects may be reduced by
using a suﬃciently weak coupling strength or a noninvasive
coupling strategy, provided the temperature is adequately
controlled. Even when the dynamics of a system are important,
these methods might be useful depending on the objective. For
example, in an implicit solvent simulation, stochastic velocity
randomization serves to model the thermal collisions of the
omitted solvent molecules, and the coupling strength can be
tuned to yield a desired solute diﬀusivity.
It is important to note that this study only examined
subnanosecond kinetics. It remains an open question if lessons
learned here also extend to longer time scale collective
dynamics. It could be that over the nanosecond or microsecond
time scales important in protein folding, velocity scaling
algorithms will also have dampening eﬀects that cannot be
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992898
detected in the simulations performed here, or alternatively that
the distortions of the fast kinetics are unimportant to much
longer dynamical processes that simply require suﬃciently large
ﬂuctuations in the local kinetic energy to activate. However,
given the absence of any statistically noticeable eﬀect of these
algorithms on short time scale kinetics, it is likely worth testing
the hypothesis that these integrators do not signiﬁcantly aﬀect
dynamics on longer time scales.
In summary, we have shown how the fundamental molecular
dynamics practice of controlling the temperature can modify
the dynamics and transport properties in simulations of simple
test systems. Improving our understanding of these eﬀects,
particularly for more complex systems, is an important step
toward correctly choosing thermostats in order to obtain the
best estimates of the properties desired from the simulation.
■ ASSOCIATED CONTENT
*S Supporting Information
Tables S1−S3 give the results for the self-diﬀusivity, rotational
correlation time, and shear viscosity of the pure TIP3P water
system shown in Figure 1. Tables S4 and S5 give the
translational diﬀusion results for the UA methane and polymer
chain systems, shown in Figure 4. Table S6 gives the results for
the decorrelation rate of the dynamic structure factor of the
polymer chain shown in Figure 5b. This material is available
free of charge via the Internet at http://pubs.acs.org.
■ AUTHOR INFORMATION
Corresponding Author
*E-mail: michael.shirts@virginia.edu.
Notes
The authors declare no competing ﬁnancial interest.
■ ACKNOWLEDGMENTS
This work was supported in part by the U.S. National Science
Foundation (Grant No. CBET-1032727). The authors
acknowledge Lillian T. Chong and Daniel M. Zuckermann
for useful discussions.
■ REFERENCES
(1) Berger, O.; Edholm, O.; Jähnig, F. Biophys. J. 1997, 72, 2002−
2013.
(2) Snow, C. D.; Nguyen, H.; Pande, V. S.; Gruebele, M. Nature
2002, 420, 102−106.
(3) Cadena, C.; Zhao, Q.; Snurr, R. Q.; Maginn, E. J. J. Phys. Chem. B
2006, 110, 2821−2832.
(4) Zhang, L.; Greenfield, M. L. J. Chem. Phys. 2007, 127, 194502.
(5) Fan, H.; Mark, A. E. Protein Sci. 2004, 13, 211−220.
(6) Hünenberger, P. H. In Advanced Computer Simulation, Holm, C.,
Kremer, K., Eds.; Adv. Polym. Sci.; Springer: Berlin, Heidelberg, 2005;
Vol. 173; pp 105−149.
(7) Shirts, M. R. J. Chem. Theory Comput. 2013, 9, 909−926.
(8) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem.
Theory Comput. 2008, 4, 435−447.
(9) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384−2393.
(10) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From
Algorithms to Applications; Academic Press: 2002.
(11) Bussi, G.; Parrinello, M. Comput. Phys. Commun. 2008, 179, 26−
29.
(12) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.;
DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690.
(13) Bussi, G.; Donadio, D.; Parrinello, M. J. Chem. Phys. 2007, 126,
014101.
(14) Nosé, S. J. Chem. Phys. 1984, 81, 511−519.
(15) Hoover, W. G. Phys. Rev. A 1985, 31, 1695−1697.
(16) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992,
97, 2635−2643.
(17) Harvey, S. C.; Tan, R. K. Z.; Cheatham, T. E. J. Comput. Chem.
1998, 19, 726−740.
(18) Mudi, A.; Chakravarty, C. Mol. Phys. 2004, 102, 681−685.
(19) Tobias, D. J.; Martyna, G. J.; Klein, M. L. J. Phys. Chem. 1993,
97, 12959−12966.
(20) Morrone, J. A.; Car, R. Phys. Rev. Lett. 2008, 101, 017801.
(21) Lingenheil, M.; Denschlag, R.; Reichold, R.; Tavan, P. J. Chem.
Theory Comput. 2008, 4, 1293−1306.
(22) Weber, W.; Hünenberger, P. H.; McCammon, J. A. J. Phys.
Chem. A 2000, 104, 3668−3675.
(23) Cuthbertson, J. M.; Bond, P. J.; Sansom, M. S. P. Biochemistry
2006, 45, 14298−14310.
(24) Tao, Y.; Rao, Z.-H.; Liu, S.-Q. J. Biomol. Struct. Dyn. 2010, 28,
143−157.
(25) Gallo, M. T.; Grant, B. J.; Teodor, M. L.; Melton, J.; Cieplak, P.;
Phillips, G. N. J.; Stec, B. Mol. Simul. 2009, 35, 349−357.
(26) Leimkuhler, B.; Noorizadeh, E.; Penrose, O. J. Stat. Phys. 2011,
143, 921−942.
(27) Groot, R. D.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423−
4435.
(28) Goga, N.; Rzepiela, A. J.; de Vries, A. H.; Marrink, S. J.;
Berendsen, H. J. C. J. Chem. Theory Comput. 2012, 8, 3637−3649.
(29) Lowe, C. P. Europhys. Lett. 1999, 47, 145−151.
(30) Pastorino, C.; Kreer, T.; Müller, M.; Binder, K. Phys. Rev. E
2007, 76, 026706.
(31) Fuzo, C. A.; Degreve, L. J. Mol. Model. 2012, 18, 2785−2794.
(32) Rosta, E.; Buchete, N.-V.; Hummer, G. J. Chem. Theory Comput.
2009, 5, 1393−1399.
(33) Spill, Y. G.; Pasquali, S.; Derreumaux, P. J. Chem. Theory
Comput. 2011, 7, 1502−1510.
(34) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys.
1977, 23, 327−341.
(35) Andersen, H. C. J. Comput. Phys. 1983, 52, 24−34.
(36) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110,
1657−1666.
(37) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569−
2577.
(38) Bjelkmar, P.; Larsson, P.; Cuendet, M. A.; Hess, B.; Lindahl, E. J.
Chem. Theory Comput. 2010, 6, 459−466.
(39) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J.
Chem. Phys. 1982, 76, 637−649.
(40) Hockney, R. W.; Goel, S. P.; Eastwood, J. J. Comput. Phys. 1974,
14, 148−158.
(41) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids;
Oxford Science Publications: 1987.
(42) Chitra, R.; Yashonath, S. J. Phys. Chem. B 1997, 101, 5437−
5445.
(43) van der Spoel, D.; van Maaren, P. J.; Berendsen, H. J. C. J. Chem.
Phys. 1998, 108, 10220−10230.
(44) Hess, B. J. Chem. Phys. 2002, 116, 209−217.
(45) Binder, K. Monte Carlo and Molecular Dynamics Simulations in
Polymer Science; Oxford University Press: 1995; pp 125−194.
(46) Harnau, L.; Winkler, R. G.; Reineker, P. J. Chem. Phys. 1996,
104, 6355−6368.
(47) Yeh, I.-C.; Hummer, G. J. Phys. Chem. B 2004, 108, 15873−
15879.
(48) Mark, P.; Nilsson, L. J. Phys. Chem. A 2001, 105, 9954−9960.
Journal of Chemical Theory and Computation Article
dx.doi.org/10.1021/ct400109a | J. Chem. Theory Comput. 2013, 9, 2887−28992899