Nose-Hoover chains: The canonical ensemble via continuous dynamics
Glenn J. Martyna and Michael L. Klein
Department o/Chemistry, University 0/Pennsylvania, Philadelphia, Pennsylvania 19104-6323
Mark Tuckermana)
Department 0/Chemistry, Columbia University, New York, New York J0027
(Received 8 April 1992; accepted 4 May 1992)
Nose has derived a set of dynamical equations that can be shown to give canonically
distributed positions and momenta provided the phase space average can be taken into the
trajectory average, i.e., the system is ergodic [So Nose, J. Chem. Phys. 81, 511 (1984),
W. G. Hoover, Phys. Rev. A 31, 1695 (1985)]. Unfortunately, the Nose-Hoover dynamics is
not ergodic for small or stiff systems. Here a modification of the dynamics is proposed
which includes not a single thermostat variable but a chain of variables, Nose-Hoover chains.
The "new" dynamics gives the canonical distribution where the simple formalism fails.
In addition, the new method is easier to use than an extension [D. Kusnezov, A. Bulgac, and
W. Bauer, Ann. Phys. 204, 155 (1990)] which also gives the canonical distribution for
stiff cases.
I. INTRODUCTION same initial conditions. The present scheme preserves the
advantages of the original Nose-Hoover approach.
Temperature control in molecular dynamics simulations
(Newtonian dynamics) is essential as most experimental
measurements are performed at constant temperature
rather than constant energy. Nose resolved this
difficulty by deriving a dynamics from an extended Hamiltonian
that can be shown to give canonically distributed
positions and momenta. l
-4 The proof, however, involves
the assumption that the dynamics itself is ergodic or, more
succinctly, that the trajectory average can be taken into the
phase space average. It has since been shown that for small
or stiff systems the dynamics is not ergodic and the correct
distributions are not generated.2
,5 The method does, however,
work extremely well in large (ergodic) systems.
In this article, a possible reason for the nonergodicity
of Nose's original formulation is discussed. A simple modification
based on this analysis is presented and tested on
model problems. The method proposed herein succeeds
precisely where the original method fails. The present
modification has the advantage that it preserves the simplicity
of the original approach.
Other methods can also be used to obtain canonical
distributions on stiff systems,6--9 some of which are
tions and generalizations of the original Nose-Hoover ap-
proach.7
-
9 In the general method of Kusnezov et ai., the
choice of functions and parameters is somewhat arbitrary
and the algorithm difficult to apply in complex situations
such as constant pressure simulations. The method of
Andersen which involves stochastic collisions (at intervals
some or all of the velocities are resampled according to the
Boltzman distribution) will also give the canonical ensemble
even in the most trivial systems.6
This method has the
disadvantage that it is not a continuous dynamics with well
defined conserved quantities. Hence the same trajectory
cannot easily be reproduced by another worker from the
"In partial fulfillment of the Ph.D. in the Department of Physics, Columbia
University.
II. METHODS
The set of dynamical equations,1-4
(1)
defines Nose-Hoover dynamics. Here Pi and qi are one
dimensional variables and the variable '1], which is decoupled
from the dynamics, is included for completeness. This
set of equations, can be shown to give the canonical distribution
in an ergodic system.1
,2 A proof due to Hoover uses
conservation of probability2 to show that the distribution
[
1 ( N P7l(p,q,P71,'1]) o::exp -kT V(q)+ 2m/2Q
is stationary
al N al. al. al. al.
at +.2: aq. qi+an.Pi+an P71+a'YI '1]1=1 1 1:'1 1:'71 .,
[
aqi api aT] ap71]
+1 -+-+-+- =0.
aqi api a'1] ap71
(2)
(3)
Therefore (if the system is ergodic), l(p,q,P71''1]) is the
static probability distribution generated by the dynamics.
The proof is necessary but not sufficient for a general system.2-4,lO
In addition, it can be shown that the quantity
J. Chern. Phys. 97 (4), 15 August 1992 0021-9606/92/162635-09$006.00 @ 1992 American Institute of Physics 2635
Downloaded 28 Dec 2011 to 128.122.251.131. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
2636 Martyna, Klein, and Tuckerman: Nose-Hoover chains
N 2 2
P; PTJ
H'(p,q,7]'PTJ) = V(q) + 2m;+2Q+NkT7] (4)
is conserved, namely,
(5)
This is why it is useful to include the variable 7] in Eq. (1).
In Appendix A, the extension due to Kusnezov et al. and
the method of Winkler are discussed.3
,7,9
The proof presented above is valid only for ergodic
systems. It does not guarantee that the system is ergodic
and that the correct limiting distribution will be generated.
Indeed, in some cases, this distribution is not produced.2-4
Let us consider how this situation can be improved. The
distribution itself has a Gaussian dependence on the particle
p, as well as thermostat momenta, PTJ' While
the GauSSian fluctuations ofP are driven with a thermostat,
there is .nothing to drive the fluctuations ofPTJ' Although
the postlons and momenta (the q's and p's) are of primary
interest, ifthe Nose-Hoover equations are ergodic, the system
covers the whole phase space, which includes the space
of the thermostat velocities. The fluctuations of the thermostat
variables which clearly occur in ergodic systems
may even be important in driving the system to fill phase
space (the dynamical equations are, of course, coupled).
This suggests thermostating PTJ and, by analogy, the thermostat
of PTJ plus its thermostat, etc., to form a chain.
Other methods of "shaking up" the.system have also been
suggested.7-9,11
The "new" dynamics, the Nose-Hoover chain method,
can be expressed as
. P;
q;=-,
m;
(6)
where M thermostats have been included. The 7]'S are
again presented for completeness. These equations can be
shown to have the phase space distribution
(7)
and the conserved quantity
N M p2
I TJi
H'(p,q,7],PTJ
) =:= V(q) + "c., -+ "c., -+NkT7]1
;=1 2m; i=1 2Q;
M
+ L kT7];. (8)
;=2
Note that even in large systems, the addition of the extra
thermostats is relatively inexpensive as they form a simple
one dimensional chain. Only the first thermostat interacts
with N particles.
A good choice for the thermostat masses will help the
dynamics achieve the canonical distribution.1
-4 If very
large masses are chosen, a distribution consistent with the
microcanonical ensemble may result. If very small masses
are chosen, the fluctuations of the momenta may be greatly
inhibited. In Appendix B, a short argument is presented
which suggest that the masses be taken to be Q1
= NkT/oi and Qj = kT/oi. This choice allows the thermostats
to be in approximate resonance with both the system
variables (thep's and q's), which are assumed to have
fundamental frequency lU, and each other. As will be
shown in the results section, the mass choice in the chain
method is much less critical than in both the simple
metllod2 and some of the newer methods.9
The method presented above increases the size of the
phase space and thus helps make the system ergodic. There
are, however, several issues that must be examined. Consider
a point in phase space (q,p) with
aV(q)
--=0
aq; ,
.P;=o. (9)
At such a point, the equations of motion for the NoseHoover
chain dynamics give dqn/ dt' = 0 for all n, which
effectively stops the dynamics of the positions and momenta
(Note this is equally true in Newtonian dynamics).
We call such points Hoover holes, after Hoover who first
began to look at the ergodicity of small systems in the
modified dynamics. For the dynamics to be ergodic, the
system must come infinitesimally close to the Hoover holes
without actually visiting them. Therefore, the dynamics
near the holes must be examined to be sure that the holes
can reasonably be assumed to be of measure zero.
Hoover holes associated with specific values of the
thermostat coordinates are fixed points of the entire NoseHoover
chain dynamics. Ifa fixed point is stable, then in its
vicinity, the equations will have exponentially damped solutions
which will drive the system into the point. This
scenario is disastrous if an ergodic system is desired. The
stabi?ty o.f a fixed p.oint can be determined by examining
the lmeanzed equatIons motion about the point.12,13 The
J. Chern. Phys., Vol. 97, No.4, 15 August 1992
Downloaded 28 Dec 2011 to 128.122.251.131. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
Martyna, Klein, and Tuckerman: Nose-Hoover chains 2637
analysis proceeds by rewriting the Nose-Hoover equations
as the convenient first order system
. Pi
qi=-'
mi
(10)
where the variables (;j = TJj have been introduced. The fixed
point condition for the system is then given by the set of
algebraic equations obtained from setting each of the time
derivatives on the left side of Eqs. (11) equal to zero. This
condition can be expressed as
aV(q)
---0
aqi - ,
Pi=O,
{;1{;2+a =0,
N{;i-a-{;2{;3=0,
{;;_l-a-{;!j+l =0,
(11 )
where the mass choice, QI = NQ and Qj = Q, has been introduced
and a=kTIQ. The linearized equations themselves
are obtained by differentiating each quantity on the
right side of Eq. (10) with respect to each variable in the
system. The stability matrix, A, is defined by the quantities
resulting from this differentiation evaluated at a fixed point
[Le., a solution of Eqs. (11)]. The eigenvalues of this matrix
give the stability condition of the point. If the eigenvalues
all have negative real parts, then the fixed point is
purely stable. Now if the trace of a matrix is zero, the sum
of the eigenvalues is zero. If the sum of the eigenvalues is
zero, it can be concluded that either the eigenvalues are all
zero or have both positive and negative real parts. Examination
of Eqs. (10) and (11) and the linearized motion
shows that the there are no solutions with all eigenvalues
equal to zero, so we conclude that the latter must be true.
A proof that the trace ofA is zero is presented in Appendix
C. Therefore, the fixed points of chain dynamics are not
purely stable. It is also significant to note that the rate of
change of the total phase space volume in the vicinity of
the fixed point is related to the trace of the stability matrix
by
(api aXl) aSj
-+- +i=1 api ax, j=1 a{;j
N M
L Aii+ L Aj+N,j+N=Tr(A). (12)
i=1 j=1
Thus if the trace of A is zero, then Liouville's theorem
(conservation of phase space volume) holds in the vicinity
of the fixed points.14 This implies that the fixed points are
neither attractors or repellorsy,15 These arguments only
apply for M> 1, for M = 1 or the usual Nose-Hoover dynamics
there are no fixed points of the total dynamics.
However, the rate of change of the total phase space volume
is not everywhere zero, which suggests that the NoseHoover
chain dynamics does not have an underlying
Hamiltonian structure. This fact could also be deduced
from the stability conditions on the fixed points.13 That is,
the Nose-Hoover dynamics admit spiral fixed points which
are forbidden in Hamiltonian systems.
The preceding analysis cannot determine whether the
system is ergodic. It does, however, give indications of
when a system may not behave well in this respect. It is
included for completeness and to indicate that NoseHoover
chain dynamics has some reasonable and desirable
properties. A more useful analysis would involve showing
that no periodic orbits exist in the dynamics, a periodic
orbit stability analysis. This at present can only be done
numerically. In the results section it will be shown that
Nose-Hoover chain dynamics fills phase space for some
reasonable values of M.
The Lyapunov exponent gives a measure of the degree
ofchaos present in a dynamical system.12,13,15,16 In general,
the more chaotic the dynamics of a system, the more
quickly it fills phase space. It is therefore important to
study this quantity in chain dynamics. The calculation of
Lyapunov exponents is based on dynamics cast in the generic
form
r(t)=F(r), (13)
where ret) refers to a point in phase space [which for
chain dynamics is (q,p'P7])]. In the first method used to
calculated the exponents (Method I), two nearby trajectories
are integrated for a small time interval T, and the
distance between them monitored. The initial separation is
determined by
(14)
where or(O) is a vector of norm E. After an interval T, the
norm of IOr(T) I is computed and saved. The vector
or(T) is then renormalized to € and the process repeated.
The Lyapunov exponent is calculated from
i 10glor/T) I.N,. j=1 E
(15)
In the second method (Method II), the linearized equations
about the trajectory r(t)
d
dt orct)=M(t)or(t), (16)
J. Chern. Phys., Vol. 97, No.4, 15 August 1992
Downloaded 28 Dec 2011 to 128.122.251.131. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
2638 Martyna, Klein, and Tuckerman: Nose-Hoover chains
Br(t) =Texp [ J:M(t')dt' ]Br(O) (17)
are solved as the trajectory ret) evolves. Here the matrix
M is defined to be Mij = (aF/ar) Ir=r(t), and T is the
time ordering operator. At time 'T, Br is renormalized and
the process repeated. The exponent is calculated from the
sum of the log of norms presented above in Eq. (15). In
general, computing the exponent by these two methods
serves as an excellent check on the results, however, it has
been shown that Method I is consistently the more reli-
able.16
Assuming the dynamics produces the correct distribution,
it is desirable to know how rapidly it samples the
distribution. A measure of this is the rate ofconvergence of
the time average of the potential and kinetic energy to their
ensemble averages (the p's and q's are the variables of
primary interest). Such a rate can be quantified by calculating
the number of time steps necessary to obtain an error
bar or error in the mean of desired tolerance. This number
of time steps can then be compared to the number of
Monte Carlo steps needed to give the same error bar where
the Monte Carlo calculation is performed by directly sampling
the distribution of interest (i.e., the perfect stochastic
calculation). An efficiency can be calculated which is the
ratio of the number of time steps to the number of Monte
Carlo steps R = SMDISMC' Thus the larger is R, the less
efficient the method where R has a lower limit of 1. The
error bars for the time averages are calculated using the
block averaging technique17
which takes into account correlations
which may exist in the data. The Monte Carlo
data is uncorrelated by definition. The error bar of the
average or error in the mean is then given analytically as
the standard deviation of the quantity (averaged over the
distribution) divided by the square root of the number of
steps. The efficiency obviously depends on the method of
integration and the time step chosen to perform the dynamics
but does give a useful measure convergence. Two
values of the efficiency are calculated, R'l' which measures
the convergence of the average potential energy and Rp,
which measures the convergence of the average kinetic en-
ergy.
The velocity Verlet integrator18
is used to integrate the
Nose-Hoover chain equations. The position equations are
deterministic and the velocity equations are solved iteratively
to a convergence level of 10-14• The model problems
studied in the next section were integrated with time steps
of at=0.OI6, at=0.OO5, at=0.OOI6, and at=0.OOO5 for
systems using Qi = 10, Qi = 1, Qi = 0.1, and Qi = 0.01, respectively.
These choices give energy (H') conservation to
a few parts in 105
• {For reference, integrating the harmonic
oscillator [m= I,w= I,q(O) = I,v(O) = 1], to the same tolerance
requires a time step ofO.OI6.} A great advantage of
the chain method is that there is little or no degradation in
energy conservation upon the addition of thermstats beyond
M= 1. All runs were 1.5X 106
time steps. In the calculation
of the Lyapunov exponents, we compared the exponent
for choices of'T ranging from lOat to l00at, and
2
> 0
-1
-2
.4
..-..
><........
Il..
.2
0
.4
..-..
>........
Il..
.2
o
-1
(b)
-4 -2
(e)
-4 -2
o
X
0
X
o
V
2 4
2 4
FIG. 1. (a) Density plot of the Nose-Hoover dynamics of a harmonic
oscillator [q(O) = I,p(O) = I,P'1(O) = I,Q = 1]. (b) Position distribution
function obtained from Nose-Hoover dynamics of a harmonic oscillator
(dotted line). The solid line is the exact result. (c) Velocity distribution
function obtained from Nose-Hoover dynamics of a harmonic
oscillator (dotted line). The solid line is the exact result.
found good agreement between Methods I and II in this
range.
III. RESULTS
Two systems were chosen to illustrate the new method.
The first, is a one dimensional harmonic oscillator (m
= l,w = I,Qi = 1) with initial condition [q(O)
= O,p(O) = I,Pl1/0) = 1] and kT= 1. In Fig. 1, a density
map and the projected distribution functions are presented
for the usual Nose-Hoover dynamics (M= 1). The dynamics
does not fill space. Also, if the initial conditions are
changed, (q = O,p = I,Pl1 = 10), the results are changed
(see Fig. 2). This is unacceptable as an invariant probability
distribution is desired. Similar results have been found
by others.2
,7 The Nose-Hoover chain dynamics, (M=2),
gives rather different results (see Fig. 3). The distribution
functions seem to be good approximations to the canonical
results and the dynamics fills space. Changes in the initial
conditions did not have an appreciable effect on the results.
J. Chern. Phys., Vol. 97, No.4, 15 August 1992
Downloaded 28 Dec 2011 to 128.122.251.131. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
Martyna, Klein, and Tuckerman: Nose-Hoover chains 2639
4
2
:> 0
-2
-4
-4 -2 0 2 4
X
,
(b)
,
.'
.6 .'.' ,.','
,", .,.......
>< .4
---Il..
.2
0
-4 -2 0 2 4
X
.6
(e)
r'
--- ,
:>
..4
---p..,
.2
0
-4 -2 0 2 4
V
FIG. 2.(a) Density plot of the Nose-Hoover dynamics of a harmonic
oscillator [q(O) = l,p(O) = I,P'1(O) = IO,Q = 1]. (b) Position distribution
function obtained from Nose-Hoover dynamics of a harmonic
oscillator (dotted line). The solid line is the exact result. (c) Velocity
distribution function obtained from Nose-Hoover dynamics of a harmonic
oscillator (dotted line). The solid line is the exact result.
In addition, the choice of thermostat mass is not critical in
this method unlike both the original2
and some of the
newer methods.9
Rather, for all the values attempted, Q
=100, M=5, Q=lO, M=2, Q=O.I, M=2, and Q=O.OI,
M =2, the canonical distribution was generated.
The Lyapunov exponents for systems containing M
=1-15 thermostats were calculated for wide variety ofinitial
conditions (Q= 1). The two methods used to obtain
the exponents (see Methods) were found to be in good
agreement as is shown in Tables I and II. In Fig. 4, the
exponents are plotted as a function of M. The exponents
become increasing large with M and around M =4,5 become
competitive with those determined by Kusnezov et
al. for their method.7
However, the Lyapunov exponents
include information about the dynamics of thermostat
variables which are not of primary interest. In Table III,
the efficiency of chain dynamics in the convergence of the
average potential and kinetic energy is calculated for a
variety of parameters (see Sec. II). The results indicate
that the parameter set Q= I,M= 5 converges most rapidly.
4
2
:> 0
'.
-2
-4
-2 0 2
X
.4
.3
........
>< .2
---Il..
.1
0
-4 -2 0 2 4
X
.4
.3
---:>
--- .2Il..
.1
0
-4 -2 0 2 4
V
FIG. 3.(a) Density plot of the Nose-Hoover chain dynamics (M=2,Q
= 1) of a harmonic oscillator [q(O) = I,p(O) = I,P'1(O) = 1]. (b)
Position distribution function obtained from Nose-Hoover dynamics of a
harmonic oscillator (dotted line). The solid line is the exact result. (c)
Velocity distribution function obtained from Nose-Hoover dynamics of a
harmonic oscillator (dotted line). The solid line is the exact result.
For this set, the chain method is found to be competitive
with the method ofKusnezov eta!. (Qp = I,Qq = I) for the
convergence of the potential energy (see Table III). The
average kinetic energy converges very rapidly in the latter
TABLE I. Comparison between Methods I and II of the convergence of
the Lyapunov exponent for M =4, and T= lO8.t.
Steps A,(Method I) A,(Method II)
100 0.2983 0.2954
1000 0.2517 0.2516
2000 0.2659 0.2660
3000 0.2290 0.2293
4000 0.2196 0.2200
5000 0.2371 0.2370
6000 0.2465 0.2460
7000 0.2391 0.2392
8000 0.2381 0.2382
9000 0.2402 0.2402
10000 0.2385 0.2384
J. Chern. Phys., Vol. 97, No.4, 15 August 1992
Downloaded 28 Dec 2011 to 128.122.251.131. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
2640 Martyna, Klein, and Tuckerman: Nose-Hoover chains
TABLE II. Comparison between Methods I and II of the convergence of
the Lyapunov exponent for M =6, and T= lOat.
Steps A(Method I) A(Method II)
100 0.4608 0.4571
1000 0.4366 0.4368
2000 0.3829 0.3832
3000 0.3359 0.3356
4000 0.3273 0.3268
5000 0.3139 0.3129
6000 0.2939 -0.2936
7000 0.3025 0.3028
8000 0.2966 0.2968
9000 0.2997 0.2998
10000 0.2890 0.2897
method which may be due to the use of the third power of
the momentum thermostat in the coupling term (see Appendix
A). The calculated efficiency is a subjective measure
that depends on the method of integration and the
time step used.
For a harmonic oscillator, the Hoover holes do not
have a pathological effect. The phase space in a harmonic
oscillator is such that p =0 can be visited when and
vice versa, which reduces the effect of the hole particularly
on the integrated distributions [P(v) and P(q)]. To complete
the analysis, pairs of trajectories within a radius of
10-8
about the hole were examined. Lypaonovexponents
similar to those reported in Fig. 4 for a given value of M
were obtained (Q= 1). The fact that a positive exponent is
obtained implies that the holes are not attractive and seem
to be of zero measure, consistent with the analysis in Sec.
II.15
The second system studied is a one dimensional free
particle (Q= 1). The Nose-Hoover chain dynamics can
only give canonically distributed velocity on the half space
.5
.4
.3
..<
.2
.1
0
0 5 10 15
M
FIG. 4. Lyapunov exponent for the chain dynamics of a harmonic oscillator
as a function of the number of thermostats, M (Q= 1).
TABLE III. Efficiency of chain dynamics for the harmonic oscillator."
Q M Rq Rp
1.0 2 230 230
1.0 5 230 110
1.0 15 230 110
10.0 5 370 370
0.1 5 730 70
0.01 5 2300 150
"For the method of Kusnezov et af. (Qp = l,Qq = 0, Rq = 230 and Rp
=3.
v>O. When v=O, the dynamics stops for all times
[dnq(t)/dt" = ofor all n]. The dynamics was never found to
stop, but neither could it pass through the Hoover hole.
(The holes seem to be of measure zero and are not attractive.)
The dynamics is symmetric about the hole since the
equations of motion do not depend on the sign of v. Thus
an identical trajectory can be generated on the other half of
phase space with no cost. In Fig. 5, the velocity distribution
of the free particle is presented for M = 1,2,3. The
canonical distribution is recovered when M =3. As discussed
in Appendix A, the functions chosen by Kusnezov
et al. for use in their method are not appropriate for this
system.
IV. DISCUSSION
The idea of thermostating the extended variable is potentially
quite powerful. In stiff complex systems such as
proteins, it is difficult to start near equilibrium. In such
cases, large unphysical oscillations in the temperature may
develop. It is expected that additional thermostats will effectively
damp such oscillations. Similarly, oscillations can
develop in the volume in constant pressure-constant tem-
.8
.6
.4
.2
o r----?,
-2 o
V
----M=3
- - - - --- M=2
---·-------M=l
---Exact
2 4
FIG. S. Velocity distribution function obtained from Nose-Hoover chain
dynamics (M= 1,2,3,Q= I) of a free particle. The solid line is the exact
result.
J. Chern. Phys_, Vol. 97, No.4, 15 August 1992
Downloaded 28 Dec 2011 to 128.122.251.131. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
Martyna, Klein, and Tuckerman: Nose-Hoover chains 2641
perature simulations. As the momentum of the extended
variable that drives volume is distributed canonically, it
can also be thermostated (the proof goes through straightforwardly).
Again, this should help damp the unphysical
oscillations resUlting in more stable simulations.
In summary, a modification of Nose-Hoover dynamics
which we call Nose-Hoover chain dynamics has been
shown to give a very good approximation to the canonical
ensemble even in pathological cases. The idea of thermostating
extended variables will likely find wide application.
ACKNOWLEDGMENTS
The research described herein was supported by the
National Science Foundation under Grant No. CHE-
8722481. One of us (G.M.) would like to acknowledge an
NSF Posdoctoral Research Associateship in Computational
Science and Engineering (ASC-91-08812). In addition,
we would like to thank Bruce Berne for many insights
incorporated into this paper. We would also like to thank
Professor P. Pechukas and Professor L. Schulman for useful
discussions on chaos theory, and Professor W. Hoover,
Professor S. Nose, and Dr. J. Jellinek for commenting on
the paper before publication.
APPENDIX A
The Nose-Hoover equations are in fact a subset of the
more general set of equations,1,7
(Al)
where Frp Grp hrp hp are arbitrary functions. These equations
have the conserved quantity
(A2)
where h/ = dg/ld'l}/. This is a powerful generalization that
can be used to generated canonical dynamics for variety
systems (i.e., Lie algebras3
,7). The Nose-Hoover equations
are generated from the specific choice Fq = hq = 0 and gp
='I};;2, hp = 'I}, Fp = p. Such a choice is more general than
it may first appear as can seen by examining the position
dependence of the functions Fq and Fp' This dependence
must be consistent with the boundary condition of a given
problem. Therefore, a general purpose method must take
these functions to be independent of position. Furthermore,
the proof that the dynamical system, Eq. (1), gives
rise to the canonical distribution relies on the identities,3,7
(A3)
kT / 8Fp (hQi) ) = / F ( . .) 8H(p,q) )
\ 8Pi \ P PI,ql 8Pi '
where the average is, itself, over the canonical distribution.
If the Fs are independent of position, Fq must be taken to
be zero. This argument leads naturally to a Nose-Hoover
form if the dependence of Gip) and hi'l}p) is chosen to be
linear, i.e., the simplest nontrivial choice. Chain dynamics
maintains this form. It is, however, possible to derive a
general set of equations from which chain dynamics
emerges as a specific choice of functions.
In studying simple bound problems, Kusnezov et al.
have advocated the choice, Fq = q3, Fp = p, hq = 71q' hp
= 71;.7
With these choices the trace ofthe stability matrix is
zero which guarantees certain properties of the dynamics
in the vicinity of the fixed points (see Sec. II). The stability
matrix of chain dynamics also has a zero trace (see Appendix
C). In addition, the fixed points of chain dynamics
occur at their natural places [Pi = 0,8V(q)/8qi = 0] while
fixed points of the more complex method do not.
The equations of motion [Eq. (1)] with the choice of
functions of Kusnezov et al. are rather stiff. For a harmonic
oscillator (m= 1,cu= 1,kT= 1) integration by velocity
VerIet did not conserve H' well. In the spirit of the
iterative scheme used to integrate the Nose-Hoover chains,
the set offirst order differentials in Eq. (1) were integrated
by iterating the general form
. . dt
b(t) =b(O) + [b(O) +b(t)] 2" (A4)
to convergence. A time step of 0.001 was needed to obtain
reasonable H' conservation for the oscillator. About five
iterations were needed for convergence.
The dynamics proposed by Winkler9
can be expressed
as
(A5)
PTJ= [ - (N-!)kT].
71 i=1 mi
This can be shown to have the distribution function
1 { 1 [ N ] }
f(p,q,PTJ,71) cc7jexp -kT V(q)+ 2m/2Q (A6)
with conserved quantity
J. Chern. Phys., Vol. 97, No.4, 15 August 1992
Downloaded 28 Dec 2011 to 128.122.251.131. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
2642 Martyna, Klein, and Tuckerman: Nose-Hoover chains
p2
H'=H(P,q)+2Q-(2N-1)kTlog Tf. (A7)
The variable Tf appears explicitly in the dynamical equations
but has an "unbounded" probability distribution.
Therefore, the dynamical equations appear to be rather
poorly behaved in some regions ofphase space, particularly
when Tf is small. These regions are not sampled by the
dynamics as the conserved quantity restricts phase space
such that Tf==exp{[H(p,q) + p;I2Q - E]/(2N
- 1)kT}. As E is a constant and H(p,q) and PT/ have
bounds (their distribution functions decay exponentially),
Tf is "bounded."{The inverse of Tf is bounded from above
by exp[(E - Vrnin )/(2N - 1)kT], where Vmin is the global
minimum of the potential energy surface.} However, if
another thermostat is added to the dynamics, either to
control Tf or some subset of the degrees of freedom in the
system then the restriction imposed by the conserved quantity
is insufficient "to bound" the thermostat positions and
avoid regions of phase space that are (numerically) unstable
for the dynamics. In fact, the instability was first noticed
in attempts to numerically integrate a chainlike ansatz
within this formalism. The method is therefore not as
useful or flexible as the usual Nose-Hoover construction.
Note that though there is another form of the Winkler
method with slightly different equations of motion, it has
the same difficulties as the more natural variant discussed
above.
APPENDIX B
In order to determine reasonable values for the thermostat
masses, a second order equation of motion is generated
for each of the iJj from the time derivative of 1Jj
(Bl)
cfliJM-l
dr {
2iJM - 2 '2
-Q [QM-3TfM_3- kT]
M-I
TfM [ '2 k } '--Q QM-2TfM-2- T] -TfM-l
M-l
'2 1 '2 }
X QM-l TfM+ Q)QM-ITfM_I-kT] ,
d
2
iJM= {2iJM- 1 [Q '2 -kT]} _.
dr QM M-2TfM-2 TfM QM .
These equations can be solved individually, in the limit
Tfi is fast compared to Tfi-l and Tfi+l> while Tfi+2 moves on
the same time scale. This permits us to take functions of
Tfi-I and Tfi+l equal to their average values.3
The result is
cfliJl . [2NkT 2kT] Ql'3
dr = -Tfl Q2 - Q2 Tfl'
cfliJj= _iJj[2kT_ 2kT]_ Qj-l
--;[iZ Qj Qj+ 1 Qj+ 1 J
(B2)
d
2
iJM . [2kT]
---;p.-=-TfM QM .
The choices Ql = NkT/ ui and Qj = kT/ ui give thermostats
I to M - I an average "frequency" of (i). This frequency
is calculated by averaging the iJ; in the third order
term over the distribution function. The Mth thermostat
oscillates with frequency 2(i). The arbitrary parameter (i) is
chosen based on the properties of potential energy surface
(phonon frequencies, etc). Several approximations have
gone into the analysis and, in fact, the choice of mass itself
violates some of the approximations. This is only meant to
give a rough estimate.
APPENDIX C
In this Appendix, a proof that the trace of the stability
matrix is zero for chain dynamics is presented. The diagonal
elements of A, the stability matrix, are
Aii=-Sl> i=I,N,
(el)
where the s's evaluated at a fixed point. The trace of A is
then
M
-Tr(A)=Nst + L Si (e2)
j=2
From the fixed point equations, Eqs. (11), it is immediately
clear that SM-l = ± jc;.. Solving for SM in terms of
SM-2 in Eqs. (11) and substituting into the trace gives
M-2 1
-Tr(A)=Nst + Sj+ jc;.+ .[ci
M-2 1
=Ns1+ (e3)
J. Chern. Phys., Vol. 97, No.4, 15 August 1992
Downloaded 28 Dec 2011 to 128.122.251.131. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
Martyna, Klein, and Tuckerman: Nose-Hoover chains 2643
Similarly, solving for in terms of substituting
into the trace and simplifying gives
M-3 1 1
l:j=2 va ava
(C4)
At this point, the other end of the equation is simplified.
The variable is replaced by Nbl = (a + b2b3)/bl and
b2 is replaced by b2= -a1bf' Substituting into the trace
and simplifying gives
1 M-3 1
- Tr(A) = -- + l: bj-L bt--3
a j=3 Va
1
+ C (C5)
ava
Now, substituting and simplifying gives
1 M-3 1
- Tr(A) = -- + l: C
a j=4 Va
1
+-=-:---r (C6)
ava
Similarly, substituting and simplifying gives
1 M-3 1
-Tr(A)=--bib5+ l: bj-Lbt--3
a j=5 va
1
+ C (C7)
ava
Each subsequent substitution cancels the next term in the
sum until we get to M -4, where
1 2 1 2
a Va
1
ava
(C8)
Inserting the usual term = a + gives
12 12 1".4
-Tr(A)=--bM_3bM_2- CbM-3+ C!>M-3'
a Va aVa
(C9)
Finally, using the fact that bM-2 = (bM-3 - a)1 ragives
1 2 2 1 2
aVa va
1
+-:::-7':': - 3 =0,
aVa
which completes the proof.
(ClO)
The zero trace condition can also be used to show that
Louiville's theorem holds in the vicinity of the fixed points.
In this region, the linearized equations hold
x;= l: A;J'j, (ell)
j
where the x are a general set of variables (i.e., x represents
the p's, q's, and Now according to Liouville's theorem,
the condition for incompressible phase space volume
is
ax;£.. --0
;=1 ax;- .
Now differentiating x; with respect to Xi gives
N
L Au=O,
i=1
(C12)
(C13)
which is equal to zero if the matrix, A is traceless. Thus the
phase space volume is conserved in the region of validity of
the equations and the fixed points are neither attractors nor
repellors.
IS. Nose, J. Chern. Phys. 81,511 (1984).
2W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
3S. Nose, Prog. Theor. Phys. Supp. 103, 1 (1991).
4W. G. Hoover, Computational Statistical Mechanics (Elsenieer, New
York, 1991).
5S. Toxvaerd and O. H. Olson, Ber. Bunsenges. Phys. Chern. 94, 274
( 1990).
6H. C. Andersen, J. Chern. Phys. 72, 2384 (1980).
7D. Kusnezov, A. Bulgac, and W. Bauer, Ann. Phys. 204, 155 (1990).
sl. P. Hamilton, Phys. Rev. A 42, 7467 (1990).
9R. G. Winkler, Phys. Rev. A 45, 2250 (1992).
IOJ. Jellinek and R. S. Berry, Phys. Rev. A 40,2816 (1989).
11 J. Jellinek, J. Phys. Chern. 92, 3163 (1988).
12A. J. Lichtenberg and M. A. Lieberman, Regular and Stochastic Motion
(Springer, New York, 1983).
13 Michael Tabor, Chaos and Integrability in Nonlinear Dynamics (Wiley,
New York, 1989).
14y. I. Arnold, Mathematical Methods of Classical Mechanics (Springer,
New York, 1978).
15David Ruelle, Chaotic Evolution and Strange Attractors (Cambridge
University, New York, 1983).
16A. Wolf, in Chaos, edited by Arun Y. Holden (Princeton University,
Princeton, 1986).
17J. Cao and B. J. Berne, J. Chern. Phys. 91, 6359 (1989).
ISW. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, J.
Chern. Phys. 76, 637 (1982).
J. Chern. Phys., Vol. 97, No.4, 15 August 1992
Downloaded 28 Dec 2011 to 128.122.251.131. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions