192 Canonical ensemble
-3 -2 -1 0 1 2 3
x-3
-2
-1
0
1
2
3
p-3 -2 -1 0 1 2 3
x-3
-2
-1
0
1
2
3
p-4 -2 0 2 4
p0
0.1
0.2
0.3
0.4
0.5
f(
p)-4 -2 0 2 4
x0
0.1
0.2
0.3
0.4
0.5
f(
x)(a) (b)(c) (d)Fig. 4.10 Phase space and distribution functions obtained by integrating the Nos ́e–Hoover
equations ̇x=p/mp ̇=−mω^2 x−(pη/Q)p,η ̇=pη/Q,p ̇η =p^2 /m−kT for a harmonic
oscillator withm= 1,ω= 1,Q= 1,kT = 1,x(0) = 0,p(0) = 1,η(0) = 0,pη(0) = 1.
(a) shows the phase spacepvs.xindependent ofηandpη, (b) shows the phase space for
pη∈[−ǫ,ǫ], whereǫ= 0.001, and (c) and (d) show distributionsf(p) andf(x) obtained from
the simulation (solid line) compared with the correct canonical distributions (dashed line).
4.10 Nos ́e–Hoover chains
The reason for the failure of the Nos ́e–Hoover equations when more than one conserva-
tion law is obeyed by the system is that the equations of motion do notcontain a suffi-
cient number of variables in the extended phase space to offset therestrictions placed
on the accessible phase space caused by multiple conservation laws.Each conservation
law restricts the accessible phase space by one dimension. In orderto counteract this
effect, more phase space dimensions must be introduced, which canbe accomplished
by introducing additional variables. But how should these variables be added so as to
give the correct distribution in the physical phase space? The answer can be gleaned
from the fact that the momentum variablepηin the Nos ́e–Hoover equations must have
a Maxwell-Boltzmann distribution, just as the physical momenta do.In order to ensure
that such a distribution is generated,pηitself can be coupled to a Nos ́e–Hoover-type
thermostat, which will bring in a new set of variables, ̃ηandpη ̃. But once this is done,
we have the problem thatp ̃ηmust also have a Maxwell-Boltzmann distribution, which
requires introducing a thermostat for this variable. We could continue in this wayad
infinitum, but the procedure must terminate at some point. If we terminateit after
the addition ofMnew thermostat variable pairsηjandpηj,j= 1,...,M, then the