Nonequilibrium molecular dynamics
Nonequilibrium Molecular Dynamics is a powerful simulation tool. Like its equilibrium cousin, nonequilibrium molecular dynamics is based on timereversible equations of motion. But unlike conventional mechanics, nonequilibrium molecular dynamics provides a consistent microscopic basis for the irr...
Gespeichert in:
Datum: | 2005 |
---|---|
Hauptverfasser: | , |
Format: | Artikel |
Sprache: | English |
Veröffentlicht: |
Інститут фізики конденсованих систем НАН України
2005
|
Schriftenreihe: | Condensed Matter Physics |
Online Zugang: | http://dspace.nbuv.gov.ua/handle/123456789/119545 |
Tags: |
Tag hinzufügen
Keine Tags, Fügen Sie den ersten Tag hinzu!
|
Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
Zitieren: | Nonequilibrium molecular dynamics / Wm.G. Hoover, C.G. Hoover // Condensed Matter Physics. — 2005. — Т. 8, № 2(42). — С. 247–260. — Бібліогр.: 10 назв. — англ. |
Institution
Digital Library of Periodicals of National Academy of Sciences of Ukraineid |
irk-123456789-119545 |
---|---|
record_format |
dspace |
spelling |
irk-123456789-1195452017-06-08T03:03:25Z Nonequilibrium molecular dynamics Hoover, Wm.G. Hoover, C.G. Nonequilibrium Molecular Dynamics is a powerful simulation tool. Like its equilibrium cousin, nonequilibrium molecular dynamics is based on timereversible equations of motion. But unlike conventional mechanics, nonequilibrium molecular dynamics provides a consistent microscopic basis for the irreversible macroscopic Second Law of Thermodynamics. We recall here how fast computers led to the development of nonequilibrium molecular dynamics from the statistical mechanics of the 1950s. Computer-based theories facilitated revolutionary breakthroughs in understanding during the 1970s and 1980s. The new idea key to the nonequilibrium development was the replacement of the external thermodynamic environment by internal control variables. The new variables can control temperature, or pressure, or energy, or stress, or heat flux. These thermostat, barostat, ergostat, ... variables can control and maintain nonequilibrium states. We illustrate the methods with a simple example well-suited to student exploration, a thermostatted harmonic oscillator exposed to a temperature gradient. Нерівноважна молекулярна динаміка є потужним комп’ютерним методом. Подібно до рівноважної, нерівноважна молекулярна динаміка базується на часово-зворотніх рівняннях руху. Але, на відміну від звичайної динаміки, нерівноважна молекулярна динаміка містить узгоджений мікроскопічний базис для незворотнього макроскопічного другого закону термодинаміки. Ми показуємо тут, як швидкодіючі комп’ютери приводять до розвитку нерівноважної молекулярної динаміки на основі статистичної механіки 1950-х. Теорії, що базується на застосуванні комп’ютерів, сприяли революційному прориву в розумінні протягом 1970-х, 1980-х. Новою ключовою ідеєю для нерівноважного розвитку було заміщення зовнішнього термодинамічного середовища на внутрішні контролюючі змінні. Нові змінні можуть контролювати температуру або тиск, енергію, напруження або тепловий потік. Такі термостатичні, баростатичні, ергостатичні,... змінні можуть контролювати і підтримувати нерівноважні стани. Для ілюстрації ми використовуємо термостатичний гармонічний осцилятор, який піддається дії температурного градієнта. 2005 Article Nonequilibrium molecular dynamics / Wm.G. Hoover, C.G. Hoover // Condensed Matter Physics. — 2005. — Т. 8, № 2(42). — С. 247–260. — Бібліогр.: 10 назв. — англ. 1607-324X PACS: 02.70.NS, 04.25.-g, 05.10.-a, 05.70.Ln DOI:10.5488/CMP.8.2.247 http://dspace.nbuv.gov.ua/handle/123456789/119545 en Condensed Matter Physics Інститут фізики конденсованих систем НАН України |
institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
collection |
DSpace DC |
language |
English |
description |
Nonequilibrium Molecular Dynamics is a powerful simulation tool. Like its
equilibrium cousin, nonequilibrium molecular dynamics is based on timereversible
equations of motion. But unlike conventional mechanics, nonequilibrium
molecular dynamics provides a consistent microscopic basis for the
irreversible macroscopic Second Law of Thermodynamics. We recall here
how fast computers led to the development of nonequilibrium molecular dynamics
from the statistical mechanics of the 1950s. Computer-based theories
facilitated revolutionary breakthroughs in understanding during the
1970s and 1980s. The new idea key to the nonequilibrium development
was the replacement of the external thermodynamic environment by internal
control variables. The new variables can control temperature, or pressure,
or energy, or stress, or heat flux. These thermostat, barostat, ergostat,
... variables can control and maintain nonequilibrium states. We illustrate
the methods with a simple example well-suited to student exploration,
a thermostatted harmonic oscillator exposed to a temperature gradient. |
format |
Article |
author |
Hoover, Wm.G. Hoover, C.G. |
spellingShingle |
Hoover, Wm.G. Hoover, C.G. Nonequilibrium molecular dynamics Condensed Matter Physics |
author_facet |
Hoover, Wm.G. Hoover, C.G. |
author_sort |
Hoover, Wm.G. |
title |
Nonequilibrium molecular dynamics |
title_short |
Nonequilibrium molecular dynamics |
title_full |
Nonequilibrium molecular dynamics |
title_fullStr |
Nonequilibrium molecular dynamics |
title_full_unstemmed |
Nonequilibrium molecular dynamics |
title_sort |
nonequilibrium molecular dynamics |
publisher |
Інститут фізики конденсованих систем НАН України |
publishDate |
2005 |
url |
http://dspace.nbuv.gov.ua/handle/123456789/119545 |
citation_txt |
Nonequilibrium molecular dynamics / Wm.G. Hoover, C.G. Hoover // Condensed Matter Physics. — 2005. — Т. 8, № 2(42). — С. 247–260. — Бібліогр.: 10 назв. — англ. |
series |
Condensed Matter Physics |
work_keys_str_mv |
AT hooverwmg nonequilibriummoleculardynamics AT hoovercg nonequilibriummoleculardynamics |
first_indexed |
2025-07-08T16:08:56Z |
last_indexed |
2025-07-08T16:08:56Z |
_version_ |
1837095656382005248 |
fulltext |
Condensed Matter Physics, 2005, Vol. 8, No. 2(42), pp. 247–260
Nonequilibrium molecular dynamics
Wm.G.Hoover, C.G.Hoover
Highway Contract 60, Box 565, Ruby Valley, Nevada 89833,
Department of Applied Science University of California
at Davis/Livermore and Lawrence Livermore National Laboratory,
Livermore, California 94551–7808.
Received October 14, 2004
Nonequilibrium Molecular Dynamics is a powerful simulation tool. Like its
equilibrium cousin, nonequilibrium molecular dynamics is based on time-
reversible equations of motion. But unlike conventional mechanics, nonequi-
librium molecular dynamics provides a consistent microscopic basis for the
irreversible macroscopic Second Law of Thermodynamics. We recall here
how fast computers led to the development of nonequilibrium molecular dy-
namics from the statistical mechanics of the 1950s. Computer-based the-
ories facilitated revolutionary breakthroughs in understanding during the
1970s and 1980s. The new idea key to the nonequilibrium development
was the replacement of the external thermodynamic environment by inter-
nal control variables. The new variables can control temperature, or pres-
sure, or energy, or stress, or heat flux. These thermostat, barostat, ergo-
stat, ... variables can control and maintain nonequilibrium states. We illus-
trate the methods with a simple example well-suited to student exploration,
a thermostatted harmonic oscillator exposed to a temperature gradient.
Key words: second law, thermostats, chaos, fractals
PACS: 02.70.NS, 04.25.-g, 05.10.-a, 05.70.Ln
1. Introduction
David Barash’s invitation to celebrate Doug Henderson’s 70th Birthday in Utah
was welcome for many reasons. Besides the social and scientific values there was
the chance to visit a part of the West near our home-building project in Ruby
Valley, Nevada. The Interstate Highway System facilitated Bill’s visit to honor Doug
at Brigham Young University in Provo while Carol tended the horses, cats, and
computers in and around our homes in California.
Doug is an “old” friend in more ways than one. His work in equilibrium statistical
mechanics, in the 1960s, was along similar lines to Bill’s. There were plenty of
opportunities for getting together for discussions: Gordon conferences, the annual
West Coast Statistical Mechanics meetings, and social and scientific visits at our
various Bay Area homes and workplaces.
c© Wm.G.Hoover, C.G.Hoover 247
Wm.G.Hoover, C.G.Hoover
The present report summarizes our progress since then, progress in understandi-
ng nonequilibrium statistical mechanics following the development of nonequilibrium
molecular dynamics. The work described here spans roughly the interval between
Doug’s 40th and 70th birthdays. We wish him another such interval!
We sketch the main ideas here in an informal way. The detailed arguments and
literature references can be found in the references, supplemented, for the most
recent work, by a diligent internet search. While our supply lasts, we would be
happy to mail a copy of Bill’s book, “Time Reversibility, Computer Simulation, and
Chaos” [1], to any interested readers desiring more details. An earlier (1991) book
[2] “Computational Statistical Mechanics”, is available in pdf form at our website:
http://williamhoover.info.
2. Equilibrium statistical mechanics
Statistical mechanics undertakes the explanation of macroscopic thermodynamic
and hydrodynamic behavior in terms of underlying microscopic models. All the equi-
librium behavior (radial distribution functions, equations of state, phase diagrams,
...) of these models follows from Gibbs’ and Boltzmann’s 1883 idea of replacing time
averages by phase-space averages: a properly weighted (q, p) phase-space average
can provide the same information as a long-time-averaged dynamical simulation.
The “materials” treated by statistical mechanics are mechanical models, pre-
scriptions for the motion of microscopic masses. The simplest models are described
by a Hamiltonian H = K + Φ with a momentum-dependent kinetic energy K and
a coordinate-dependent potential Φ. Hamilton’s equations of motion describe the
corresponding dynamics in terms of first-order differential equations of motion.
H(q, p) = Φ(q) + K(p) −→ {q̇ = +∂H/∂p = ∂K/∂p = p/m} ;
{ṗ = −∂H/∂q = −∂Φ/∂q = F (q)} .
Numerical solutions of Hamilton’s equations can be obtained using the Runge-Kutta
method. In numerical work the product of the number of (q, p) degrees of freedom
and the number of timesteps is limited to about 1012.
A theoretical analysis of Hamilton’s equations has a fundamental consequence
for the equilibrium statistical mechanics of isolated systems. In this case the prob-
ability density f(q, p, t) is necessarily stationary at any fixed phase-space location:
∂f/∂t = 0, implying that the density must also be stationary following the flow,
df/dt = 0. Liouville’s Theorem (that a phase-space probability density governed
by Hamiltonian mechanics flows through the space with unchanging density, like an
incompressible fluid) implies that the stationary solution must necessarily have the
same constant value throughout all the mutually-accessible parts of phase space:
df/dt = 0 −→ f(q, p) = constant.
The time derivative following the flow, df/dt, is necessarily zero as a direct conse-
quence of the motion equations, leading to the conclusion that f has the same value
at any state accessible from the initial conditions placed on the flow.
248
Nonequilibrium molecular dynamics
The resulting constant-density distribution of states is Gibbs’ “microcanonical”
ensemble. There is an important special case. In the event that a part of the system
is an ideal-gas thermometer the equilibrium phase-space distribution is a product,
f(q, p) = fideal × fother .
In this special case, the entropies, S = −k ln f are additive, and entropy maximiza-
tion (the definition of the equilibrium situation) establishes that the temperatures,
Tideal = 〈p2/mk〉ideal , Tother = 〈p2/mk〉other
and pressures are equal too, and that the distribution fother is Gibbs’ “canonical”
distribution:
fother ∝ exp[−Hother(q, p)/kT ], T = Tideal = Tother .
The identification of the macroscopic thermodynamic concepts of entropy and tem-
perature with microscopic properties of a stationary distribution is the most signif-
icant result of Gibbs’ and Boltzmann’s work.
The canonical distribution then leads to all the usual textbook statistical me-
chanics, relating all the macroscopic equilibrium properties (including the Green-
Kubo expressions for the transport coefficients) to microscopic averages. The main
catch is that the microscopic expressions, integrals in 6N -dimensional (q, p) “phase
space” or 3N -dimensional (q) “coordinate space”, are just as hard to evaluate as are
many-body dynamical averages.
From the 1930s through the 1950s a widespread theoretical effort sought ways
to simplify the evaluation of microscopic phase-space averages. The Mayers’ virial
series provided density expansions around the low-density (ideal gas) limit. Einstein
models and cell models provided higher-density approximations. Integral equations
(for the pair distribution function) provided approximate approaches of uncertain
validity over a wide density range.
In the 1950s a new approach began to render these older ways obsolete. Fermi,
Alder, and Vineyard, along with many coworkers, developed molecular dynamics,
going back to the idea of time averages and using fast computers to solve the many-
body dynamics. Bill Wood and his coworkers at Los Alamos developed Metropolis’,
the Rosenbluths’, and the Tellers’ Monte Carlo method for phase-space averages.
Comparisons showed that the old Gibbs-Boltzmann approach was indeed correct.
Time averages and phase averages agreed.
Even though it was possible to solve equilibrium problems by direct simulati-
on(molecular dynamics and Monte Carlo were about equally costly) it was desirable
to save time/money by relating new problems to old solutions. “Perturbation the-
ory”, based on the computer results (free energy and radial distribution function)
for “reference systems”, allowed accurate estimates of thermodynamic properties.
Such theories were developed by several different groups. John Barker, Doug Hen-
derson, and George Stell were in the forefront of this effort. The result of their labors
was to finish off the problem of finding thermodynamic properties using equilibrium
249
Wm.G.Hoover, C.G.Hoover
statistical mechanics. A fringe benefit was the opening up of interest in nonequi-
librium problems, a subject which has kept us occupied for the past 30 years. The
high spots of this nonequilibrium effort are what we emphasize here in this paper
honoring Doug and the happy occasion of his 70th birthday.
3. Motion equations for nonequilibrium molecular dynamics
The equilibrium identification of macroscopic temperature with the microscopic
kinetic energy made it natural to constrain or control the temperature of a many-
body system by “rescaling” the kinetic energy. If the temperature was too low or
too high it could be adjusted by multiplying all the velocities by a common velocity-
rescaling multiplier. Though it was not realized at the time, a continuous version of
this rescaling procedure satisfies the simple differential equations of motion:
{mq̈ = ṗ = F (q) − ζp = F (q) − ζmq̇},
ζ =
∑
F · p/
∑
p2 =
∑
F · q̇/
∑
mq̇2,
−→ dK/dt =
∑
p · ṗ/m =
∑
p · (F − ζp)/m ≡ 0.
The “friction coefficient” (or “control variable”) ζ is positive whenever the forces
{F} are acting to increase the kinetic energy, and negative whenever the forces would
cause the kinetic energy to decrease. ζ keeps the kinetic energy (temperature) of the
system strictly constant.
Of course there are other possibilities for equations of motion controlling kinetic
energy. A constraint force varying as p3 or p5 could be used. Why should the “right”
constraint force be linear in the momenta? There is a long-standing principle in
mechanics, Gauss’ Principle of Least Constraint, stating that constraints should be
applied with the smallest possible (in the least-squared sense) constraint forces:
{ṗ = F + Fc},
∑
F 2
c /m a minimum.
When this least-constraint principle is applied to the problem of constraining the
kinetic temperature the solution,
{Fc = −ζp}, ζ =
∑
F · p/
∑
p2
is identical to the motion equations obtained above with velocity rescaling.
In a remarkable generalization of these simple ideas, Shuichi Nosé [3] showed
that a friction coefficient based on integral feedback:
ζ̇ =
∑
[(
p2/mkT
)
− 1
]
/#τ 2 →
ζ(t) = ζ(0) +
∫ t
0
∑
[(
p2/mkT
)
− 1
]
dt′/#τ 2 ,
was fully consistent with Gibbs’ canonical distribution for # thermostatted degrees
of freedom, each with mass m. Here τ is a relaxation time, the characteristic time
250
Nonequilibrium molecular dynamics
required to control temperature fluctuations. In this approach the kinetic energy
fluctuates, and in just such a way as to trace out the canonical distribution over
time. The friction coefficient consistent with Gibbs’ distribution has a Gaussian
distribution:
f(q, p, ζ) ∝ exp(−H/kT ) exp(−#ζ2τ 2/2).
Nosé based his development on Hamiltonian mechanics, through two special Hamil-
tonians which lead directly to the friction-coefficient equations of motion. This
Hamiltonian basis guarantees their time-reversible nature. Time reversibility im-
plies that any movie of a motion satisfying the equations will also satisfy them if
the frames are projected in reversed order. It is an interesting excercise to show that
the “Nosé-Hoover” equations mq̈ = F − ζmq̇ result from the better of Nosé’s two
Hamiltonians:
HNosé = (K/s) = s
[
Φ + (p2
s/2M) + #kT ln s
]
= 0,
where the “mass” M associated with the action variable ps is #kTτ 2.
It is somewhat simpler to develop Nosé’s ideas (as well as their generalizations to
constant-stress or constant-strain ensembles) by (i) assuming a constraint force linear
in the momenta, FC = −ζp and then (ii) using the corresponding form of Liouville’s
Theorem [4] to evaluate the change of f in time for the resulting compressible phase-
space flow:
∂f/∂t =
∑
[−∂(f q̇)/∂q − ∂(fṗ)/∂p] − ∂(f ζ̇)/∂ζ,
d ln f/dt =
∑
ζ.
Adding the further requirement (iii) that the resulting f have Gibbs’ equilibrium
form results again in the “Nosé-Hoover” form of Nosé’s feedback equation:
ζ̇ =
∑
[
(p2/mkT ) − 1
]
/#τ 2.
The new motion equations, with their control variables, made it possible to
constrain temperature or pressure or energy for selected degrees of freedom. The
constraints could also be generalized to apply locally, in either space or time. In the
equilibrium case, with temperature and pressure uniform, all the results of statistical
mechanics were reproduced and recovered. We turn next to the consequences of the
new dynamics away from equilibrium.
4. Results from nonequilibrium molecular dynamics
Many-body diffusive flows, viscous flows, and heat flows were all simulated in
the 1970s. By 1990 million-particle systems had been simulated. A sample appears
in figure 1. In that work we simulated the plastic indentation of silicon (using the
Stillinger-Weber many-body potential) and found a yield stress in good agreement
with experimental data for silicon [5].
251
Wm.G.Hoover, C.G.Hoover
Figure 1. Indentation of silicon (as represented by the Stillinger-Weber potential)
using thermostatted nonequilibrium molecular dynamics. Results were compared
using two different types of tetrahedral indentors, one that was perfectly smooth,
the other atomistic.
Few-particle systems, small enough for a thorough analysis of their phase-space
distributions, were simulated too. This work was stimulated by Nosé’s breakthrough
of 1984 [3]. The phase-space structures the small-system simulations revealed were all
multifractal, with the phase-space probability density singular everywhere, varying
as a fractional power of separation. These wild nonequilibrium distributions are
qualitatively unlike the smooth distributions known to Gibbs, Boltzmann, and to
the generations of students of equilibrium statistical mechanics brought up on their
ideas.
Consider the simple Galton-Board problem shown in figure 2. A mass point par-
ticle falls (under the effect of a fixed gravitational field) through a triangular-lattice
array of scatterers. The dynamics is constrained to occur at fixed kinetic energy, by
using Gauss’ Principle. A sample nonequilibrium distribution, describing the colli-
sions which occur in a nonequilibrium Galton Board, is also shown in figure 2 [6].
Each point in the figure indicates the location of a collision (on the abscissa) and the
tangential velocity (the ordinate) for a mass point scattering from the array of hard
scatterers. Evidently the distribution is fractal rather than continuous. The density
252
Nonequilibrium molecular dynamics
Figure 2. A finite portion of the Galton Board scattering system is shown at top
right. Below is shown the phase plane for the Galton Board problem. The 106
attractor points (α, sin β) shown here represent one million successive collisions of
a point particle in an infinite periodic array of hard-disk scatterers. This structure
is a fractal “strange attractor”. Its time-reversed cousin, the unobservable “even
stranger repellor”, is shown to the right. For details see [6].
of points varies as a fractional power of the separation everywhere in the distribu-
tion. At equilibrium, where the field vanishes and Gauss’ constraint is unnecessary,
the corresponding distribution of collisions would instead cover the plot smoothly
and uniformly.
At first, the time-reversibility properties of these fractal distributions seemed
paradoxical. For every trajectory going forward in time a reversed trajectory, with
the coordinate frames played backwards, was also a solution of the motion equa-
tions. How could this symmetric situation be consistent with the Second Law of
Thermodynamics, which restricts observable time developments to those for which
entropy increases? We explore the details of the explanation for a simple mechanical
model in section 5.
The resolution of (Loschmidt’s) Reversibility Paradox using Nosé-Hoover me-
chanics is simplest for nonequilibrium steady states [1,7] such as the Galton Board
problem discussed above. The phase-space flow in such cases invariably seeks out a
fractal strange attractor. The dynamics on this attractor is (Lyapunov) unstable in
the sense that two nearby points tend to separate, exponentially fast, in time. The
attractor dynamics is stable in the sense that the flow contracts (on the average),
eventually occupying zero phase volume! The corresponding phase-space density
f(q, p) is singular everywhere! The time-reversed version of the attractor dynamics
(which defines the “repellor”) has also reversed stability properties. In the neigh-
253
Wm.G.Hoover, C.G.Hoover
borhood of the repellor the time-averaged flow expands rather than contracts (on
the average) and is, as a consequence, too unstable to be observed. In fact, the
only way in which the unobservable phase-space repellors can be generated is by
time-reversing the observed attractor states.
5. Example problem: nonequilibrium heat
conducting oscillator
The concepts of fractal distributions, time-reversible irreversibility, and Lya-
punov instability are all best assimilated through simple examples. The harmonic
oscillator is arguably the simplest system with a smooth phase-space distribution.
When Nosé-Hoover mechanics is applied to it the motion takes place in a three-
dimensional (q, p, ζ) phase space:
q̇ = p/m, ṗ = −κq − ζp, ζ̇ =
[
(p2/mkT ) − 1
]
/τ 2.
These flow equations are consistent with Gibbs’ canonical distribution for the osci-
llator:
f(q, p, ζ) ∝ fq × fp × fζ ,
fq = exp
(
−κq2
2kT
)(
2πkT
κ
)
−1/2
,
fp = exp
(
−p2
2mkT
)
(2πmkT )−1/2 ,
fζ = exp
(
−ζ2τ 2
ζ
2
)
(
2π
τ 2
ζ
)
−1/2
.
Numerical investigation [8,9] showed that this system is not sufficiently mixing
in phase space to access the full canonical distribution. The oscillator phase-space
is broken up into infinitely-many noncommunicating regions. There is a “chaotic”
region, in which nearby trajectories depart exponentially fast from one another.
This exponential rate of trajectory separation is “Lyapunov instability”, the defining
feature of chaos. In addition to the unstable region there are an infinite number of
noncommunicating stable regions. The lack of ergodicity in the model suggests that
the Nosé-Hoover oscillator is not a good model for equilibrium statistical mechanics.
The results obtained from it would depend upon the initial conditions.
On the other hand, adding a second control variable, to control also the fourth
moment of the velocity distribution does achieve the full distribution [10]. The new
set of differential equations,
q̇ = p/m, ṗ = −κq − ζp − ξp(p2/mkT ),
ζ̇ = [(p2/mkT ) − 1]/τ 2
ζ ,
ξ̇ = [(p2/mkT )2 − 3(p2/mkT )]/τ 2
ξ
254
Nonequilibrium molecular dynamics
can readily be integrated using the fourth-order Runge-Kutta method. With a
timestep of 0.001 the integration error that results is of the same order as the
roundoff error in double-precision arithmentic. The distribution which results from
these equilibrium equations of motion is again an extended canonical distribution:
f(q, p, ζ, ξ) = fq × fp × fζ × fξ ,
where the ξ distribution is also Gaussian:
fξ = exp
(
−ξ2τ 2
ξ /2
) (
2π/τ 2
ξ
)
−1/2
.
Figure 3 shows the time development of the distribution in the (q, p) subspace
of the phase space for the special case
m = k = T = 1, τζ = τξ = 3.7.
The relaxation times are a little more than half the unconstrained oscillator period
of 2π. Numerical investigations suggest that this oscillator accesses the entire phase
space. The Lyapunov spectrum of exponents is
{λeq} = {+0.102, 0.000, 0.000,−0.102}.
In this equilibrium case there is no strange attractor. The full phase space is explored
by the motion, with the Gaussian distribution seen in figure 3.
Figure 3. Development of the {q, p} distribution function at equilibrium, for the
doubly-thermostatted oscillator. The final distribution is a Gaussian in all four
variables, with the rms values of ζ and ξ smaller than those of q and p by a factor
of τ = 3.7 Each of the three plots shows 200,000 points. The time separations
between successive points are 0.01, 0.10, and 1.0, where the Runge-Kutta timestep
is 0.001.
A nonequilibrium version of the same oscillator can be achieved by making the
temperature T a function of the coordinate q:
T
T0
= 1 + ε tanh
( q
h
)
.
255
Wm.G.Hoover, C.G.Hoover
In the special case
T0 = ε = h = 1, τζ = τξ = 3.7,
the motion is the limit cycle shown in figure 4. The dimensionality of this limit-cycle
attractor is unity, rather than four. The Lyapunov exponents for the limit cycle are
{λ3.7,3.7} = {0.000,−0.018,−0.079,−0.122}.
The lack of any positive Lyapunov exponents indicates that the one-dimensional
limit-cycle motion is completely stable.
Figure 4. Limit cycle for the nonequilibrium oscillator with both thermostat
variables set equal to 3.7. The reversed cycle is also shown. In both cases half a
million points, separated by one thousand timesteps, are shown. (q, p) and (ζ, ξ)
pairs of points are shown separately. The (ζ, ξ) pairs are those closer to the origin.
Figure 5. Strange attractor for the nonequilibrium oscillator with thermostat
relaxation times (τζ , τξ) = (3.8, 3.6). The reversed attractor (the repellor) is also
shown. In both cases half a million points, separated by one thousand timesteps,
are shown.
The sum total of the Lyapunov exponents, −0.219, is the time-averaged rate-of-
change of the comoving phase volume ⊗ in the neighborhood of the limit cycle:
〈
d ln⊗
dt
〉
= −
〈
d ln f
dt
〉
=
∑
λ = −0.219.
256
Nonequilibrium molecular dynamics
The reversed motion, shown also in figure 4, has a reversed Lyapunov spectrum
(with three positive exponents and no negative ones), and, though this motion is
formally a solution of the motion equations, it is unobservable.
A slight change in the two thermostat relaxation times to
τζ = 3.8, τξ = 3.6,
gives rise to the “strange attractor” shown in figure 5. The dimensionality of the
attractor can be estimated from the Lyapunov spectrum
{λ3.8,3.6} = {+0.010, 0.000,−0.056,−0.159}.
by using Kaplan and Yorke’s conjecture that the dimension corresponds to the
number of Lyapunov exponents whose sum is zero. Interpolating between the two-
exponent sum and the three-exponent sum:
λ1 + λ2 = 0.010, λ1 + λ2 + λ3 = −0.046,
gives a dimensionality of DKY = 2 + 10/56 = 2.18.
Figure 6. Histograms showing the densities of the time-averaged heat transfer
〈dQ/dt〉 (lighter curve) and entropy production 〈dSint/dt〉 (heavier curve) as
functions of the oscillator coordinate q for the chaotic oscillator away from equi-
librium. The total heat transfer is exactly zero while the entropy production is
necessarily positive (for microscopic stability and for consistency with Clausius’
form of the Second Law of Thermodynamics). Both total heat and total entropy
production are obtained by integrating the densities shown with respect to the
coordinate q.
In this model system the overall time-averaged heat transfer is zero, but, on
the average, the heat is taken in at a higher temperature than that at which heat is
extracted, corresponding to an entropy increase dSext/dt in the heat reservoirs inter-
acting with the oscillator. Figure 6 shows the time-averaged distributions of 〈dQ/dt〉
257
Wm.G.Hoover, C.G.Hoover
and 〈dSint/dt〉 in space corresponding to the 2.18-dimensional strange attractor. The
distributions for the limit cycle are scarcely different.
For either of the models, the total entropy change is
−
〈
d(Sext/k)
dt
〉
=
〈
d(Sint/k)
dt
〉
=
〈
(dQ/dt)
kT
〉
=
〈[
−ζp − ξ
p2
mkT
p
]
q̇
T
〉
=
〈
−ζ − 3ξ
p2
mkT
〉
.
The last equality follows from two time averages:
〈ζζ̇〉 ≡ 0 →
〈
ζ
p2
mkT
〉
= 〈ζ〉,
〈ξξ̇〉 ≡ 0 →
〈
ξ
(
p2
mkT
)2
〉
=
〈
3ξ
(
p2
mkT
)〉
.
Notice that the time-averaged expression
〈
d(Sext/k)
dt
〉
=
〈
ζ + 3ξ
p2
mkT
〉
is identically equal to the time-averaged rate at which the phase volume collapses
onto the strange attractor:
−
〈
d ln⊗
dt
〉
=
〈
d ln f
dt
〉
= −
〈
∂q̇
∂q
〉
−
〈
∂ṗ
∂p
〉
−
〈
∂ζ̇
∂ζ
〉
−
〈
∂ξ̇
∂ξ
〉
= 0 +
〈
ζ + 3ξ
p2
mkT
〉
+ 0 + 0
−→
〈
d(Sext/k)
dt
〉
= −
〈
d(Sint/k)
dt
〉
=
〈
d ln f
dt
〉
=
〈
ζ + 3ξ
p2
mkT
〉
= +0.205.
This correspondence is no coincidence. Here, as well as in the simpler Nosé-Hoover
case, the time-rate-of-change of the phase volume (equivalent to the time-rate-of-
change of the internal Gibbs’ entropy in equilibrium statatistical mechanics) is ex-
actly equal to the time-averaged rate of the external entropy increase in the heat
reservoirs represented by the control variables ζ and ξ.
The Second Law statement that the overall entropy change in the external reser-
voirs must be positive is precisely equivalent to the microscopic statement that phase
volume must collapse (onto an attractor) rather than diverge. Recall Clausius’ form
for the Second Law of Thermodynamics,
∮
1
T
dQ
dt
dt′ < 0,
258
Nonequilibrium molecular dynamics
where the integral is over any cyclic process which returns the system to its initial
state, (dQ/dt) is the rate at which heat is taken in by the system during the pro-
cess, and T is the temperature of the reservoir with which the heat transfer takes
place at the time t′. We see that the microscopic Second Law of Thermodynamics,
from nonequilibrium molecular dynamics, has exactly Clausius’ form, but with an
additional averaging operation:
〈
∮
1
kT
dQ
dt
〉
< 0.
In the microscopic case the average is over (many repetitions of) a nonequilibrium
cyclic process. The requirement that the overall entropy increase (of the reservoirs!)
be positive, is a direct consequence of the microscopic requirement that the occupied
phase space cannot diverge in any stable microscopic nonequilibrium flow cycle.
Either by consciously violating Gauss’ Principle of Least Constraint or by making
an otherwise unwise choice of thermostat variables, the equivalence of the rate of
microscopic phase-volume collapse with the rate of thermodynamic entropy increase,
can be lost. For examples, see the references listed on Rainer Klages’ website. It is
evident that thermostats and nonequilibrium motion equations which maintain the
equivalence of the dynamics with the Second Law are best. The zero-volume nature
of the nonequilibrium attractors required for this stability also shows just how rare
the nonequilibrium states are.
6. Conclusions
The same dynamical tools which solved the equilibrium many-body problem
of statistical mechanics have been generalized to make possible the simulation of
nonequilibrium problems. Velocity rescaling, Gauss’ Principle of Least Constraint,
and nonequilibrium Nosé-Hoover mechanics all lead to exactly the same form for
the nonequilibrium motion equations,
{ṗ = F − ζp}.
Nosé-Hoover mechanics (and the generalization of it in the example detailed above)
has also an intimate connection with the Second Law of Thermodynamics. Anal-
yses of the strange attractors generated by that mechanics have shown that the
irreversible Second Law has its roots in the chaotic nature of the time-reversible
strange attractors generated by constrained flows. The microscopic version of the
Second Law is the (time-averaged) statement that phase volume away from equilib-
rium must shrink (and ultimately vanish) rather than grow (and ultimately diverge).
Despite this new understanding there is considerable work yet to be done in char-
acterizing the attractors. These are ideal problems for undergraduate and graduate
research. The dimensionality of these small-system phase spaces is just right for the
numerical investigation of Lyapunov Spectra, the various fractal dimensionalities,
and conjectures like that proposed so long ago by Kaplan and Yorke.
259
Wm.G.Hoover, C.G.Hoover
This work was carried out under the auspices of the United States Department of
Energy at the Lawrence Livermore National Laboratory under Contract W–7405–
Eng–48.
References
1. Hoover Wm.G. Time Reversibility, Computer Simulation, and Chaos. World Scientific,
Singapore, 1999 and 2001.
2. Hoover W.G. Computational Statistical Mechanics. Elsevier, New York, 1991.
3. Nosé S., J. Chem. Phys., 1984, 81, 511.
4. Hoover Wm.G., J. Chem. Phys., 1998, 109, 4164.
5. Kallman J.S., Hoover W.G., Hoover C.G., De Groot A.J., Lee S.M., Wooten F., Phys.
Rev. B, 1993, 47, 7705.
6. Moran B. , Hoover W.G., Bestiale S., J. Stat. Phys., 1987, 48, 709.
7. Holian B.L., Hoover W.G., Posch H.A., Phys. Rev. Letts., 1987, 59, 10.
8. Hoover W.G., Phys. Rev. A, 1985, 31, 1685.
9. Posch H.A., Hoover W.G., Vesely F.J., Phys. Rev. A, 1986, 33, 4253.
10. Hoover Wm.G., Holian B.L., Phys. Lett. A, 1996, 211, 253.
Нерівноважна молекулярна динаміка
У.Г.Гувер, К.Г.Гувер
Університет Каліфорнії, Девіс, Каліфорнія, США
Отримано 14 жовтня 2004 р.
Нерівноважна молекулярна динаміка є потужним комп’ютерним
методом. Подібно до рівноважної, нерівноважна молекулярна дина-
міка базується на часово-зворотніх рівняннях руху. Але, на відміну
від звичайної динаміки, нерівноважна молекулярна динаміка мі-
стить узгоджений мікроскопічний базис для незворотнього ма-
кроскопічного другого закону термодинаміки. Ми показуємо тут,
як швидкодіючі комп’ютери приводять до розвитку нерівноважної
молекулярної динаміки на основі статистичної механіки 1950-
х. Теорії, що базується на застосуванні комп’ютерів, сприяли
революційному прориву в розумінні протягом 1970-х, 1980-х. Новою
ключовою ідеєю для нерівноважного розвитку було заміщення зов-
нішнього термодинамічного середовища на внутрішні контролю-
ючі змінні. Нові змінні можуть контролювати температуру або тиск,
енергію, напруження або тепловий потік. Такі термостатичні, баро-
статичні, ергостатичні,... змінні можуть контролювати і підтримувати
нерівноважні стани. Для ілюстрації ми використовуємо термоста-
тичний гармонічний осцилятор, який піддається дії температурного
градієнта.
Ключові слова: другий закон, термостати, хаос, фрактали
PACS: 02.70.NS, 04.25.-g, 05.10.-a, 05.70.Ln
260
|