Adv. Polym. Sci. (2005) 173:105–149
DOI:10.1007/b99427
© Springer-Verlag Berlin Heidelberg 2005
Thermostat Algorithms for Molecular Dynamics
Simulations
Philippe H. Hünenberger
Laboratorium für Physikalische Chemie, ETH Zürich, CH-8093 Zürich, Switzerland
phil@igc.phys.chem.ethz.ch
Abstract Molecular dynamics simulations rely on integrating the classical (Newtonian)
equations of motion for a molecular system and thus, sample a microcanonical (constantenergy)
ensemble by default. However, for compatibility with experiment, it is often desirable
to sample conﬁgurations from a canonical (constant-temperature) ensemble instead. A modiﬁcation
of the basic molecular dynamics scheme with the purpose of maintaining the temperature
constant (on average) is called a thermostat algorithm. The present article reviews the
various thermostat algorithms proposed to date, their physical basis, their advantages and their
shortcomings.
Keywords Computer simulation, Molecular dynamics, Canonical ensemble, Thermostat al-
gorithm
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
2 Ensembles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
3 Thermostat Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
3.1 Temperature in the Monte Carlo Algorithm. . . . . . . . . . . . . . . . . . . . . . . . 118
3.2 Temperature Relaxation by Stochastic Dynamics . . . . . . . . . . . . . . . . . . . 120
3.3 Temperature Relaxation by Stochastic Coupling . . . . . . . . . . . . . . . . . . . 122
3.4 Temperature Constraining . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
3.5 Temperature Relaxation by Weak Coupling . . . . . . . . . . . . . . . . . . . . . . . 127
3.6 Temperature Relaxation by the Extended-System Method . . . . . . . . . . . 129
3.7 Generalizations of the Previous Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 137
4 Practical Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
5 Appendix: Phase-Space Probability Distributions . . . . . . . . . . . . . . . . 138
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
List of Abbreviations and Symbols
CM center of mass
MC Monte Carlo
106 Philippe H. Hünenberger
MD molecular dynamics
NMR nuclear magnetic resonance
SD stochastic dynamics
δ Kronecker delta symbol or Dirac delta function
h Heaviside step function
kB Boltzmann’s constant
U instantaneous potential energy
K instantaneous kinetic energy
H Hamiltonian
E energy (thermodynamical)
H enthalpy
L Hill energy
R Ray enthalpy
T temperature (instantaneous)
T temperature (thermodynamical)
To reference temperature (heat bath)
β β = (kBTo)−1
V volume (instantaneous)
V volume (thermodynamical)
P pressure (instantaneous)
P pressure (thermodynamical)
NN number of particles (n-species, instantaneous)
N number of particles (n-species, thermodynamical)
νν chemical potential (n-species, instantaneous)
µµ chemical potential (n-species, thermodynamical)
psys system linear momentum
Lsys system angular momentum
pbox box linear momentum
Lbox box angular momentum
Nd f number of internal degrees of freedom
Nc number of geometrical constraints
Nr number of external degrees of freedom
mi mass of atom i
˙ro
i real velocity of atom i
˙ri peculiar velocity of atom i
Fi force on atom i
pi momentum of atom i
Ri stochastic force on atom i
γi friction coefﬁcient of atom i
λ velocity scaling factor
t timestep
ζT temperature relaxation time
α collision frequency (Andersen thermostat)
τB relaxation time (Berendsen thermostat)
Q “mass” of the time-scaling coordinate (Nosé-Hoover thermostat)
Thermostat Algorithms 107
τN H effective relaxation time (Nosé-Hoover thermostat)
Le extended-system Lagrangian (Nosé-Hoover thermostat)
He extended-system Hamiltonian (Nosé-Hoover thermostat)
Ee extended-system energy (Nosé-Hoover thermostat)
1
Introduction
Classical atomistic simulations, and in particular molecular dynamics (MD) simulations,
have nowadays become a common tool for investigating the properties of
polymer [1] and (bio-)molecular systems [2, 3, 4, 5, 6, 7]. Due to their remarkable
resolution in space (single atom), time (femtosecond), and energy, they represent a
powerful complement to experimental techniques, providing mechanistic insight into
experimentally observed processes. However, direct comparison with experiment requires
that the boundary conditions imposed on the simulated system are in adequation
with the experimental conditions. The term boundary condition is used here to
denote any geometrical or thermodynamical constraint enforced within the whole
system during the simulation. One may distinguish between hard and soft boundary
conditions. A hard boundary condition represents a constraint on a given instantaneous
observable, i.e. it is satisﬁed exactly at any timepoint during the simulation.
A soft boundary condition represents a constraint on the average value of an observable,
i.e. the corresponding instantaneous value is allowed to ﬂuctuate around the
speciﬁed average. The deﬁnition of a soft boundary condition generally also requires
the speciﬁcation of a timescale for which the average observable should match the
speciﬁed value. There exist four main types of boundary conditions in simulations:
1. Spatial boundary conditions include the deﬁnition of the shape of the simulated
system and the nature of its surroundings. In molecular simulations, one typically
uses either: (i) vacuum boundary conditions (solute molecule surrounded
by vacuum); (ii) ﬁxed boundary conditions (solute-solvent system surrounded
by vacuum, e.g. droplet [8, 9, 10, 11, 12, 13, 14, 15]); (iii) periodic boundary
conditions (solute-solvent system in a space-ﬁlling box, surrounded by an inﬁnite
array of periodic copies of itself [16, 17]). In the two former cases, the effect
of a surrounding solvent can be reintroduced in an implicit fashion by a modiﬁcation
of the system Hamiltonian. Typical modiﬁcations are the inclusion of: (i)
solvation forces accounting for the mean effect of the solvent [18, 19, 20, 21];
(ii) stochastic and frictional forces accounting for the effect of collisions with
solvent molecules [22, 23, 24, 25, 2]; (iii) forces at the system boundary to mimick
a system-solvent interface [8, 9, 10, 11, 12, 13, 14, 15]. Spatial boundary
conditions are hard boundary conditions, because they apply strictly to all conﬁgurations
during a simulation.
2. Thermodynamical boundary conditions include the deﬁnition of the n + 2
thermodynamical quantities characterizing the macroscopic state of a (monoplastic)
n-component system (for systems under vacuum boundary conditions,
only n + 1 quantities are required because the volume is not deﬁned while
108 Philippe H. Hünenberger
the thermodynamical pressure is zero). These quantities can be selected form
pairs of extensive and intensive quantities including: (i) the number of particles
(N ≡ {Ni | i = 1...n}) or chemical potential (µµ ≡ {µi | i = 1...n}) of all
species; (ii) the volume V or pressure P; (iii) the energy E (or a related extensive
thermodynamical potential) or temperature T . The selected set of n + 2
quantities, together with their reference (macroscopic) values, deﬁne the thermodynamical
ensemble that is sampled during a simulation (Table 1). By default,
MD simulations sample microstates in the microcanonical (NV E) ensemble. By
applying speciﬁc modiﬁcations to the system Hamiltonian or equations of motion,
it is possible to maintain instead a constant temperature, pressure or chemical
potential for the different species (or any combination of these changes). The
thermodynamical boundary conditions involving extensive quantities should be
treated as hard boundary conditions, while those involving intensive quantities
should be soft.
3. Experimentally derived boundary conditions are used to explicitly enforce agreement
between a simulation and some experimental result. These may be applied
to enforce, e.g., the reproduction of (average) electron density maps from X-ray
crystallography [27, 28, 29, 30], or the agreement with (average) interatomic
distances and J-coupling constants from NMR measurements [31, 32, 33, 30].
Since experiments always provide averages over a given time and number of
molecules, experimentally derived boundary conditions should be handled as
soft boundary conditions.
4. Geometrical constraints can also be considered as boundary conditions. A typical
example is the use of bond-length constraints in simulations [34, 35, 36,
37, 38, 39], which represent a better approximation to the quantum-mechanical
behavior of high-frequency oscillators (hν kBT) compared to the classical
treatment [40]. Since they are satisﬁed exactly at every timepoint during a simulation,
geometrical constraints represent hard boundary conditions.
The present article is concerned with one speciﬁc type of thermodynamical
boundary condition, namely the imposition of a constant (average) temperature during
MD simulations by means of thermostat algorithms. The simultaneous enforcement
of a constant (average) pressure [51, 52, 53, 54, 55, 46, 56, 57, 58, 59, 60,
61, 53, 62, 63, 64, 65, 66, 67, 68, 69, 70] or chemical potential [71, 72, 73, 74]
will not be considered here. The discussion is also restricted to systems under either
vacuum or periodic boundary conditions, i.e., isolated systems. This implies
that the Hamiltonian is time-independent, and invariant upon translation or rotation
of the whole system. This Hamiltonian may contain terms accounting for the mean
effect of the environment (e.g., implicit-solvation term), as long as it still satisﬁes the
above conditions. The only exception considered here (whenever explicitly stated) is
the possible inclusion of stochastic and frictional forces as applied in stochastic dynamics
(SD) simulations, or of random collisional forces as applied in the stochasticcoupling
(Andersen) thermostat. Finally, it should be stressed that the inclusion of
geometrical constraints during a simulation affects the statistical mechanics of the
sampled microstates [75]. This is mainly because in the presence of such constraints,
Thermostat Algorithms 109
Table 1. The eight thermodynamical ensembles, and the corresponding independent and dependent
variables. Intensive variables are the chemical potential for all n species (µµ ≡ {µi |
i = 1...n}), the pressure (P) and the temperature (T ). Extensive variables are the number of
particles for all species (N ≡ {Ni | i = 1...n}), the volume (V ), and the energy (E), enthalpy
(H = E + PV ), Hill energy (L = E − µi Ni ), or Ray enthalpy (R = E + PV − µi Ni ).
Note that grand-ensembles may be open with respect to a subset of species only (e.g., semigrand-canonical
ensemble). The generalized ensemble is not a physical ensemble, because its
size is not speciﬁed (no independent extensive variable). Isothermal ensembles are discussed
in many standard textbooks. Speciﬁc references are given for the (less common) adiabatic
ensembles
Independent Dependent Ensemble
NV E µµPT Microcanonical [41, 42, 43, 44, 45]
NV T µµPE Canonical
NP H µµV T Isoenthalpic-isobaric [46, 47, 48, 45]
NPT µµV H Isothermal-isobaric (Gibbs)
µµV L NPT Grand-microcanonical [49, 45]
µµV T NPL Grand-canonical
µµP R NV T Grand-isothermal-isobaric [50, 45]
µµPT NV R Generalized
the kinetic energy of the system cannot be written in a conﬁguration-independent
way (unless the constraints are exclusively involved in fully-rigid atom groups, e.g.,
rigid molecules). This restriction limits the validity of a number of equations presented
in this article. However, many results are expected to remain approximately
valid for systems involving a small proportion of constrainted degrees of freedom,
and no attempt is made here to derive forms including explicitly the effect of geometrical
constraints.
2
Ensembles
An isolated system is characterized by a time-independent, translationally invariant
and rotationally invariant Hamiltonian. Integration of the classical equations of motion
for such a system leads, in the limit of inﬁnite sampling, to a trajectory mapping
a microcanonical (NV E) ensemble of microstates1. Assuming an inﬁnite numerical
precision, this is also what a standard MD simulation will deliver.
The laws of classical mechanics also lead to two additional conserved quantities,
namely the linear momentum psys of the system, and the angular momentum Lsys of
1 Thermodynamical ensembles are generally deﬁned without the constraint of Hamiltonian
translational and rotational invariance, in which case the previous statement is not entirely
correct. In the present article, however, the terminology of Table 1 will be (loosely) retained
to encompass ensembles where this invariance is enforced. The statistical mechanics of
these latter ensembles must be adapted accordingly [76, 77, 78, 79, 80, 81]. This requires
in particular the introduction of a modiﬁed deﬁnition for the instantaneous temperature,
relying solely on internal degrees of freedom and kinetic energy (Sect. 3).
110 Philippe H. Hünenberger
the system around its center of mass (CM). In simulations under periodic boundary
conditions, the two quantities refer to the inﬁnite periodic system. However, in this
case, if the linear momentum pbox of the computational box is also conserved, the
corresponding angular momentum Lbox is not. This is because correlated rotational
motion in two adjacent boxes exert friction on each other, leading to an exchange
of kinetic energy with the other (internal) degrees of freedom of the system. Note
that the physical properties of a molecular system are independent of psys. However,
they depend on Lsys, because the rotation of the system leads to centrifugal forces.
For this reason, Lsys should be added to the list of independent variables deﬁning the
ensemble sampled. Whenever Lsys is not given, it generally implicitly means that
Lsys = 0. The use of Lsys = 0 in simulations under periodic boundary conditions
(overall uniform rotation of the inﬁnite periodic system) is actually impossible, because
it would lead to non-periodic centrifugal forces. Finally, it should be speciﬁed
that the total energy E of the system is deﬁned here so as to exclude the kinetic energy
contributions corresponding to the overall translation and rotation of the system
(so that E is independent of psys and Lsys).
Because the independent variables of the microcanonical ensemble are all extensive,
they should be strictly conserved (i.e., time-independent) during the course of
a simulation. The corresponding dependent variables, namely the chemical potential
µµ, the pressure P, and the temperature T, are not conserved. In a non-equilibrium
simulation, these quantities may undergo a systematic drift. In an equilibrium simulation,
the corresponding instantaneous observables (denoted by νν, P, and T ) will
ﬂuctuate around well-deﬁned average values µµ, P, and T . Two important comments
should be made concerning the previous statement. First, the instantaneous observables
νν, P, and T are not uniquely deﬁned. The instantaneous temperature is generally
related to the total kinetic energy of the system (Eq. (8)), and the instantaneous
pressure to the total virial and kinetic energy. However, alternative deﬁnitions are
available (differing from the above by any quantity with a vanishing equilibrium
average), leading to identical average values in equilibrium situations, but to different
ﬂuctuations. Second, a microcanonical ensemble at equilibrium could equally
well be speciﬁed by stating that N, V , and E are conserved, and giving the values
of µµ instead of N, P instead of V, or T instead of E (as long as at least one
extensive variable is speciﬁed). However, such a speciﬁcation would be rather unnatural
as well as inapplicable to non-equilibrium situations. Furthermore, only the
natural variables for deﬁning a given thermodynamical ensemble (Table 1) are either
time-independent or characterized by vanishing ﬂuctuations in the limit of a
macroscopic system. Finally, it should be stressed that computer simulations cannot
be performed at inﬁnite numerical precision. As a consequence, quantities which
are formally time-independent in classical mechanics may still undergo a numerical
drift in simulations. In microcanonical simulations, this is typically the case for E,
as well as psys and Lsys (vacuum boundary conditions), or pbox (periodic boundary
conditions).
Unfortunately, the microcanonical ensemble that comes out of a standard MD
simulation does not correspond to the conditions under which most experiments are
Thermostat Algorithms 111
carried out. For comparison with experiment, the following ensembles are more useful,
which involve one or more intensive independent variables (Table 1):
1. In the canonical ensemble (NV T), the temperature has a speciﬁed average
(macroscopic) value, while the instantaneous observable representing the total
energy of the system (i.e., the Hamiltonian H) can ﬂuctuate. At equilibrium, the
root-mean-square ﬂuctuations σE of the Hamiltonian around its average value E
are related to the system isochoric heat capacity, cV , through [16]
σ2
E = H2
NV T
− H 2
NVT = kBT 2
cV . (1)
The ﬂuctuations σT of the instantaneous temperature T (deﬁned by Eq. (8)) in a
canonical ensemble are given by [16]
σ2
T = T 2
NV T
− T 2
NVT = 2N−1
d f T 2
, (2)
where Nd f is the number of internal degrees of freedom in the system (Eq. (9)).
These ﬂuctuations vanish in the limit of a macroscopic system, but are often
non-negligible for the system sizes typically considered in simulations.
2. In the isothermal-isobaric (Gibbs) ensemble (NPT ), the pressure has (just as
the temperature) a speciﬁed average value, while the instantaneous volume V of
the system can ﬂuctuate. At equilibrium, the root-mean-square ﬂuctuations σV
of the instantaneous volume around its average value V are related to the system
isothermal compressibility, βT , through [16]
σ2
V = V2
NPT − V 2
NPT = VkBTβT . (3)
The root-mean-square ﬂuctuations σH of the instantaneous enthalpy H + PV
around its average value H are related to the system isobaric heat capacity, cP,
through [16]
σ2
H = (H + PV)2
NPT − H + PV 2
NPT = kBT2
cP . (4)
Both the instantaneous temperature T and the instantaneous pressure P will
ﬂuctuate around their corresponding macroscopic values, the magnitude of these
ﬂuctuations vanishing in the limit of a macroscopic system.
3. The grand-canonical ensemble (µµV T ) has a constant volume and temperature
(as the canonical ensemble), but is open for exchanging particles with a surrounding
bath. In this case, the chemical potential of the different species has
a speciﬁed average, while the instantaneous value NN of the number of particles
can ﬂuctuate. For a one-component system at equilibrium, the ﬂuctuations σN
of the instantaneous number of particles around its average value N are related
to the system isothermal compressibility, βT , through [16]
σ2
N = N2
µV T
− N 2
µV T = N2
V−1
kBTβT . (5)
112 Philippe H. Hünenberger
The root-mean-square ﬂuctuations σL of the instantaneous Hill energy H − µN
around its average value L are given by [16]
σ2
L = (H − µN)2
µV T
− H − µN 2
µV T =kBT2 ∂L(µ, V, T )
∂T µV
. (6)
Three other combinations of variables are possible (Table 1), but the corresponding
ensembles [45] are of more limited practical relevance. The last combination
(generalized ensemble) is not physical, because its size is not speciﬁed (no independent
extensive variable). Note that although MD samples the microcanonical ensemble
by default, the basic Monte Carlo (MC; [82, 83, 84]) and stochastic dynamics
(SD; [22, 23, 24, 25, 2]) algorithms sample the canonical ensemble.
Performing a MD simulation in an other ensemble than microcanonical requires a
means to keep at least one intensive quantity constant (on average) during the simulation.
This can be done either in a hard or in a soft manner. Applying a hard boundary
condition on an intensive macroscopic variable means constraining a corresponding
instantaneous observable to its speciﬁed macroscopic value at every timepoint during
the simulation (constraint method). Remember, however, that the choice of this
instantaneous observable is not unique. In contrast, the use of a soft boundary condition
allows for ﬂuctuations in the instantaneous observable, only requiring its average
to remain equal to the macroscopic value (on a given timescale). Typical methods
for applying soft boundary conditions are the penalty-function, weak-coupling,
extended-system and stochastic-coupling methods [85]. These methods will be discussed
in the following sections in the context of constant-temperature simulations.
Although there are many ways to ensure that the average of an instantaneous quantity
takes a speciﬁed value, ensuring that the simulation actually samples the correct
ensemble (and in particular provides the correct ﬂuctuations for the speciﬁc instantaneous
observable in the given ensemble) is much more difﬁcult.
3
Thermostat Algorithms
A modiﬁcation of the Newtonian MD scheme with the purpose of generating a thermodynamical
ensemble at constant temperature is called a thermostat algorithm. The
use of a thermostat can be motivated by one (or a number) of the following reasons:
(i) to match experimental conditions (most condensed-phase experiments are performed
on thermostatized rather than isolated systems); (ii) to study temperaturedependent
processes (e.g., determination of thermal coefﬁcients, investigation of
temperature-dependentconformationalor phase transitions); (iii) to evacuate the heat
in dissipative non-equilibrium MD simulations (e.g., computation of transport coefﬁcients
by viscous-ﬂow or heat-ﬂow simulations); (iv) to enhance the efﬁciency of
a conformational search (e.g., high-temperature dynamics, simulated annealing); (v)
Thermostat Algorithms 113
to avoid steady energy drifts caused by the accumulation of numerical errors during
MD simulations2.
The use of a thermostat requires the deﬁnition of an instantaneous temperature.
This temperature will be compared to the reference temperature To of the heat bath to
which the system is coupled. Following from the equipartition theorem, the average
internal kinetic energy K of a system is related to its macroscopic temperature T
through
K = K =
1
2
kB Nd f T (7)
where kB is Boltzmann’s constant, Nd f the number of internal degrees of freedom
of the system, and K its instantaneous internal kinetic energy. Deﬁning the instantaneous
temperature T at any timepoint as
T =
2
kB Nd f
K , (8)
one ensures that the average temperature T is identical to the macroscopic temperature
T . This deﬁnition is commonly adopted, but by no means unique. For example,
the instantaneous temperature could be deﬁned based on the equipartition principle
for only a subset of the internal degrees of freedom. It may also be deﬁned purely on
the basis of conﬁguration, without any reference to the kinetic energy [87, 88].
In the absence of stochastic and frictional forces (see below; Eq. (17)), a few
degrees of freedom are not coupled (i.e., do not exchange kinetic energy) with the
internal degrees of freedom of the system. These external degrees of freedom correspond
to the system rigid-body translation and, under vacuum boundary conditions,
rigid-body rotation. Because the kinetic energy associated with these external degrees
of freedom can take an arbitrary (constant) value determined by the initial
atomic velocities, they must be removed from the deﬁnition of the system internal
temperature. Consequently, the number of internal degrees of freedom is calculated
as three times the total number N of atoms in the system, minus the number Nc of
geometrical constraints, i.e.
Nd f = 3N − Nc − Nr . (9)
The subtraction of constrained degrees of freedom is necessary because geometrical
constraints are characterized by a time-independent generalized coordinate associated
with a vanishing generalized momentum (i.e., no kinetic energy). A more formal
statistical-mechanical justiﬁcation for the subtraction of the external degrees of
2 A thermostat algorithm (involving explicit reference to a heat-bath temperature To) will
avoid systematic energy drifts, because if the instantaneous temperature is forced to ﬂuctuate
within a limited range around To, the energy will also ﬂuctuate within a limited range
around its corresponding equilibrium value. To perform long microcanonical simulations
(no thermostat), it is also advisable to employ an algorithm that will constrain the energy
to its reference value Eo (ergostat algorithm [86]).
114 Philippe H. Hünenberger
freedom in the case of periodic boundary conditions can be found elsewhere [45].
A corresponding derivation for vacuum boundary conditions has, to our knowledge,
never been reported. When stochastic and frictional forces are applied, as in SD,
these forces will couple the rigid-body translational and rotational degrees of freedom
with the internal ones. In this case all degrees of freedom are considered internal
to the system. Thus, Eq. (9) is to be used with Nr = 0 in the presence of stochastic
and frictional forces, and otherwise with Nr = 3 under periodic boundary conditions
or Nr = 6 under vacuum boundary conditions. Similarly, the instantaneous internal
kinetic energy is deﬁned as
K =
1
2
N
i=1
mi ˙r2
i , (10)
where the internal (also called peculiar) velocities ˙ri are obtained from the real
atomic velocities ˙ro
i by excluding any component along the external degrees of freedom3.
These corrected velocities are calculated as
˙ri =
⎧
⎨
⎩
˙ro
i if Nr = 0
˙ro
i − ˙ro
C M if Nr = 3
˙ro
i − ˙ro
C M − I−1
C M(ro) Lo
C M × (ro
i − ro
C M) if Nr = 6
, (11)
where ro
C M is the coordinate vector of the system center of mass (CM), Lo
C M the
system angular momentum about the CM, and IC M is the (conﬁguration-dependent)
inertia tensor of the system relative to the CM. The latter quantity is deﬁned as
IC M(r) =
N
i=1
mi (ri − rC M) ⊗ (ri − rC M) , (12)
where a ⊗ b denotes the tensor with elements µ, ν equal to aµbν. Application of
Eq. (11) ensures that
N
i=1
mi ˙ri = 0 for Nr = 3 or 6 (13)
and (irrespective of the origin of the coordinate system)
N
i=1
mi ri × ˙ri = 0 for Nr = 6. (14)
Equation (13) is a straightforward consequence of the deﬁnition of ro
C M. Equation
(14) is proved by using ωωo
C M × (ro
i − ro
C M) = ˙ro
i − ˙ro
C M where ωωo
C M =
I−1
C M(ro) Lo
C M is the angular velocity vector about the CM. Thus, the linear and angular
momenta of the internal velocities vanish, as expected.
3 It is assumed that the velocities ˙ro
i are already exempt of any component along possible
geometrical constraints.
Thermostat Algorithms 115
Because the instantaneous temperature is directly related to the atomic internal
velocities (Eqs. (8) and (10)), maintaining the temperature constant (on average)
in MD simulations requires imposing some control on the rate of change of these
velocities. For this reason, thermostat algorithms require a modiﬁcation of Newton’s
second law4
¨ri(t) = m−1
i Fi(t) . (15)
In the present context, this equation (and the thermostatized analogs discussed below)
should be viewed as providing the time-derivative of the internal velocity ˙ri
deﬁned by Eq. (11). In turn, ˙ri is related to the real atomic velocity ˙ro
i through the
inverse of Eq. (11), namely
˙ro
i =
⎧
⎨
⎩
˙ri if Nr = 0
˙ri + ˙r∗
C M if Nr = 3
˙ri + ˙r∗
C M + I−1
C M(ro) L∗
C M × (ro
i − ro
C M) if Nr = 6
, (16)
where r∗
C M and L∗
C M are constant parameters determined by the initial velocities
˙ro
i (0). This distinction between real and internal velocities is often ignored in standard
simulation programs. Many programs completely disregard the problem, while
others only remove the velocity component along the external degrees of freedom for
the computation of the temperature (but do not use internal velocities in the equations
of motion). However, as discussed in Sect. 4, this can have very unpleasant consequences
in practice. In the following discussion, it is assumed that the equation of
motion (Eq. (15) or any thermostatized modiﬁcation) is applied to the internal velocities
deﬁned by Eq. (11), while the atomic coordinates are propagated simultaneously
in time using the real velocities ˙ro
i deﬁned by Eq. (16).
The prototype of most isothermal equations of motion is the Langevin equation
(as used in SD; see Sect. 3.2), i.e.
¨ri(t) = m−1
i Fi(t) − γi(t)˙ri (t) + m−1
i Ri(t) , (17)
where Ri is a stochastic force and γi a (positive) atomic friction coefﬁcient. Many
thermostats avoid the stochastic force in Eq. (17) and use a single friction coefﬁcient
for all atoms. This leads to the simpliﬁed form
¨ri(t) = m−1
i Fi(t) − γ (t)˙ri (t) . (18)
In this case, γ loses its physical meaning of a friction coefﬁcient and is no longer
restricted to positive values. A positive value indicates that heat ﬂows from the system
to the heat bath. A negative value indicates a heat ﬂow in the opposite direction.
Note that if Eq. (18) was applied to the real velocities ˙ro
i (as often done in simulation
programs) instead of the internal velocities ˙ri, the linear and angular momenta of the
system would not be conserved (unless they exactly vanish).
4 It is assumed that the forces Fi are exempt of any component along possible geometrical
constraints.
116 Philippe H. Hünenberger
Any algorithm relying on the equation of motion given by Eq. (18) is smooth
(i.e., generates a continuous velocity trajectory) and deterministic5. It is also timereversible
if γ is antisymmetric with respect to time-reversal6.
Practical implementations of Eq. (18) often rely on the stepwise integration of
Newton’s second law (Eq. (15)), altered by the scaling of the atomic velocities after
each iteration step. In the context of the leap-frog integrator7 [89], this can be written
˙ri t +
t
2
= λ(t; t) ˙ri t +
t
2
= λ(t; t) ˙ri(t −
t
2
) + m−1
i Fi(t) t , (20)
where λ(t; t) is a time- and timestep-dependent velocity scaling factor. Imposing
the constraint8 λ(t; 0) = 1, one recovers Eq. (18) in the limit of an inﬁnitesimal
timestep t, with
γ (t) = − lim
t→0
λ(t; t) − 1
t
= −
∂λ(t; t)
∂( t)
| t=0
. (21)
Note that for a given equation of motion, i.e., a speciﬁed form of γ (t), Eq. (21) does
not uniquely specify the scaling factor λ(t; t). It can be shown that Eq. (20) retains
the original accuracy of the leap-frog algorithm if the velocity-scaling factor applied
to atom i is chosen as [90]
λi(t; t) = 1 − γ (t) t +
γ 2(t)
2
+
γ (t)Fi (t)
2mi ˙ri(t)
( t)2
. (22)
From a thermodynamical point of view, some thermostats can be proved to generate
(at constant volume and number of atoms) a canonical ensemble in the limit of
inﬁnite sampling times (and within the usual statistical-mechanical assumptions of
equal a priori probabilities and ergodicity). More precisely, some thermostats lead to
a canonical ensemble of microstates, i.e., microstates are sampled with a statistical
5 The advantages of deterministic algorithms are that (i) the results can be exactly reproduced
(in the absence of numerical errors), and (ii) there are well-deﬁned conserved
quantities (constants of the motion). In the case of Eq. (18), the constant of the motion is
C = K(t) + U(t) + 2
t
0
dt K(t)γ (t) . (19)
6 Considering a given microstate, time-reversibility is achieved if the change dt → −dt
(leading in particular to r → r, ˙r → −˙r, and ¨r → ¨r) leaves the equation of motion for the
coordinates unaltered (while the velocities are reversed). Clearly, this condition is satisﬁed
for Eq. (18) only if the corresponding change for γ is γ → −γ .
7 The implementation of thermostats will only be discussed here in the context of the leapfrog
integrator. However, implementation with other integrators is generally straightfor-
ward.
8 An algorithm with λ(t; 0) = 1 would involve a Dirac delta function in its equation of
motion.
Thermostat Algorithms 117
weight proportional to e−βH where β = (kBTo)−1. In this case and in the absence
of geometrical constraints, expressing the Hamiltonian in Cartesian coordinates as
H(r, p) = U(r) + K(p), U being the potential energy, the probability distribution
of microstates may be written
ρ(r, p) =
e−βH(r,p)
dr dp e−βH(r,p)
=
e−βU(r)
dr e−βU(r)
e−βK(p)
dp e−βK(p)
. (23)
Integrating this expression over either momenta or coordinates shows that the distribution
is also canonical in both conﬁgurations (i.e., conﬁgurations are sampled
with a statistical weight proportional to e−βU) and momenta (i.e., momenta are sampled
with a statistical weight proportional to e−βK ). In Cartesian coordinates, such a
canonical distribution of momenta reads
ρp(p) =
e−βK(p)
dp e−βK(p)
=
3N
iµ
e−β(2mi)−1 p2
iµ
dpiµ e−β(2mi)−1 p2
iµ
=
3N
iµ
p(piµ) , (24)
where Eq. (10) was used together with pi = mi ˙ri. Noting that p(˙riµ) = mi p(piµ)
and evaluating the required Gaussian integral, this result shows that internal velocities
obey a Maxwell-Boltzmann distribution, i.e., the velocity components ˙riµ appear
with the probability
p(˙riµ) =
βmi
2π
1/2
e−(1/2)βmi ˙r2
iµ . (25)
Note that the above statements do not formally hold in the presence of geometrical
constraints, but are generally assumed to provide a good approximation in this case.
Some other thermostats only generate a canonical ensemble of conﬁgurations, but
not of microstates and momenta. This is generally not a serious disadvantage for the
computation of thermodynamical properties, because the contribution of momenta to
thermodynamical quantities can be calculated analytically (ideal-gas contribution).
Finally, there also exists thermostats that generate distributions that are canonical
neither in conﬁgurations nor in momenta.
From a dynamical point of view, assessing the relative merits of different thermostats
is somewhat subjective9. Clearly, the conﬁgurational dynamics of a system
will be affected by the timescale of its instantaneous temperature ﬂuctuations, and a
good thermostat should reproduce this timescale at least qualitatively. However, the
direct comparison between experimental thermostats (e.g., a heat bath surrounding
9 An objective question, however, is whether the thermostat is able to produce correct timecorrelation
functions (at least in the limit of a macroscopic system). Since transport coefﬁcients
(e.g., the diffusion constant) can be calculated either as ensemble averages (Einstein
formulation) or as integrals of a time-correlation function (Green-Kubo formulation), at
least such integrals should be correct if the thermostat leads to a canonical ensemble. When
this is the case, it has been shown that the correlation functions themselves are also correct
at least for some thermostats [91, 86].
118 Philippe H. Hünenberger
a macroscopic system, or the bulk medium around a microscopic sample of matter)
and thermostats used in simulations is not straightforward. The reason is that experimental
thermostats, because they involve the progressive diffusion of heat from the
system surface towards its center (or inversely), lead to inhomogeneities in the spatial
temperature distribution within the sample. On the other hand, the thermostats
used in simulations generally modify instantaneously and simultaneously the velocities
of all atoms irrespective of their locations, and should lead to an essentially
homogeneous temperature distribution.
One may nevertheless try to quantify the timescale of the temperature ﬂuctuations
to be expected in a thermostatized simulation. This timescale can be estimated based
on a semi-macroscopic approach [55]. Consider a system characterized by an average
temperature T , in contact with a heat bath at a different temperature To. By average
temperature, it is meant that the quantity T is spacially-averaged over the entire
system and time-averaged over an interval that is short compared to the experimental
timescale, but long compared to the time separating atomic collisions. The difference
between T and To may result, e.g., from a natural ﬂuctuation of T within the system.
From macroscopic principles, the rate of heat transfer from the heat bath to the the
system should be proportional to the temperature difference To−T and to the thermal
conductivity κ of the system. Thus, the rate of change in the average temperature can
be written (at constant volume)
˙T (t) = c−1
v
˙E(t) = ζ−1
T [To − T (t)] (26)
with the deﬁnition
ζT = ξ−1
V−1/3
cvκ−1
, (27)
where cv is the system isochoric heat capacity, V the system volume, and ξ a dimensionless
constant depending on the system shape and on the temperature inhomogeneity
within the system. For a given system geometry (e.g., spherical) and initial
temperature distribution (i.e., T (x, 0)), a reasonable value for ξ could in principle be
estimated by solving simultaneously the ﬂux equation
J(x, t) = −κ∇T(x, t), (28)
where J(x, t) is the energy ﬂux through a surface element perpendicular to the direction
of the vector, and the conservation equation
∂T(x, t)
∂t
= −Vc−1
v ∇ · J(x, t) . (29)
Eq. (26) implies that, at equilibrium, the natural ﬂuctuations of T away from To
decay exponentially with a temperature relaxation time ζT , i.e.
T (t) = To + [T(0) − To] e−ζ−1
T t
. (30)
Note that on a very short timescale (i.e., of the order of the time separating atomic
collisions), the instantaneous temperature T (t) is also affected by important stochastic
variations (see Sect. 3.2). Only on an intermediate timescale does the mean effect
Thermostat Algorithms 119
of these stochastic ﬂuctuations result in an exponential relaxation for T (t). However,
because stochastic variations contribute signiﬁcantly to the instantaneous temperature
ﬂuctuations, a thermostat based solely on an exponential relaxation for T (t)
leads to incorrect (underestimated) temperature ﬂuctuations (see Sect. 3.5).
To summarize, although assessing whether one thermostat leads to a better description
of the dynamics compared to another one is largely subjective, it seems
reasonable to assume that: (i) thermostats permitting temperature ﬂuctuations are
more likely to represent the dynamics correctly compared to thermostats constraining
the temperature at a ﬁxed value; (ii) thermostats with temperature ﬂuctuations
are more likely to represent the dynamics correctly when these ﬂuctuations occur
at a timescale (measured in a simulation, e.g., as the decay time of the temperature
autocorrelation function) of the order of ζT (Eq. (26)), and when the dynamics is
smooth (continuous velocity trajectory). These differences will be more signiﬁcant
for small systems, where the temperature ﬂuctuations are of of larger magnitudes (the
corresponding root-mean-square ﬂuctuations scale as N−1/2, see Eq. (2)) and higher
frequencies (the corresponding relaxation times scale as N−1/3, see Eq. (27)).
A summary of the common thermostats used in MD simulations, together with
their main properties, is given in Table 2. The various algorithms are detailed in the
following sections.
3.1
Temperature in the Monte Carlo Algorithm
Although the present discussion mainly focuses on thermostatized MD, the simplest
way to generate a thermodynamical ensemble at constant temperature is to use the
MC algorithm [82, 83, 84]. This algorithm does not involve atomic velocities or
kinetic energy. Random trial moves are generated, and accepted with a probability
p = min{e−β U
, 1} (31)
depending on the potential energy change U associated with the move and on
the reference temperature To. Following this criterion, moves involving rigid-body
translation and, under vacuum boundary conditions, rigid-body rotation are always
accepted because they do not change the potential energy. For this reason, the corresponding
degrees of freedom are external to the system. Note also that under vacuum
boundary conditions, the centrifugal forces due to the rigid-body rotation of the system,
which would be included in a MD simulation, are absent in the MC procedure.
Therefore, MC samples by default an ensemble at zero angular momentum. It can be
shown that the ensemble generated by the MC procedure represents (at constant volume)
a canonical distribution of conﬁgurations. The modiﬁcation of the MC scheme
to sample other isothermal ensembles (including the grand-canonical ensemble [92])
is possible. Modiﬁcations permitting the sampling of adiabatic ensembles (e.g., the
microcanonical ensemble [92, 93, 94]) have also been devised. The MC procedure is
non-smooth, non-deterministic, time-irreversible, and does not provide any dynamical
information.
120 Philippe H. Hünenberger
Table 2. Characteristics of the main thermostat algorithms used in MD simulations. MD:
molecular dynamics (generates a microcanonical ensemble, only shown for comparison); MC:
Monte Carlo (Sect. 3.1); SD: stochastic dynamics (with γi > 0 for at least one atom; Sect. 3.2);
A: MD with Andersen thermostat (with α > 0; Sect. 3.3); HE: MD with Hoover-Evans thermostat
(Sect. 3.4); W: MD with Woodcock thermostat (Sect. 3.4); HG: MD with Haile-Gupta
thermostat (Sect. 3.4); B: MD with Berendsen thermostat (with t < τB < ∞; Sect. 3.5);
NH: MD with Nosé-Hoover thermostat (with 0 < Q < ∞; Sect. 3.6). MD is a limiting case
of SD (with γi = 0 for all atoms), A (with α = 0), B (with τB → ∞), and NH (with Q → ∞,
γ (0) = 0). HE/W is a limiting cases of B (with τB = t) and is a constrained form of NH. HG
is also a constrained form of NH. Deterministic: trajectory is deterministic; Time-reversible:
equation of motion is time-reversible; Smooth: velocity trajectory is available and continuous.
Energy drift: possible energy (and temperature) drift due to accumulation of numerical errors;
Oscillations: possible oscillatory behavior of the temperature dynamics; External d.o.f.: some
external degrees of freedom (rigid-body translation and, under vacuum boundary conditions,
rotation) are not coupled with the internal degrees of freedom. Constrained K: no kinetic energy
ﬂuctuations; Canonical in H: generates a canonical distribution of microstates; Canonical
in U: generates a canonical distribution of conﬁgurations. Dynamics: dynamical information
on the system is either absent (−−) or likely to be unrealistic (−; constrained temperature
or non-smooth trajectory), moderately realistic (+; smooth trajectory, but temperature ﬂuctuations
of incorrect magnitude), or realistic (++; smooth trajectory, correct magnitude of
the temperature ﬂuctuations). The latter appreciation is rather subjective and depends on an
adequate choice of the adjustable parameters of the thermostat
MD MC SD A HE W HG B NH
Deterministic + − − − + + + + +
Time−reversible + − − − + + + − +
Smooth + − + − + + + + +
Energy drift + − − − + − − − −
Oscillations − − − − − − − − +
External d.o.f. + + − − + + + + +
Constrained K − − − − + + + − −
Canonical in H − − + + − − − − +
Canonical in U − + + + + + − − +
Dynamics ++ − − ++ − − − − + ++
Eqn. of motion 15 17 41 46 51 52 57 78,79
3.2
Temperature Relaxation by Stochastic Dynamics
The SD algorithm relies on the integration of the Langevin equation of motion [22,
23, 95, 96, 97, 98, 99, 100] as given by Eq. (17). The stochastic forces Ri (t) have the
following properties10: (i) they are uncorrelated with the velocities ˙r(t ) and systematic
forces Fi(t ) at previous times t < t; (ii) their time-averages are zero; (iii) their
mean-square components evaluate to 2miγikBTo; (iv) the force component Riµ(t)
10 More complex SD schemes can be used, which incorporate time or space correlations
in the stochastic forces. It is also assumed here that the friction coefﬁcients γi are time-
independent.
Thermostat Algorithms 121
along the Cartesian axis µ is uncorrelated with any component Rjν(t ) along axis ν
unless i = j, µ = ν, and t = t. The two last conditions can be combined into the
relation
Riµ(t)Rjν(t ) = 2miγikBToδij δµνδ(t − t) . (32)
It can be shown that a trajectory generated by integrating the Langevin equation
of motion (with at least one non-vanishing atomic friction coefﬁcient γi) maps (at
constant volume) a canonical distribution of microstates at temperature To.
The Langevin equation of motion is smooth, non-deterministic and timeirreversible.
Under vacuum boundary conditions and aiming at reproducing bulk
properties, it may produce a reasonable picture of the dynamics if the mean effect
of the surrounding solvent is incorporated into the systematic forces, and if the friction
coefﬁcients are representative of the solvent viscosity (possibly weighted by
the solvent accessibility). If SD is merely used as a thermostat in explicit-solvent
simulations, as is the case, e.g., when applying stochastic boundary conditions to a
simulated system [101, 11, 12], some care must be taken in the choice of the atomic
friction coefﬁcients γi. On the one hand, too small values (loose coupling) may cause
a poor temperature control. Indeed, the limiting case of SD where all friction coefﬁcients
(and thus the stochastic forces) are set to zero is MD, which generates a
microcanonical ensemble. However, arbitrarily small atomic friction coefﬁcients (or
even a non-vanishing coefﬁcient for a single atom) are sufﬁcient to guarantee in
principle the generation of a canonical ensemble. But if the friction coefﬁcients are
chosen too low, the canonical distribution will only be obtained after very long simulation
times. In this case, systematic energy drifts due to accumulation of numerical
errors may interfere with the thermostatization. On the other hand, too large values
of the friction coefﬁcients (tight coupling) may cause the large stochastic and frictional
forces to perturb the dynamics of the system. In principle, the perturbation of
the dynamics due to stochastic forces will be minimal when the atomic friction coefﬁcients
γi are made proportional to mi . In this case, Eqs. (17) and (32) show that
the root-mean-square acceleration due to stochastic forces is identical for all atoms.
In practice, however, it is often more convenient to set the friction coefﬁcients to a
common value γ . The limiting case of SD for very large friction coefﬁcients (i.e.,
when the acceleration ¨ri can be neglected compared to the other terms in Eq. (17))
is Brownian dynamics (BD), with the equation of motion
˙ri(t) = γ −1
i m−1
i [Fi(t) + Ri(t)] . (33)
Although the magnitude of the temperature ﬂuctuations is in principle not affected
by the values of the friction coefﬁcients (unless they are all zero), the timescale
of these ﬂuctuations strongly depends on the γi coefﬁcients. In fact, it can be
shown [54] that there is a close relationship between the friction coefﬁcients in SD
(used as a mere thermostat) and the temperature relaxation time ζT in Eq. (26). Consider
the case where all coefﬁcients γi are set to a common value γ . Following from
Eqs. (8) and (10), the change T of the instantaneous temperature over a time interval
from t = 0 to τ can be written
122 Philippe H. Hünenberger
T =
2
kB Nd f
τ
0
dt ˙K(t) =
2
kB Nd f
N
i=1
mi
τ
0
dt ¨ri(t) · ˙ri(t) . (34)
Inserting Eq. (17), this can be rewritten
T =
2
kB Nd f
N
i=1
τ
0
dt Fi(t) − γ mi ˙ri(t) · ˙ri(t) +
τ
0
dt (35)
× Ri(t) · ˙ri(0) +
t
0
dt m−1
i Fi(t ) − γ ˙ri (t ) + m−1
i Ri(t ) .
Using Eq. (32) and the fact that the stochastic force is uncorrelated with the velocities
and systematic forces at previous times, this simpliﬁes to
T =
2
kB Nd f
N
i=1
{
τ
0
dt [Fi(t) · ˙ri(t) − γ mi ˙r2
i (t)]
+ m−1
i
τ
0
dt Ri(t) ·
t
0
dt Ri(t )}
=
2
kB Nd f
N
i=1
τ
0
dt [Fi(t) · ˙ri(t) − γ mi ˙r2
i (t)]
+ 6N−1
d f Nγ To τ . (36)
This expression can be rewritten11
T
τ
=
2
kB Nd f
N
i=1
Fi · ˙ri + 2γ [To − T ] , (37)
where Fi · ˙ri and T stand for averages over the interval τ. The ﬁrst term represents
the temperature change caused by the effect of the systematic forces, and would be
unaltered in the absence of thermostat (Newtonian MD simulation). Thus, the second
term can be identiﬁed with a temperature change arising from the coupling to a heat
bath. This means that on an intermediate timescale (as deﬁned at the end of Sect. 3),
the mean effect of thermostatization can be written
˙T (t) = 2γ [To − T (t)] . (38)
Comparing with Eq. (26) allows to identify 2γ with the inverse of the temperature
relaxation time ζT in Eq. (26), i.e., to suggest γ = (1/2)ζ−1
T as an appropriate value
for simulations. This discussion also shows that the semi-macroscopic expression
of Eq. (26) is only valid on an intermediate timescale, when the stochastic ﬂuctuations
occuring on a shorter timescale (i.e., of the order of the time separating atomic
collisions) are averaged out and only their mean effect is retained.
11 In the absence of constraints Nd f = 3N due to Eq. (9) with Nc = 0 and Nr = 0 (as
appropriate for SD). In the presence of constraints, the derivation should include constraint
forces.
Thermostat Algorithms 123
3.3
Temperature Relaxation by Stochastic Coupling
The stochastic-coupling method was proposed by Andersen [55]. In this approach,
Newton’s equation of motion (Eq. (15)) is integrated in time, with the modiﬁcation
that at each timestep, the velocity of all atoms are conditionally reassigned from a
Maxwell-Boltzmann distribution. More precisely, if an atom i is selected for a velocity
reassignment, each Cartesian component µ of the new velocity is selected at
random according to the Maxwell-Boltzmann probability distribution of Eq. (25).
The selection procedure is such that the time intervals τ between two successive
velocity reassignments of a given atom are selected at random according to a probability
p(τ) = αe−ατ , where α is a constant reassignment frequency. In principle, one
can select at random and for each atom a series of successive τ-values (obeying the
speciﬁed probability distribution) before starting the simulation. This series is then
used to determine when the particle is to undergo a velocity reassignment. In practice,
a simpler procedure can be used when t α−1 (infrequent reassignments).
At each timestep and for each atom in turn, one generates a random number between
0 and 1. If this number is larger than α t for a given atom, this atom undergoes a
velocity reassignment. This procedure leads to a probability distribution
p(τ) t = (1 − α t)τ/ t
α t (39)
for the intervals τ without velocity reassignment. This implies
ln α−1
p(τ) = (τ/ t) ln(1 − α t) = −ατ + O[(α t)2
] . (40)
Thus, when t α−1, p(τ) = αe−ατ , as expected. If the condition is not satisﬁed,
this second method will not work because the probability of multiple reassignments
within the same timestep becomes non-negligible.
The equation of motion for the Andersen thermostat can formally be written
¨ri(t) = m−1
i Fi(t) +
∞
n=1
δ t −
n
m=1
τi,m ˙r∗
i,n(t) − ˙ri(t) , (41)
where {τi,n | n = 1, 2, . . .} is the series of intervals without reassignment for particle
i, and ˙r∗
i,n the randomly-reassigned velocity after the nth interval. This approach
mimicks the effect irregularly-occurring stochastic collisions of randomly chosen
atoms with a bath of ﬁctitious particles at a temperature To. Because, the system
evolves at constant energy between the collisions, this method generates a succession
of microcanonical simulations, interrupted by small energy jumps corresponding to
each collision.
It can be shown [55] that the Andersen thermostat with non-zero collision frequency
α leads to a canonical distribution of microstates. The proof [55] involves
similar arguments as the derivation of the probability distribution generated by the
MC procedure. It is based on the fact that the Andersen algorithm generates a Markov
chain of microstates in phase space. The only required assumption is that every microstate
is accessible from every other one within a ﬁnite time (ergodicity). Note
124 Philippe H. Hünenberger
also that the system total linear and angular momenta are affected by the velocity
reassignments, so that these degrees of freedom are internal to the system, as in SD.
The Andersen algorithm is non-deterministic and time-irreversible. Moreover, it
has the disadvantage of being non-smooth, i.e., generating a discontinuous velocity
trajectory where the randomly-occurring collisions may interfere with the natural
dynamics of the system.
Some care must be taken in the choice of the collision frequency α [55, 102]. On
the one hand, too small values (loose coupling) may cause a poor temperature control.
The same observations apply here as those made for SD. The limiting case of
the Andersen thermostat with a vanishing collision frequency is MD, which generates
a microcanonical ensemble. Arbitrarily small collision frequencies are sufﬁcient
to guarantee in principle the generation of a canonical ensemble. But if the collision
frequency is too low, the canonical distribution will only be obtained after very long
simulation times. In this case, systematic energy drifts due to accumulation of numerical
errors may interfere with the thermostatization. On the other hand, too large
values for the collision frequency (tight coupling) may cause the velocity reassignments
to perturb heavily the dynamics of the system. Although the magnitude of
the temperature ﬂuctuations is in principle not affected by the value of the collision
frequency (unless it is zero), the timescale of these ﬂuctuations strongly depends on
this parameter. In fact, it can be shown [55] that there is a close relationship between
the collision frequency and the temperature relaxation time ζT in Eq. (26). Each collision
changes the kinetic energy of the system by (3/2)kB[To − T (t)] on average,
and there are Nα such collisions per unit of time. Thus, one expects
˙T (t) = c−1
v
˙E(t) = (3/2)c−1
v NαkB[To − T (t)] , (42)
where T and E stand for averages over an intermediate timescale (as deﬁned at the
end of Sect. 3). Comparing with Eq. (26) allows to identify (3/2)c−1
v NαkB with
the inverse of the temperature relaxation time ζT in Eq. (26), i.e., to suggest α =
(2/3)(NkB)−1cvζ−1
T as an appropriate value for simulations. Note that because ζT
scales as N−1/3 (Eq. (27)), the collision frequency for any particle scales as N−2/3,
so that the time each particle spends without reassignment increases with the system
size. On the other hand, the collision frequency for the whole system, Nα, scales as
N1/3, so that the length of each microcanonical sub-simulation decreases with the
system size.
3.4
Temperature Constraining
Temperature constraining aims at ﬁxing the instantaneous temperature T to the reference
heat-bath value To without allowing for any ﬂuctuations. In this sense, temperature
constraining represents a hard boundary condition, in constrast to the soft
boundary conditions employed by all other thermostats mentioned in this article.
Note that constraining the temperature, i.e., enforcing the relation T (t) = To (or
˙T (t) = 0) represents a non-holonomic constraint. Holonomic constraints are those
Thermostat Algorithms 125
which only involve generalized coordinates and time, excluding any dependence on
the generalized velocities. Two main temperature-constraining algorithms have been
proposed. The ﬁrst one is due to Woodcock [103], and the second one was simultaneously
proposed by Hoover [104] and Evans [105].
In the Hoover-Evans algorithm [104, 52, 51, 105, 106], the quantity λ(t; t) in
Eq. (20) is found by imposing temperature conservation in the form T (t + t
2 ) =
T (t − t
2 ). Using Eqs. (8) and (10), this leads to the condition
λ2(t; t)
kB Nd f
N
i=1
mi ˙ri t −
t
2
+ m−1
i Fi(t) t
2
=
1
kB Nd f
N
i=1
mi ˙r2
i t −
t
2
. (43)
Solving for λ(t; t) gives
λ(t; t) =
⎧
⎪⎨
⎪⎩
N
i=1 mi ˙r2
i (t − t
2 )
N
i=1 mi ˙ri(t − t
2 ) + m−1
i Fi(t) t
2
⎫
⎪⎬
⎪⎭
1/2
=
T (t − t
2 )
T (t + t
2 )
1/2
, (44)
where T (t − t
2 ) and T (t + t
2 ) are the instantaneous temperatures computed based
on the velocities ˙ri(t− t
2 ) and ˙ri(t+ t
2 ), see Eq. (20). Because this quantity satisﬁes
λ(t; 0) = 1, applying Eq. (21) gives
γ (t) = [Nd f kBT (t)]−1
N
i=1
˙ri (t) · Fi(t) . (45)
Inserting into Eq. (18) shows that the equation of motion corresponding to the
Hoover-Evans thermostat is
¨ri(t) = m−1
i Fi(t) − Nd f kBT (t)
−1
N
i=1
˙ri (t) · Fi(t) ˙ri(t) . (46)
This equation of motion should sample an isothermal trajectory at a temperature
determined by the initial internal velocities. Eq. (46) can also be derived directly
from Eq. (18) by imposing ˙T (t) = 0. Using Eqs. (8) and (10) this becomes
˙T (t) =
d
dt
1
kB Nd f
N
i=1
mi ˙r2
i (t) =
2
kB Nd f
N
i=1
mi ˙ri(t) · ¨ri(t) = 0 . (47)
Inserting Eq. (18) and solving for γ (t) leads to Eq. (46).
126 Philippe H. Hünenberger
The Hoover-Evans algorithm should in principle ensure a constant temperature at
all time points. However, this condition is only enforced by zeroing the temperature
derivative. Because the reference temperature To does not appear explicitly in the
scaling factor of Eq. (44), numerical inaccuracies will inevitably prevent temperature
conservation, and cause the temperature to actually drift in simulations.
In the Woodcock algorithm [103], the quantity λ(t; t) in Eq. (20) is found by
imposing temperature conservation in the form T (t + t
2 ) = g
Nd f
To, thereby making
explicit use of the reference temperature. Although g = Nd f seems to be the obvious
choice, it turns out that g = Nd f − 1 is the approptiate choice for the algorithm to
generate a canonical ensemble of conﬁgurations at temperature To (see below). Using
Eqs. (8) and (10), this leads to the condition
λ2(t; t)
kB Nd f
N
i=1
mi ˙ri t −
t
2
+ m−1
i Fi (t) t
2
=
g
Nd f
To . (48)
Solving for λ(t; t) gives12
λ(t; t) =
⎧
⎪⎨
⎪⎩
gkBTo
N
i=1 mi ˙ri t − t
2 + m−1
i Fi(t) t
2
⎫
⎪⎬
⎪⎭
1/2
=
g
Nd f
To
T t + t
2
1/2
. (49)
If T (t − t
2 ) = g
Nd f
To (i.e., if the simulation was started with internal velocities
corresponding to the reference temperature, or otherwise, after a ﬁrst equilibration
timestep), this quantity satisﬁes λ(t; 0) = 1. In this case, applying Eq. (21) gives
γ (t) = (gkBTo)−1
N
i=1
˙ri(t) · Fi (t) . (50)
Inserting into Eq. (18) shows that the equation of motion corresponding to the Woodcock
thermostat is
¨ri(t) = m−1
i Fi(t) − (gkBTo)−1
N
i=1
˙ri(t) · Fi(t) ˙ri(t) . (51)
This equation of motion is rigorously equivalent to the Hoover-Evans equation of
motion (Eq. (46)), provided that the initial internal velocities are appropriate for
12 Note that some simulation programs do not apply the scaling of the velocities by λ(t; t) at
every timestep, but perform the scaling on a periodic basis, or when the difference between
the instantaneous and reference temperatures is larger than a given tolerance [107].
Thermostat Algorithms 127
the temperature To, i.e., that T (0) = g
Nd f
To. However, even in this case, the corresponding
algorithms differ numerically. Because the Woodcock algorithm explicitly
involves the reference temperature To in the calculation of the scaling factor of
Eq. (49), its application removes the risk of a temperature drift.
Because maintaining the temperature constant represents a single constraint
equation involving a total of Nd f velocity variables, it should not be surprising that
numerous other choices of equations of motion lead to an isothermal dynamics (also
satisfying the two constraints that the system linear and angular momenta are constants
of the motion). For example, Haile and Gupta [108] have shown how to construct
two general classes of isothermal equations of motion based on generalized
forces or generalized potentials. An example of the former class is the Hoover-Evans
thermostat. An example of the second class is a thermostat similar (but, contrary
to the claim of the authors [108], not identical) to the Woodcock thermostat. The
equations of motion of this Haile-Gupta thermostat are
˙ri(t) =
g
Nd f
To
T (t)
1/2
˙ri(t) and ¨ri(t) = m−1
i Fi (t) . (52)
In words, the auxiliary velocities ˙ri are propagated independently in time according
to Newton’s second law, and the true velocities obtained by multiplying these by the
appropriate scaling factor at each timestep. In contrast, in the Woodcock thermostat,
the auxiliary velocities ˙ri are obtained at each timestep by increasing the true
velocities ˙ri by m−1
i Fi t.
It can be shown that the ensemble generated by the (identical) Woodcock and
Hoover-Evans equations of motion represents a canonical distribution of conﬁgurations
(though obviously not of momenta) at temperature To, provided that one sets
g = Nd f − 1 ([104, 53]; see Appendix). This may seem surprizing at ﬁrst sight, but
canonical sampling of conﬁgurations is only achieved with T (t) =
Nd f −1
Nd f
To = To,
i.e., when simulating at a slightly lower internal temperature. The reason is that constraining
the temperature effectively removes one degree of freedom from the system.
A more consistent approach would be to alter the deﬁnition of the instantaneous
temperature (Eq. (8)) by changing Nd f to Nd f − 1 in this case. However, since other
thermostats may involve different values for the factor g (see Sect. 3.6), it is more
convenient here to stick to a single deﬁnition of T . On the other hand (and countrary
to the author’s initial claim [108]), the Haile-Gupta thermostat does not generate a
canonical ensemble of conﬁgurations ([53, 109]; see Appendix). The equations of
motion of temperature constraining are smooth, deterministic and time-reversible.
However, the absence of kinetic energy ﬂuctuations may lead to inaccurate dynamics,
especially in the context of the microscopic systems typically considered in sim-
ulations.
128 Philippe H. Hünenberger
3.5
Temperature Relaxation by Weak Coupling
The idea of a thermostat based on a ﬁrst-order relaxation equation is due to Berendsen
[54]. As discussed at the end of Sect. 3, when a system at a given average temperature
T is in contact with a heat bath at a different temperature To, the rate of
temperature change is given by Eq. (26). As discussed in Sect. 3.2, this equation is
only valid when the average temperature is calculated on an intermediate timescale
(short compared to the experimental timescale, but long compared to the time separating
atomic collisions). On this timescale, only the mean effect of the stochastic
forces acting in SD needs to be considered, leading to the ﬁrst-order temperature
relaxation law of Eq. (26).
The idea behind the Berendsen thermostat is to modify the Langevin equation of
motion (Eq. (17)) in the sense of removing the local temperature coupling through
stochastic collisions (random noise), while retaining the global coupling (principle
of least local perturbation). This prescription is equivalent to assuming that Eq. (26)
also applies to the instantaneous temperature T , i.e., that
˙T (t) = τ−1
B [To − T (t)] , (53)
where the appropriate value for τB should be the temperature relaxation time ζT .
In this case, the quantity λ(t; t) in Eq. (20) is found by imposing T (t + t
2 ) =
T (t − t
2 ) + τ−1
B t g
Nd f
[To − T (t − t
2 )], where in principle g = Nd f . Using
Eqs. (8) and (10), this leads to the condition
λ2
(t; t)T t +
t
2
= T t −
t
2
+ τ−1
B t
g
Nd f
To − T t −
t
2
.
(54)
Solving for λ(t; t) gives
λ(t; t) =
T t − t
2
T t + t
2
+ τ−1
B t
g
Nd f
To − T t − t
2
T t + t
2
1/2
≈ 1 + τ−1
B t
g
Nd f
To
T t + t
2
− 1
1/2
. (55)
In general, the algorithm is implemented following the second (approximate) expression.
For either of the two expressions, Eq. (21) gives
γ (t) =
1
2
τ−1
B
g
Nd f
To
T (t)
− 1 . (56)
Inserting into Eq. (18) shows that the equation of motion corresponding to the
Berendsen thermostat is
¨ri(t) = m−1
i Fi(t) −
1
2
τ−1
B
g
Nd f
To
T (t)
− 1 ˙ri(t) . (57)
Thermostat Algorithms 129
In practice, τB is used as an empirical parameter to adjust the strength of the coupling.
Its value should be chosen in a appropriate range. On the one hand, a too large
value (loose coupling) may cause a systematic temperature drift. Indeed, in the limit
τB → ∞, the Berendsen thermostat is inactive leading to the MD equation of motion,
which samples a microcanonical ensemble. Thus, the temperature ﬂuctuations
will increase with τB until they reach the appropriate value for a microcanonical ensemble.
However, they will never reach the appropriate value for a canonical ensemble,
which are larger. For large values of τB, a systematic energy (and thus temperature)
drift due to numerical errors may also occur, just as in MD. On the other hand,
a too small value (tight coupling) will cause unrealistically low temperature ﬂuctuations.
Indeed, the special case of the Berendsen algorithm (Eq. (55)) with τB = t
is the Woodcock thermostat (Eq. (49)), which does not allow for temperature ﬂuctuations.
This shows that the limiting case of the Berendsen equation of motion for
τB → 0 is the Woodcock/Hoover-Evans equation of motion. Values of τB ≈ 0.1 ps
are typically used in MD simulations of condensed-phase systems. Note, however,
that this choice generally leads to ﬂuctuations close to those of the microcanonical
ensemble. With this choice, the Berendsen thermostat merely removes energy drifts
from a MD simulation, without signiﬁcantly altering the ensemble sampled (and thus
rather plays the role of an ergostat algorithm).
The Berendsen equation of motion is smooth and deterministic, but timeirreversible.
The ensemble generated by the Berendsen equations of motion is not
a canonical ensemble ([109]; see Appendix). Only in the limit τB → 0 (or in practice
τB = t), when the Berendsen equation of motion becomes identical to the
Woodcock/Hoover-Evans equation of motion, does it generate a canonical distribution
of conﬁgurations. In the limit τB → ∞, the microcanonical ensemble is recovered.
All intermediate situations correspond to the sampling of an unusual “weakcoupling”
ensemble13, which is neither canonical nor microcanonical [109]. The
reason why the Berendsen thermostat systematically (for all values of τB) underestimates
temperature ﬂuctuations (and thus does not give the correct thermodynamical
ensemble) resides in the transition from Eq. (26) to Eq. (53), corresponding to
the neglect of the stochastic contribution to these ﬂuctuations on the microscopic
timescale.
3.6
Temperature Relaxation by the Extended-System Method
The idea of a thermostat based on an extended-system method is due to Nosé [61]. A
simpler formulation of the equations of motion was later proposed simultaneously14
13 Assuming a relationship of the form of Eq. (115), it is possible to derive the conﬁgurational
partition function of the weak-coupling ensemble as a function of α ([109]; see Appendix)
The limiting cases α = 0 (τB → 0; canonical) and α = 1 (τB → ∞; microcanonical) are
reproduced. Note that the Haile-Gupta thermostat generates conﬁgurations with the same
probability distribution as the Berendsen thermostat with α = 1/2.
14 Eqs. (2.24) and (2.25) in [53] are equivalent to Eq. (6) in [63], provided that one identiﬁes
ζ = ˙s /s . These equations are Eqs. (78) and (79) of the present article.
130 Philippe H. Hünenberger
by Nosé [53] and Hoover [63], so that this algorithm is generally referred to as the
Nosé-Hoover thermostat.
The idea behind the original Nosé [61] algorithm is to extend the real system by
addition of an artiﬁcial (Nd f + 1)th dynamical variable ˜s (associated with a “mass"
Q > 0, with actual units of energy×(time)2, as well as a velocity ˙˜s, and satisfying
˜s > 0) that plays the role of a time-scaling parameter15. More precisely, the timescale
in the extended system is stretched by the factor ˜s, i.e., an inﬁnitesimal time interval
d ˜t at time ˜t in the extended system corresponds to a time interval dt = ˜s−1(˜t) d ˜t
in the real system16. Consequently, although the atomic coordinates are identical in
both systems, the extended-system velocities are ampliﬁed by a factor ˜s−1 compared
to the real-system velocities, i.e.
˜r = r , ˙˜r = ˜s−1
˙r , ˜s = s , and ˙˜s = ˜s−1
˙s . (58)
The Lagrangian for the extended system is chosen to be
Le(˜r, ˙˜r, ˜s, ˙˜s) =
1
2
N
i=1
mi ˜s2˙˜r2
i − U(˜r) +
1
2
Q˙˜s2
− gkBTo ln ˜s , (59)
where g is equal to the number of degrees of freedom Nd f in the real system, possibly
increased by one (see below). The ﬁrst two terms of the Lagrangian represent the
kinetic energy minus the potential energy of the real system (the extended-system
velocities are multiplied by ˜s to recover the real-system ones). The third and fourth
terms represent the kinetic energy minus the potential energy associated with the ˜svariable.
The form of the last term is chosen to ensure that the algorithm produces a
canonical ensemble of microstates (see below). The Lagrangian equations of motion
derived from Eq. (59) read
¨˜ri = m−1
i ˜s−2 ˜Fi − 2˜s−1 ˙˜s ˙˜ri (60)
for the physical variables, and17
15 All extended-system variables will be noted with a tilde overscript, to distinguish them
from the real-system variables (the real-system variable corresponding to ˜s is noted s).
16 To simplify the notation, explicit dependence of the different functions on time is generally
omitted in this section. The time-dependent functions are ˜F(˜t), ˜r(˜t), ˜p(˜t), ˜s(˜t) and ˜ps(˜t)
(together with their ﬁrst and second time derivatives) for the extended system, and F(t),
r(t), p(t), s(t), ps(t), γ (t) and T (t) (together with their ﬁrst and second time derivatives)
for the real system. The dot overscripts indicate differentiation with respect to the extendedsystem
time ˜t for the extended-system variables, and with respect to the real-system time t
for the real-system variables.
17 Because the time-average of the time-derivative of a bounded quantity (for example, ˙˜s)
vanishes, this equation implies
˜s−1
N
i=1
mi ˜s2˙˜r2
i
e,v
= ˜s−1
N
i=1
mi ˙r2
i
e,v
= ˜s−1
e,v
gkBTo , (61)
Thermostat Algorithms 131
¨˜s = Q−1
˜s−1
N
i=1
mi ˜s2˙˜r2
i − gkBTo (63)
for the ˜s-variable. These two second-order differential equations can be discretized
(based on a timestep ˜t in the extended system) and integrated simultaneously during
the simulation. The successive values of r = ˜r and ˙r = ˜s˙˜r describe the evolution
of the atomic coordinates and velocities in the real system at successive time points
separated by t = ˜s−1(˜t) ˜t. This means that the algorithm implemented in this
form (referred to as virtual-time sampling) leads to sampling of the real-system trajectory
at uneven time intervals. A trajectory with real-time sampling can be achieved
either by interpolation at evenly spaced real-time points of the coordinates and velocities
issued from virtual-time sampling, or by rewriting the equations of motion in
terms of the real-system variables (Nosé-Hoover formulation [53, 63]; see below).
As an alternative to Eqs. (60) and (63), the Nosé equations of motion can be
equivalently formulated using a Hamiltonian formalism. In this case, the extendedsystem
conjugate momenta ˜pi and ˜ps associated with the physical degrees of freedom
and with the ˜s-variable are deﬁned as
˜pi =
∂Le(˜r, ˙˜r, ˜s, ˙˜s)
∂˙˜r
= mi ˜s2˙˜ri and ˜ps =
∂Le(˜r, ˙˜r, s, ˙˜s)
∂ ˙˜s
= Q˙˜s . (64)
Comparison with the corresponding real-system momenta18, deﬁned as
pi = mi ˙ri and ps = Qs−2
˙s , (66)
shows that the extended-system momenta are ampliﬁed by a factor ˜s compared to
the real-system momenta. The extended-system Hamiltonian corresponding to the
Lagrangian of Eq. (59) can now be written
where ... e,v denotes ensemble averaging over the extended system (with virtual-time
sampling). Considering Eqs. (8) and (10), this result already suggests that the average
temperature of the real system coincides with To. Using Eq. (101) and pi = mi ˙ri, the
above equation can indeed be rewritten
T e,r =
1
kB Nd f
N
i=1
mi ˙ri
2
e,r
=
g
Nd f
To , (62)
where ... e,r denotes ensemble averaging over the extended system (with real-time sampling).
From Eq. (102), this latter ensemble average can be identiﬁed with a canonical one
when g = Nd f .
18 In the absence of thermostat (i.e., when the variable s is uncoupled from the system), the
real-system Lagrangian may be written
L(r, ˙r, s, ˙s) =
1
2
N
i=1
mi ˙r2
i − U(r) +
1
2
Qs−2 ˙s2 . (65)
The momenta of Eq. (66) are derived from this Lagrangian.
132 Philippe H. Hünenberger
He(˜r, ˜p, ˜s, ˜ps) =
1
2
N
i=1
m−1
i ˜s−2
˜p2
i + U(˜r) +
1
2
Q−1
˜p2
s + gkBTo ln ˜s . (67)
This function is a constant of the motion and evaluates to Ee, the total energy of the
extended system. The corresponding Hamiltonian equations of motion read
˙˜pi = ˜Fi and ˙˜ri = m−1
i ˜s−2
˜pi (68)
for the physical variables, and
˙˜ps = ˜s−1
N
i=1
m−1
i ˜s−2
˜p2
i − gkBTo and ˙˜s = Q−1
˜ps (69)
for the ˜s-variable.
The Nosé equations of motion sample a microcanonical ensemble in the extended
system (˜r, ˜p, ˜t), with a constant total energy Ee. However, the energy of the real system
(r = ˜r, p = ˜s−1 ˜p, t = ˜s−1d ˜t) is not constant. Accompanying the ﬂuctuations
of ˜s, heat transfers occur between the system and a heat bath, which regulate the
system temperature. As will be seen below (Eq. (78) where γ = ˜γ can be identiﬁed
with ˙˜s), the sign of ˙˜s determines the direction of the heat ﬂow. When ˙˜s < 0, heat
ﬂows into the real system. When ˙˜s > 0, heat ﬂows out of the real system. It can be
proved ([61]; see Appendix) that the Nosé equations of motion sample a canonical
ensemble of microstates in the real system, provided that g = Nd f + 1 (virtual-time
sampling) or g = Nd f (real-time sampling), and that Q is ﬁnite, this irrespective of
the actual values of Q and Ee. If the potential energy U(˜r) does not involve terms
giving rise to external forces, the total linear and angular momenta associated with
the physical degrees of freedom in the extended system, namely
N
i=1
˜pi =
N
i=1
mi ˜s2˙˜ri and
N
i=1
˜ri × ˜pi =
N
i=1
mi ˜s2
˜ri × ˙˜ri , (70)
are also conserved. Because ˙˜ri = ˜s−1˙ri (Eq. (58)), this implies that the total linear
and angular momenta of the real system are linearly related to ˜s−1 and thus not
conserved, unless they vanish. This should be the case if the components of the real
velocities ˙ro
i along the external degrees of freedom have been removed by application
of Eq. (11).
The Nosé equations of motion are smooth, deterministic and time-reversible.
However, because the time-evolution of the variable ˜s is described by a second-order
equation (Eq. (63)), heat may ﬂow in and out of the system in an oscillatory fashion
[110], leading to nearly-periodic temperature ﬂuctuations. However, from the
discussion of Sects. 3 and 3.2, the dynamics of the temperature evolution should not
be oscillatory, but rather result from a combination of stochastic ﬂuctuations and exponential
relaxation. At equilibrium, the approximate frequency of these oscillations
can be estimated in the following way [61]. Consider small deviations δ˜s of ˜s away
from the equilibrium value ˜s e,r , where ... e,r denotes ensemble averaging over the
Thermostat Algorithms 133
extended system with real-time sampling. Assuming that the interatomic forces have
a weak effect on the temperature dynamics (as, e.g., in a perfect gas), the quantity
N
i=1 m−1
i ˜p2
i , which is solely altered by the action of the forces (Eq. (68)), can be
assumed nearly constant. In this case, one may write (with g = Nd f )
1
2
N
i=1
mi ˜s2˙˜r2
i =
1
2
N
i=1
m−1
i ˜s−2
˜p2
i ≈
1
2
˜s 2
e,r ˜s−2
Nd f kBTo . (71)
Using this result, Eq. (63) may be written (for small δ˜s)
δ¨˜s = Nd f kBTo Q−1
˜s−1
( ˜s 2
e,r ˜s−2
− 1) ≈ −2Nd f kBTo Q−1
˜s −2
e,r δ˜s . (72)
This corresponds to a harmonic oscillator with frequency
ν = (2π)−1
(2Nd f kBTo)1/2
Q−1/2
˜s −1
e,r , (73)
where19
˜s e,r =
Nd f
Nd f + 1
1/2
exp Nd f kBTo
−1
Ee − H (r, p) , (75)
... denoting a canonical ensemble average. Comparison of the approximate oscillation
frequency ν with the inverse of the temperature relaxation time ζT (Eq. (26))
may guide the choice of parameters Q and Ee leading to a realistic timescale of
temperature ﬂuctuations. If the number of degrees of freedom is large and Ee is
close to H(r, p) , the average canonical energy corresponding to the real system,
Eq. (75) becomes ˜s e,r = 1. The latter condition will be satisﬁed if the algorithm is
initiated using real-system velocities taken from a Maxwell-Boltzmann distribution
(Eq. (25)), together with ˜s(0) = 1 and ˙˜s(0) = 0. In this case, comparing Eq. (73)
with Eq. (26) allows to identify20 1.2(2Nd f kBTo)−1/2Q1/2 with the temperature relaxation
time ζT , i.e., to suggest Q ≈ 1.4Nd f kBToζ2
T as an appropriate value for
simulations. In fact, it may make sense to use an effective relaxation time
τN H = (Nd f kBTo)−1/2
Q1/2
(76)
instead of the (less intuitive) effective mass Q to characterize the strength of the
coupling to the heat bath.
19 Considering Eqs. (98) and (101), one has
˜s e,r =
dp dr d ˜ps d˜s ˜sNd f +1δ[˜s − ˜so]
dp dr d ˜ps d˜s ˜sNd f δ[˜s − ˜so]
. (74)
Inserting Eq. (94), integrating over ˜ps for the numerator and denominator, and setting
g = Nd f leads to Eq. (75).
20 Because cos(1.2) ≈ exp(−1), it is assumed here that the exponential relaxation time is
approximately given by 1.2/(2π)ν−1.
134 Philippe H. Hünenberger
The use of an extended system with a stretched timescale is not very intuitive, and
the sampling of a trajectory at uneven time intervals is rather impractical for the investigation
of the dynamical properties of a system. However, as shown by Nosé [53]
and Hoover [63], the Nosé equations of motion can be reformulated in terms of
real-system variables (together with real-time sampling) so as to avoid these problems.
The transformation from extended-system to real-system variables is achieved
through
s = ˜s , ˙s = ˜s˙˜s , ¨s = ˜s2 ¨˜s + ˜s˙˜s2
,
r = ˜r , ˙r = ˜s˙˜r , ¨r = ˜s2¨˜r + ˜s˙˜s˙˜r ,
ps = ˜s−1
˜ps , ˙ps = ˙˜ps − Q−1
˜s−1
˜p2
s , (77)
p = ˜s−1
˜p , ˙p = ˙˜p − Q−1
˜s−1
˜ps ˜p ,
and F = ˜F .
Because dt = ˜s−1 d ˜t, these equations are derived using d/dt = ˜s d/d ˜t, together
with the deﬁnition of the real-system (Eq. (66)) and extended-system (Eq. (64)) momenta.
Based on these expressions, and deﬁning the quantity γ = s−1 ˙s = Q−1sps,
the Lagrangian equations of motion (Eqs. (60) and (63)) can be rewritten
¨ri = m−1
i Fi − γ ˙ri (78)
and
˙γ = −kB Nd f Q−1
T
g
Nd f
To
T
− 1 . (79)
Note that the variable γ is a real-system variable (i.e., a function of the real-system
time t). The equivalent extended-system variable is ˜γ = γ = ˙˜s. If the effective
coupling time τN H (Eq. (76)) is used instead of Q, Eq. (79) becomes
˙γ = −τ−2
N H
T
To
g
Nd f
To
T
− 1 . (80)
In a similar way, the Hamiltonian equations of motion (Eqs. (68) and (69)) can be
rewritten
˙pi = Fi − γ pi and ˙ri = m−1
i pi (81)
and
˙ps = −kB Nd f T s−1 g
Nd f
To
T
− 1 − γ ps and ˙s = γ s . (82)
Equation (81) is easily identiﬁed with Eq. (78), and Eq. (82) with Eq. (79). The variable
s is absent from the ﬁrst set of equations, i.e., its dynamics has been decoupled.
In the second set of equations, the evolution of the real-system variables is independent
of the actual value of s, i.e., any choice of the initial value s(0) will lead to the
Thermostat Algorithms 135
same dynamics. Finally, it should be stressed that these equations of motion are no
longer Hamiltonian, although the quantity (Eq. (67))
H(r, p, s, ps) =
1
2
N
i=1
m−1
i p2
i + U(r) +
1
2
Q−1
s2
p2
s + gkBTo ln s (83)
is still a constant of the motion (evaluating to Ee). On the other hand, Eqs. (78) and
(79) are still Lagrangian, the corresponding Lagrangian being
L(r, ˙r, s, ˙s) = s [ Le(r, s−1
˙r, s, s−1
˙s) + Ee ] . (84)
To obtain the corresponding Lagrangian equations of motion, Ee is initially treated
as a constant and later expanded using Eq. (83). It appears that Eq. (78) has exactly
the form of Eq. (18), i.e., the Nosé-Hoover thermostat has one equation of motion
in common with both the Woodcock/Hoover-Evans and the Berendsen thermostats.
However, in contrast to these other thermostats where the value of γ was uniquely
determined by the instantaneous microstate of the system (compare Eq. (79) with
Eqs. (45), (50), and (56)), γ is here a dynamical variable which derivative (Eq. (79))
is determined by this instantaneous microstate. Accompanying the ﬂuctuations of
γ , heat transfers occur between the system and a heat bath, which regulate the system
temperature. Because γ = s−1 ˙s = ˜γ = ˙˜s (Eq. (77)), the variable γ in the
Nosé-Hoover formulation plays the same role as ˙˜s in the Nosé formulation. When γ
(or ˙˜s) is negative, heat ﬂows from the heat bath into the system due to Eq. (78) (or
Eq. (60)). When the system temperature increases above To, the time derivative of
γ (or ˙˜s) becomes positive due to Eq. (79) (or Eq. (63)) and the heat ﬂow is progressively
reduced (feedback mechanism). Conversely, when γ (or ˙˜s) is positive, heat is
removed from the system until the system temperature decreases below To and the
heat transfer is slowed down.
The second- and ﬁrst-order Eqs. (78) and (79) can be discretized (based on a
timestep t in the real system) and integrated simultaneously during the simulation.
Note that the Nosé thermostat with g = Nd f + 1 and virtual-time sampling and
the Nosé-Hoover thermostat with g = Nd f formally sample the same trajectory. In
practice, however, the trajectories are sampled at different real-system time points
and will numerically diverge for ﬁnite timestep sizes. It can be proved ([63]; see
Appendix) that the Nosé-Hoover equations of motion sample a canonical ensemble.
in the real system provided that g = Nd f and that Q is ﬁnite, this irrespective of the
actual values of Q and Ee. Though interesting, such a proof is not really necessary
since the Nosé and Nosé-Hoover formalisms are equivalent. As a by-product of this
proof, it is shown that the probability distribution of the γ variable is a Gaussian
of width determined by the parameter Q (Eq. (112)). The Nosé-Hoover equations
of motion are smooth, deterministic and time-reversible. However, just as the Nosé
algorithm, Nosé-Hoover dynamics may lead to temperature oscillations.
In both algorithms, Some care must be taken in the choice of the ﬁctitious mass
Q and extended-system energy Ee. On the one hand, too large values of Q (loose
coupling) may cause a poor temperature control. Indeed, the limiting case of the
136 Philippe H. Hünenberger
Nosé-Hoover thermostat with Q → ∞ and γ (0) = 0 is MD, which generates a microcanonical
ensemble. Although any ﬁnite (positive) mass is sufﬁcient to guarantee
in principle the generation of a canonical ensemble, if Q is too large, the canonical
distribution will only be obtained after very long simulation times. In this case, systematic
energy drifts due to accumulation of numerical errors may interfere with the
thermostatization. On the other hand, too small values (tight coupling) may cause
high-frequency temperature oscillations (Eq. (73)) leading to the same effect. This
is because if the ˜s variable oscillates at a very high frequency, it will tend to be
off-resonance with the characteristic frequencies of the real system, and effectively
decouple from the physical degrees of freedom (slow exchange of kinetic energy).
The choice of the parameters Q and Ee can be guided by comparison of the frequency
ν (Eq. (73)) with the inverse of the temperature relaxation time ζT (Eq. (26)).
Note that if a simulation is initiated with s(0) = 1 and γ (0) = 0, which seems the
most reasonable choice, the value of Ee will match the initial energy of the real system.
The numerical integration of the Nosé and Nosé-Hoover equations will not be
discussed here. A number of alternative schemes have been proposed in the literature
[111, 112, 113, 114, 90].
The constant-temperature Woodcock/Hoover-Evans equation of motion can be
retrieved from the the Nosé-Hoover formalism by a slight modiﬁcation of the
extended-system Hamiltonian. This is done by introducing the constraints
˜s = (gkBTo)−1/2
N
i=1
m−1
i ˜p2
i
1/2
and ˜ps = 0 (85)
into Eq. (67), leading to the modiﬁed Hamiltonian
Hc(˜r, ˜p) =
1
2
gkBTo + U(˜r) +
1
2
gkBTo ln (gkBTo)−1
N
i=1
m−1
i ˜p2
i . (86)
The corresponding Hamiltonian equations of motion for the physical variables in the
extended system are still given by Eq. (68). It can be proved ([61]; see Appendix) that
these equations of motion sample a canonical ensemble of conﬁgurations in the real
system (r = ˜r, ˙r = ˜s−1 ˜p, t = ˜s−1d ˜t) with ˜s given by Eq. (85), provided that g =
Nd f (virtual-time sampling) or g = Nd f −1 (real-time sampling), irrespective of the
constant value Ee of Hc. To show that this situation matches the Woodcock/HooverEvans
equation of motion, the equation of motion must be rewritten in terms of
real-system variables. Evaluating ˙˜s based on Eq. (85) gives
˙˜s = (gkBTo)−1
˜s−1
N
i=1
m−1
i ˜pi · ˙˜pi . (87)
Applying the transformations of Eq. (77), setting γ = s−1 ˙s, and inserting Eq. (68)
for ˙˜p, it is easily seen that the equation of motion in terms of the real-system variables
is Eq. (78) together with
Thermostat Algorithms 137
γ = (gkBTo)−1
N
i=1
˙ri · Fi , (88)
which is identical to the corresponding factor for the Woodcock (Eq. (50)) and
Hoover-Evans (Eq. (45), when T (0) = g
Nd f
To) equations of motion. This also shows
that the appropriate choice for g so as to generate a canonical ensemble of conﬁgurations
using either of these two thermostats is g = Nd f − 1.
If, in addition to incorporating the constraint of Eq. (85), the Hamiltonian is
changed to
Hh(˜r, ˜p) = gkBTo
N
i=1
m−1
i ˜p2
i
1/2
+ U(˜r) , (89)
the corresponding Hamiltonian equations of motion for the physical variables in the
extended system become
˙˜pi = ˜Fi and ˙˜ri = m−1
i ˜s−1
˜pi . (90)
Setting ˙ri = ˙˜ri and pi = mi ˙ri, one recovers (using Eq. (85) and identifying ˙ri =
m−1
i ˜pi) the Haile-Gupta equations of motion (Eq. (52)). Note that ˜s no longer acts
as a scaling parameter and the extended-system Lagragian has been changed, so that
Eqs. (58) and (64) no longer apply. It can be proved ([53]; see Appendix) that this
equation of motion does not sample a canonical ensemble of conﬁgurations in the
real system, irrespective of the choice of g.
3.7
Generalizations of the Previous Methods
An interesting extension of the SD thermostat (Sect. 3.2) is the so-called dissipativeparticle-dynamics
(DPD) thermostat [115, 116, 117, 118]. This scheme retains a key
advantage of the SD thermostat (shared with the Andersen thermostat), namely that
it couples atomic velocities to the heat bath on a local basis (as opposed to the global
coupling applied by all other thermostats discussed in this article). Local coupling
leads to an intrinsically more efﬁcient thermostatization and in turn, permits the use
of longer timesteps to integrate the equations of motion (in applications where this
timestep is not further limited by the curvature of the interaction function, i.e., when
using soft intermolecular potentials). On the other hand, the DPD scheme alleviates
two drawbacks of the SD scheme, namely (i) the non-conservation of the system
linear and angular momentum, and (ii) the loss of local hydronamic correlations between
particles. In practice, this is achieved by replacing the frictional and stochastic
forces acting on individual atoms in SD (Eq. (17)), by corresponding central pairwise
forces acting on atom pairs within a given cutoff distance.
A number of extensions or generalizations of the Nosé or Nosé-Hoover approaches
(Sect. 3.6) have also been reported in the literature, following three main
directions.
138 Philippe H. Hünenberger
First, more general extended Hamiltonians (including the Nosé Hamiltonian as
a particular case) have been shown to also produce canonical phase-space sampling
[119, 120, 121, 122, 123, 124, 125, 126, 127]. The additional ﬂexibility offered
by these new schemes may be used to overcome the non-ergodic behaviour of
the Nosé-Hoover thermostat in the context of small or stiff systems (e.g., single harmonic
oscillator) or systems at low temperatures [63, 128, 129, 121, 130, 131, 132].
A number of these variants can indeed produce the correct canonical distribution for
a single harmonic oscillator [121, 122, 123, 125, 126]. The most popuar of these
schemes is probably the Nosé-Hoover chain thermostat [123], where the single thermostatting
variable γ of the Nosé-Hoover scheme is replaced by a series of variables
thermostatting each other in sequence.
Second, alternatives have been proposed for the Nosé-Hoover scheme, which are
phase-space conserving [133] or even symplectic [127]. One of the latter schemes,
referred to as the Nosé-Poincaré thermostat, leads to the same phase-space trajectory
as the Nosé-Hoover thermostat with real time sampling, but has the advantage of
being Hamiltonian.
Third, generalized equations of motion have been proposed to sample arbitrary
(i.e., not necessarily canonical) probability distributions [134, 135, 136, 137]. Such
methods can be used, e.g., to optimize the efﬁciency of conformationalsearches [134,
135, 137] or for generating Tsallis distributions of microstates [136].
4
Practical Considerations
This section brieﬂy mentions two practical aspects related to the use of thermostats
in MD simulations of (bio-)molecular systems.
The ﬁrst problem is encountered when simulating molecular systems involving
distinct sets of degrees of freedom with either (i) very different characteristic frequencies
or (ii) very different heating rates caused by algorithmic noise. In this case,
the joint coupling of all degrees of freedom to a thermostat may lead to different
effective temperatures for the distinct subsets of degrees of freedom, due to a too
slow exchange of kinetic energy between them. A typical example is the so-called
“hot solvent – cold solute problem” in simulations of macromolecules. Because the
solvent is more signiﬁcantly affected by algorithmic noise (e.g., due to the use of
an electrostatic cutoff), the coupling of the whole system to a single thermostat may
cause the average solute temperature to be signiﬁcantly lower than the average solvent
temperature. A solution to this problem is to couple separately the solute and
solvent degrees of freedom to two different thermostats.
The second problem is encountered when using a simulation program that (incorrectly)
applies the thermostatization directly to the atomic velocities rather than to the
internal (peculiar) velocities (Sect. 3). In this case, the system linear and (under vacuum
boundary conditions) angular momenta are not conserved, unless they exactly
vanish. However, even if these quantities are set to zero at the beginning of a simulation,
numerical errors will unavoidably alter these initial values, permitting a ﬂow
Thermostat Algorithms 139
of kinetic energy between internal and external degrees of freedom. Unfortunately,
it appears that at least some thermostats tend to pump kinetic energy from highfrequency
to low-frequency degrees of freedom (thereby violating equipartition). In
this case, uniform translational and (under vacuum boundary conditions) rotational
motion will tend to build up. For simulations started with vanishing overall momenta,
one generally observes a very slow initial rise (often taking nanoseconds) followed
by a very sudden burst of translational and rotational kinetic energy. The accumulation
of kinetic energy in these degrees of freedom will effectively cool down the
internal ones, giving rise to the so-called “ﬂying ice cube effect” [138, 139]. The
most obvious remedy to this problem is to remove the overall center of mass motion
from the atomic velocities at regular interval during the simulation. However, the
application of thermostatization on the basis of internal velocities (Sect. 3) should
probably be preferred, because it is more consistent and avoids the indrotuction of
discontinuities in the generated velocity trajectory.
5
Appendix: Phase-Space Probability Distributions
Here, the phase-space probability distributions are derived that correspond to the
Woodcock/Hoover-Evans (Sect. 3.4), Haile-Gupta (Sect. 3.4), Berendsen
(Sect. 3.5), Nosé (Sect. 3.6) and Nosé-Hoover (Sect. 3.6) thermostats. The derivations
are given for all but the Berendsen thermostat, for which the result is merely
quoted.
The proof that the Nosé thermostat samples a canonical ensemble of microstates,
provided that g = Nd f +1 (virtual-time sampling) or g = Nd f (real-time sampling),
is as follows [53]. The partition function of the microcanonical ensemble generated
for the extended system using virtual-time sampling (i.e., using the natural time evolution
of the extended system) reads
Ze,v = C d ˜p d˜r d ˜ps d ˜s δ[He(˜r, ˜p, ˜s, ˜ps) − Ee] , (91)
where C is a normalization constant, Ee the (constant) extended-system energy, δ the
Dirac delta function, and He the extended-system Hamiltonian (Eq. (67)). Recasting
this expression in terms of the real-system momenta pi = ˜s−1 ˜pi and substituting
r = ˜r leads to
Ze,v = C dp dr d ˜ps d ˜s ˜sNd f
× δ H(r, p) +
1
2
Q−1
˜p2
s + gβ−1
ln ˜s − Ee , (92)
where H is the real-system Hamiltonian, i.e.
H(r, p) =
1
2
N
i=1
m−1
i p2
i + U(r) . (93)
140 Philippe H. Hünenberger
The argument of the delta function in Eq. (92) has a single zero with respect to the ˜s
variable, namely
˜so = exp −g−1
β H(r, p) +
1
2
Q−1
˜p2
s − Ee . (94)
Using the relationship δ[ f (˜s)] =| f (˜so) |−1 δ(˜s − ˜so), one obtains
Ze,v = Cg−1
β dp dr d ˜ps d ˜s ˜sNd f +1
δ[˜s − ˜so] (95)
= Cg−1
β dp dr d ˜ps
× exp −(Nd f + 1)g−1
β H(r, p) +
1
2
Q−1
˜p2
s − Ee .
Integrating with respect to the variable ˜ps and using the appropriate Gaussian integral
gives
Ze,v = C dp dr exp[−(Nd f + 1)g−1
βH(r, p)] (96)
with
C = C[(Nd f + 1)β]−1/2
(2πgQ)1/2
exp[(Nd f + 1)g−1
βEe] . (97)
Equation (96) shows that the virtual-time extended-system ensemble average of any
quantity A depending on the real-system coordinates r = ˜r and momenta p = ˜s−1 ˜p
(and also possibly on ˜s), deﬁned as
A(r, p) e,v =
dp dr d ˜ps d ˜s ˜sNd f +1 A(r, p)δ[˜s − ˜so]
dp dr d ˜ps d ˜s ˜sNd f +1δ[˜s − ˜so]
, (98)
is equivalent to a canonical ensemble average, i.e.
A(r, p) e,v = A(r, p) when g = Nd f + 1. (99)
If real-time sampling is used instead, the probability of any microstate in the ensemble
is ampliﬁed by a factor ˜s−1 due to the contraction of the timescale. For example,
ten microstates at 1 ps interval in the extended system represent 10 ps of
real-system time if ˜s = 1 but only 5 ps if ˜s = 2. Thus, the larger ˜s, the lower the
real-system weight. Consequently, the real-time ensemble average of the quantity A,
deﬁned from Eqs. (98) and (101) as
A(r, p) e,r =
dp dr d ˜ps d ˜s ˜sNd f A(r, p)δ[˜s − ˜so]
dp dr d ˜ps d ˜s ˜sNd f δ[˜s − ˜so]
, (100)
satisﬁes
A(r, p) e,r = ˜s−1 −1
e,v ˜s−1
A(r, p) e,v , (101)
Thermostat Algorithms 141
irrespective of the value of g. A straightforward consequence of Eqs. (98) and (101)
is that the real-time extended-system ensemble average of any quantity A is equivalent
to a canonical ensemble average, i.e.
A(r, p) e,r = A(r, p) when g = Nd f . (102)
From Eq. (96), the real-system phase-space probability density for the Nosé thermostat
(virtual-time sampling) can be written
ρv(r, p) =
exp[−(Nd f + 1)g−1βH(r, p)]
dp dr exp[−(Nd f + 1)g−1βH(r, p)]
. (103)
The proof that the Woodcock/Hoover-Evans thermostat samples a canonical ensemble
of conﬁgurations provided that g = Nd f −1 follows similar lines [53]. It has
been seen that this thermostat is identical to the Nosé thermostat with the constraints
of Eq. (85) and the Hamiltonian of Eq. (86). In this case, the analog of Eq. (92) reads
Ze,v = C dp dr d ˜s ˜sNd f δ
⎡
⎣g−1/2
β1/2
˜s
N
i=1
m−1
i p2
i
1/2
− ˜s
⎤
⎦
×δ Hc r, ˜sp − Ee ,
= Cgβ−1
dp δ
1
2
N
i=1
m−1
i p2
i −
1
2
gβ−1
× dr d ˜s ˜sNd f −1
δ
1
2
gβ−1
+ U(r) + gβ−1
ln ˜s − Ee
= C dp δ
1
2
N
i=1
m−1
i p2
i −
1
2
gβ−1
× dr d ˜s ˜sNd f δ ˜s − exp −g−1
β
1
2
gβ−1
+ U(r) − Ee
= C dpδ
1
2
N
i=1
m−1
i p2
i −
1
2
gβ−1
dr exp −Nd f g−1
βU(r) .(104)
The second equality follows from the relationship δ[ f (x)] =| f (xo) |−1 δ(x − xo),
with x = 1
2
N
i=1 m−1
i p2
i and xo = 1
2 gβ−1 and inserting Eqs. (85) and (86). The
third equality follows from the relationship δ[ f (˜s)] =| f (˜so) |−1 δ(˜s − ˜so). This
partition function is canonical in the conﬁgurations if g = Nd f (virtual-time sampling).
Considering Eq. (101), Eq. (104) shows that the real-time extended-system
ensemble average of any quantity A depending solely on the coordinates is equivalent
to a canonical ensemble average if g = Nd f −1. From Eq. (104), the real-system
phase-space probability density for the Woodcock/Hoover-Evans thermostat can be
written (real-time sampling)
ρr (r, p) =
δ(1
2
N
i=1 m−1
i p2
i − 1
2 gβ−1) exp[−(Nd f − 1)g−1βU(r)]
dr exp[−(Nd f − 1)g−1βU(r)]
. (105)
142 Philippe H. Hünenberger
The derivation of the phase-space probability distribution for the Haile-Gupta
thermostat follows similar lines [53]. It has been seen that this thermostat is identical
to the Nosé thermostat with the constraints of Eq. (85) and the Hamiltonian of
Eq. (89). In this case, the analog of Eq. (104) reads
Ze,v = C dp dr d ˜s ˜sNd f δ
⎡
⎣g−1/2
β1/2
˜s
N
i=1
m−1
i p2
i
1/2
− ˜s
⎤
⎦
×δ[Hh(r, ˜sp) − Ee]
= Cgβ−1
dp δ
1
2
N
i=1
m−1
i p2
i −
1
2
gβ−1
× dr d ˜s ˜sNd f −1
δ gβ−1
˜s + U(r) − Ee
= C dp δ
1
2
N
i=1
m−1
i p2
i −
1
2
gβ−1
× dr{g−1
β[Ee − U(r)]}Nd f −1
h[Ee − U(r)] , (106)
where h is the Heaviside function. This function arises because ˜s ≥ 0, so that Ee <
U(r) leads to no solution for ˜s. Thus, the real-system phase-space probability density
corresponding to the Haile-Gupta thermostat (virtual-time sampling) is
ρv (r, p) =
δ(1
2
N
i=1 m−1
i p2
i − 1
2 gβ−1){g−1β[Ee − U(r)]}Nd f −1h[Ee − U(r)]
dr {g−1β[Ee − U(r)]}Nd f −1h[Ee − U(r)]
.
(107)
Using Eq. (106), it is easily seen that
˜s e,v = g−1
β[Ee − U(r)] e,v . (108)
Thus, Ee in Eq. (107) can be evaluated as
Ee = gβ−1
s e,v + U(r) e,v . (109)
This distribution function is not canonical, irrespective of the value of g (an alternative
derivation of this result can be found in [109]).
The proof that the Nosé-Hoover thermostat samples a canonical ensemble of
microstates provided that g = Nd f is as follows [63]. Consider the Nosé-Hoover
equations of motion, Eqs. (79) and Eq. (81). Because the variables r, p, and γ are
independent, the ﬂow of the (2Nd f + 1)-dimensional probability density ρ(r, p, γ )
is given by the generalized (non-Hamiltonian) analog of the Liouville equation21
∂ρ
∂t
= −
∂ρ
∂r
· ˙r −
∂ρ
∂p
· ˙p −
∂ρ
∂γ
· ˙γ − ρ
∂
∂r
· ˙r +
∂
∂p
· ˙p +
∂
∂γ
· ˙γ . (111)
21 The generalized Liouville equation states the conservation of the total number of systems
in an ensmble. If ΓΓ = (q, p), this conservation law can be written
Thermostat Algorithms 143
One can postulate the following extended-system phase-space distribution function
ρe,r (r, p, γ ) = C exp −β U(r) +
1
2
N
i=1
m−1
i p2
i +
1
2
Qγ 2
, (112)
where C is a normalization factor. Using Eqs. (58), (79), and (81) the derivatives
involved in Eq. (111) are
−
∂ρ
∂r
· ˙r = βρ
N
i=1
m−1
i Fi · pi
−
∂ρ
∂p
· ˙p = −βρ
N
i=1
m−1
i pi · (Fi − γ pi)
−
∂ρ
∂γ
· γ = βργ kB Nd f T
g
Nd f
To
T
− 1
ρ
∂
∂r
· ˙r = 0
ρ
∂
∂p
· ˙p = −ρNd f γ
ρ
∂
∂γ
· ˙γ = 0 . (113)
Using these results and the deﬁnition of T , it is easily shown that ∂ρ/∂t = 0 in
Eq. (111) provided that g = Nd f . This shows that the extended-system phase-space
density ρe,r (r, p, γ ) is a stationary (equilibrium) solution of Eq. (111) corresponding
to the Nosé-Hoover equations of motion. Integrating out the γ variable leads to the
real-system phase-space probability density for the Nosé-Hoover thermostat (realtime
sampling)
ρr (r, p) =
exp[−Nd f g−1βH(r, p)]
dp dr exp[−Nd f g−1βH(r, p)]
. (114)
which is the canonical probability density if g = Nd f . The Nosé-Hoover equations of
motion are unique in leading to the stationary extended-system phase-space density
of Eq. (112). However, they are not unique in leading to the real-system phase-space
density, because the γ distribution (Gaussian in Eq. (112)) is irrelevant here. Finally,
it should be mentionned that because the Nosé-Hoover thermostat can be derived
from the Nosé thermostat with real-time sampling, the above proof is not really necessary.
It is given anyway as a nice illustration of the use of the generalized Liouville
equation to derive probability distribution functions for thermodynamical ensembles.
∂ρ
∂t
= −
∂
∂ΓΓ
· (ρ ˙ΓΓ ) or
dρ
dt
= −ρ
∂
∂ΓΓ
· ˙ΓΓ . (110)
The two forms can be interconverted by expressing dρ/dt as a total derivative. If the
equations of motion are Hamiltonian, one shows easily that dρ/dt = 0. If ρ(ΓΓ ) is a
stationary (equilibrium) solution, one has ∂ρ/∂t = 0.
144 Philippe H. Hünenberger
The derivation of the phase-space probability distribution for the Berendsen thermostat
follows again similar lines [109]. The ﬁnal expression relies on the assumption
of a relationship
[ K2
− K 2
]1/2
= α(τB)[ U2
− U 2
]1/2
, (115)
between the ﬂuctuations in the kinetic and potential energies in simulations with the
Berendsen thermostat. Clearly, such a relationship exists for any system. However,
it is unclear is whether a common α(τB) applies to all systems, irrespective of their
composition and size. The derivation is then based on ﬁnding a stationary solution
for the generalized Liouville equation (Eq. (111)). The ﬁnal (approximate) result is
(with g = Nd f )
ρb(r, p)=
ρp(p) exp{−β[U(r) − N−1
d f αβ2[U(r)2 − U(r) 2
b]]}
dp ρp(p) dr exp{−β[U(r) − N−1
d f αβ2[U(r)2 − U(r) 2
b]]}
,
(116)
where ρp(p) is the (unknown) momentum probability distribution. Note that the
Haile-Gupta thermostat generates conﬁgurations with the same probability distribution
as the Berendsen thermostat with α = 1/2 [109].
References
1. Baschnagel J, Binder K, Doruker P, Gusev AA, Hahn O, Kremer K, Mattice WL, MüllerPlathe
F, Murat M, Paul W, Santos S, Suter UW, Tries W (2000) Adv Polymer Sci
152:41–156 106
2. van Gunsteren WF, Berendsen HJC (1990) Angew Chem Int Ed Engl 29:992–1023 106,
107, 111
3. Karplus M, McCammon JA (2002) Nature Struct Biol 9:646–652 106
4. Daura X, Glättli A, Gee P, Peter C, van Gunsteren WF (2002) Adv Prot Chem 62:341–
360 106
5. Cheatham III TE, Kollman PA (2000) Annu Rev Phys Chem 51:435–471 106
6. Saiz L, Bandyopadhyay S, Klein ML (2002) Bioscience Reports 22:151–173 106
7. Engelsen SB, Monteiro C, Hervé de Penhoat C, Pérez S (2001) Biophys Chem 93:103–
127 106
8. Warshel A, King G (1985) Chem Phys Lett 121:124–129 107
9. King G, Warshel A (1989) J Chem Phys 91:3647–3661 107
10. Juffer AH, Berendsen HJC (1993) Mol Phys 79:623–644 107
11. Brooks III CL, Karplus M (1983) J Chem Phys 79:6312–6325 107, 120
12. Brünger A, Brooks III CL, Karplus M (1984) Chem Phys Lett 105:495–500 107, 120
13. Beutler TC, Mark AE, van Schaik R, Gerber PR, van Gunsteren WF (1994) Chem Phys
Lett 222:529–539 107
14. Essex JW, Jorgensen WL (1995) J Comput Chem 16:951–972 107
15. Im W, Berbèche S, Roux B (2001) J Chem Phys 117:2924–2937 107
16. Allen MP, Tildesley DJ (1987) Computer simulation of liquids. Oxford University Press,
New York 107, 110, 111
Thermostat Algorithms 145
17. Bekker H (1997) J Comput Chem 18:1930–1942 107
18. Tomasi J, Persico M (1994) Chem Rev 94:2027–2094 107
19. Smith PE, Pettitt BM (1994) J Phys Chem 98:9700–9711 107
20. Roux B, Simonson T (1999) Biophys Chem 78:1–20 107
21. Orozco M, Luque FJ (2000) Chem Rev 100:4187–4225 107
22. Schneider T, Stoll E (1978) Phys Rev B 17:1302–1322 107, 111, 120
23. Schneider T, Stoll E (1978) Phys Rev B 18:6468–6482 107, 111, 120
24. van Gunsteren WF, Berendsen HJC, Rullmann JAC (1981) Mol Phys 44:69–95 107, 111
25. Åkesson T, Jönsson B (1985) Mol Phys 54:369–381 107, 111
26. Börjesson U, Hünenberger PH (2001) J Chem Phys 114:9706–9719
27. Gros P, van Gunsteren WF, Hol WGJ (1990) Science 249:1149–1152 108
28. Gros P, van Gunsteren WF (1993) Mol Simul 10:377–395 108
29. Schiffer CA, Gros P, van Gunsteren WF (1995) Acta Crystallogr D D51:85–92 108
30. van Gunsteren WF, Brunne RM, Gros P, van Schaik RC, Schiffer CA, Torda AE (1994)
Accounting for molecular mobility in structure determination based on nuclear magnetic
resonance spectroscopic and X-ray diffraction data. In: James TL, Oppenheimer NJ (eds)
Methods in enzymology: nuclear magnetic resonance, Vol. 239, Academic Press, New
York, pp 619–654 108
31. Torda AE, Scheek RM, van Gunsteren WF (1989) Chem Phys Lett 157:289–294 108
32. Torda AE, Scheek RM, van Gunsteren WF (1990) J Mol Biol 214:223–235 108
33. Torda AE, Brunne RM, Huber T, Kessler H, van Gunsteren WF (1993) J Biomol NMR
3:55–66 108
34. van Gunsteren WF, Berendsen HJC (1977) Mol Phys 34:1311–1327 108
35. Ryckaert JP, Ciccotti G, Berendsen HJC (1977) J Comput Phys 23:327–341 108
36. Andersen HC (1983) J Comput Phys 52:24–34 108
37. Miyamoto S, Kollman PA (1992) J Comput Chem 13:952–962 108
38. Hess B, Bekker H, Berendsen HJC, Fraaije JGEM (1997) J Comput Chem 18:1463–1472
108
39. Kräutler V, van Gunsteren WF, Hünenberger PH (2001) J Comput Chem 22:501–508
108
40. Tironi IG, Brunne RM, van Gunsteren WF (1996) Chem Phys Lett 250:19–24 108
41. Lebowitz JL, Percus JK, Verlet L (1967) Phys Rev 153:250–254 108
42. Cheung PSY (1977) Mol Phys 33:519–526 108
43. Ray JR, Graben HW (1981) Mol Phys 43:1293–1297 108
44. Ray JR, Graben HW (1991) Phys Rev A 44:6905–6908 108
45. Graben HW, Ray JR (1991) Phys Rev A 43:4100–4103 108, 111, 113
46. Haile JM, Graben HW (1980) J Chem Phys 73:2412–2419 108
47. Haile JM, Graben HW (1980) Mol Phys 40:1433–1439 108
48. Ray JR, Graben HW (1986) Phys Rev A 34:2517–2519 108
49. Ray JR, Graben HW, Haile JM (1981) J Chem Phys 75:4077–4079 108
50. Ray JR, Graben HW (1990) J Chem Phys 93:4296–4298 108
51. Evans DJ, Morriss GP (1983) Chem Phys 77:63–66 108, 124
52. Evans DJ, Morriss GP (1983) Phys Lett 98A 98A:433–436 108, 124
53. Nosé S (1984) J Chem Phys 81:511–519 108, 127, 129, 131, 133, 136, 138, 140, 141
54. Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR (1984) J Chem
Phys 81:3684–3690 108, 121, 127
55. Andersen HC (1980) J Chem Phys 72:2384–2393 108, 117, 122, 123
56. Parrinello M, Rahman A (1980) Phys Rev Lett 45:1196–1199 108
57. Parrinello M, Rahman A (1982) J Phys Chem 76:2662–2666 108
146 Philippe H. Hünenberger
58. Ryckaert JP, Ciccotti G (1983) J Chem Phys 78:7368–7374 108
59. Nosé S, Klein ML (1983) Mol Phys 50:1055–1076 108
60. Heyes DM (1983) Chem Phys 82:285–301 108
61. Nosé S (1984) Mol Phys 52:255–268 108, 129, 132, 136
62. Brown D, Clarke JHR (1984) Mol Phys 51:1243–1252 108
63. Hoover WG (1985) Phys Rev A 31:1695–1697 108, 129, 131, 133, 135, 137, 142
64. Hoover WG (1986) Phys Rev A 34:2499–2500 108
65. Ferrario M, Ryckaert JP (1985) Mol Phys 54:587–603 108
66. Ray JR, Rahman A (1985) J Chem Phys 82:4243–4247 108
67. Melchionna S, Ciccotti G, Holian BL (1993) Mol Phys 78:533–544 108
68. Martyna GJ, Tobias DJ, Klein ML (1994) J Chem Phys 101:4177–4189 108
69. Melchionna S, Ciccotti G (1997) J Chem Phys 106:195–199 108
70. Hünenberger PH (2002) J Chem Phys 116:6880–6897 108
71. Çaˇgin T, Pettitt BM (1991) Mol Phys 72:169–175 108
72. Ji J, Çaˇgin T, Pettitt BM (1992) J Chem Phys 96:1333–1342 108
73. Lo C, Palmer B (1995) J Chem Phys 102:925–931 108
74. Lynch GC, Pettitt BM (1997) J Chem Phys 107:8594–8610 108
75. Fixman M (1974) Proc Natl Acad Sci USA 71:3050–3053 109
76. Lado F (1981) J Chem Phys 75:5461–5463 109
77. Wallace DC, Straub GK (1983) Phys Rev A 27:2201–2205 109
78. Çaˇgin T, Ray JR (1988) Phys Rev A 37:247–251 109
79. Çaˇgin T, Ray JR (1988) Phys Rev A 37:4510–4513 109
80. Lustig R (1994) J Chem Phys 100:3048–3059 109
81. Ray JR, Zhang H (1999) Phys Rev E 59:4781–4785 109
82. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E (1953) J Chem Phys
21:1087–1092 111, 119
83. Valleau JP, Whittington SG (1977) A guide to Monde Carlo for statistical mechanics:
2. Byways. In: Berne BJ (ed) Modern theoretical chemistry, Vol. 5, Plenum Press, New
York, pp 137–168 111, 119
84. Frenkel D (1993) Monte Carlo simulations: A primer. In: van Gunsteren, WF, Weiner PK,
Wilkinson AJ (eds) Computer simulation of biomolecular systems, theoretical and experimental
applications, Vol. 2, ESCOM Science Publishers, B.V., Leiden, The Netherlands,
pp 37–66 111, 119
85. van Gunsteren WF, Nanzer AP, Torda AE (1995) Molecular simulation methods for generating
ensembles or trajectories consistent with experimental data. In: Binder K, Ciccotti
G (eds) Monte Carlo and molecular dynamics of condensed matter systems, Proceedings
of the Euroconference, Vol. 49, SIF, Bologna, Italy, pp 777–788 112
86. Evans DJ, Sarman S (1993) Phys Rev E 48:65–70 112, 117
87. Rugh HH (1997) Phys Rev Lett 78:772–774 113
88. Butler BD, Ayton G, Jepps OG, Evans DJ (1998) J Chem Phys 109:6519–6522 113
89. Hockney RW (1970) Methods Comput Phys 9:136–211 115
90. Morishita T (2003) Mol Simul 29:63–69 116, 135
91. Holian BL, Evans DJ (1983) J Chem Phys 78:5147–5150 117
92. Adams DJ (1975) Mol Phys 29:307–311 120
93. Creutz M (1983) Phys Rev Lett 50:1411–1414 120
94. Ray JR (1991) Phys Rev A 44:4061–4064 120
95. Ciccotti G, Ferrario M, Ryckaert JP (1982) Mol Phys 46:875–889 120
96. Berkowitz M, Morgan JD, McCammon JA (1983) J Chem Phys 78:3256–3261 120
97. Vesely FJ (1984) Mol Phys 53:505–524 120
Thermostat Algorithms 147
98. Guàrdia E, Padró JA (1985) J Chem Phys 83:1917–1920 120
99. van Gunsteren WF, Berendsen HJC (1988) Mol Simul 1:173–185 120
100. Wang W, Skeel RD (2003) Mol Phys 101:2149–2156 120
101. Berkowitz M, McCammon JA (1982) Chem Phys Lett 90:215–217 120
102. Tanaka H, Nakanishi K, Watanabe N (1983) J Chem Phys 78:2626–2634 123
103. Woodcock LV (1971) Chem Phys Lett 10:257–261 124, 125
104. Hoover WG, Ladd AJC, Moran B (1982) Phys Rev Lett 48:1818–1820 124, 127
105. Evans DJ (1983) J Chem Phys 78:3297–3302 124
106. Evans DJ, Hoover WG, Failor BH, Moran B, Ladd AJC (1983) Phys Rev A 28:1016–
1020 124
107. Broughton JQ, Gilmer GH, Weeks JD (1981) J Chem Phys 75:5128–5132 126
108. Haile JM, Gupta S (1983) J Chem Phys 79:3067–3076 126, 127
109. Morishita T (2000) J Chem Phys 113:2976–2982 127, 129, 142, 143, 144
110. Nosé S (1993) Phys Rev E 47:164–177 132
111. Holian BL, De Groot AJ, Hoover WG, Hoover CG (1990) Phys Rev A 41:4552–4553
135
112. Toxvaerd S (1991) Mol Phys 72:159–168 135
113. Martyna GJ, Tuckerman ME, Tobias DJ, Klein ML (1996) Mol Phys 87:1117–1157 135
114. Jang S, Voth GA (1997) J Chem Phys 107:9514–9526 135
115. Hoogerbrugge PJ, Koelman JMVA (1992) Europhys Lett 19:155–160 137
116. Español P, Warren P (1995) Europhys Lett 30:191–196 137
117. Groot RD, Warren PB (1997) J Chem Phys 107:4423–4435 137
118. Soddemann T, Dünweg B, Kremer K (2003) Phys Rev E 68:046702(1)–046702(8) 137
119. Jellinek J, Berry SR (1988) Phys Rev A 38:3069–3072 137
120. Hoover WG (1989) Phys Rev A 40:2814–2815 137
121. Hamilton IP (1990) Phys Rev A 42:7467–7470 137
122. Winkler RG (1992) Phys Rev A 45:2250–2255 137
123. Martyna GJ, Klein ML, Tuckerman M (1992) J Chem Phys 97:2635–2643 137
124. Holian BL, Posch HA, Hoover WG (1993) Phys Rev E 47:3852–3861 137
125. Bra´nka AC, Wojciechowski KW (2000) Phys Rev E 62:3281–3292 137
126. Bra´nka AC, Kowalik M, Wojciechowski KW (2003) J Chem Phys 119:1929–1936 137
127. Laird BB, Leimkuhler BJ (2003) Phys Rev E 68:016704(1)–016704(6) 137
128. Posch HA, Hoover WG, Vesely FJ (1986) Phys Rev A 33:4253–4265 137
129. Jellinek J, Berry SR (1989) Phys Rev A 40:2816–2818 137
130. Toxvaerd S (1990) Ber Bunsenges Phys Chem 94:274–278 137
131. Calvo F, Galindez J, Gadéa FX (2002) J Phys Chem A 106:4145–4152 137
132. D’Alessandro M, Tenenbaum A, Amadei A (2002) J Phys Chem B 106:5050–5057 137
133. Winkler RG, Kraus V, Reineker P (1995) J Chem Phys 102:9018–9025 137
134. Pak Y, Wang S (1999) J Chem Phys 111:4359–4361 137
135. Fukuda I (2001) Phys Rev E 64:016203(1)–016203(8) 137
136. Fukuda I, Nakamura H (2002) Phys Rev E 65:026105(1)–026105(5) 137
137. Barth EJ, Laird BB, Leimkuhler BJ (2003) J Chem Phys 118:5759–5768 137
138. Harvey SC, Tan RKZ, Cheatham III TE (1998) J Comput Chem 19:726–740 138
139. Chiu SW, Clark M, Subramaniam S, Jakobsson E (2000) J Comput Chem 21:121–131
138
Index
boundary condition, 107
hard, 107
soft, 107
spatial, 107
thermodynamical, 108
Ensemble, 108
ensemble
canonical, 108
generalized, 108
grand-canonical, 108
grand-isothermal-isobaric, 108
grand-microcanonical, 108
isoenthalpic-isobaric, 108
isothermal-isobaric, 108
microcanonical, 108
thermodynamical, 108
force
frictional, 107, 108
stochastic, 107, 108
Hamiltonian
rotational-invariance, 108
time-independence, 108
translational-invariance, 108
molecular dynamics, 107