PHYSICAL REVIEW 8VOLUME 43, NUMBER 41FEBRUARY 1991
Tempe a u e and p essu e cons ain s nea he eezing poin
Juan J.Mo ales and Ma iaaia J.Nue o
Depa amen o de FIsica, Facul ad de Ciencias, Uni e sidad de Ex emadu a, 06011Badajoz, Spain
Luis F.Rull
Depa amen o de Fssica Teo ica, Uni e sidad de Se illa, 41080Se illa, Spain
(Recei ed 26 Feb ua y 1990; e ised manusc ip ecei ed 3May 1990)
The iso he mal-isoba ic ensemble molecula -dynamics me hod MD(T,p, N) p oposed by Nose
and Hoo e is used o s udy he Auc ua ions in a wo-dimensional Lenna d-Jones Auid, close o he
eezing poin . The Tand pcons ain s in his me hod do no a ec he dynamical beha io o he
sys em, since spon aneous Auc ua ions in he densi y allow he sys em o eeze and mel jus as do
he Tand pAuc ua ions in he mic ocanonical ensemble MD(E, V, N) close o he mel ing zone.
I. INTRODUCTION
The dynamics o a omic sys ems wi h con inuous in-
e ac ions be ween he a oms equi es nume ical in eg a-
ion o Hamil on's equa ion o mo ion, '
dH
dpi
d ;
=- '= =P
d 'I
dp;
d PI= dU( )
d ,
whe e M=T(p)+ U( ) is he Hamil onian o he X-body
sys em, T(p )and U( )being he kine ic and po en ial
ene gies, espec i ely; New onian mechanics implies ha
he ene gy and momen um a e he conse ed a iables o
he mo ion. In adi ional molecula -dynamics (MD) ex-
pe imen s, he o al ene gy E, he numbe o pa icles N,
and he olume Va e conse ed as he dynamics o he
sys ems e ol es in ime. The ime a e age o any p ope -
y is an app oxima e measu e o he mic ocanonical en-
semble MD(E, V, X). Fo ce ain applica ions, i may be
desi able o pe o m dynamical simula ions a cons an
empe a u e and/o p essu e leading o he canonical en-
semble MD( T, V, N) and he iso he mal-isoba ic ensemble
MD(T,p, X) Howe e so. me o hese echniques need he
in oduc ion o a i icial p ocesses in o de o make his
p ac ical. Such schemes may be classi ied as non-
New onian MD, since hey a e no based on Hamil oni-
ans o mo ion o eal sys ems, bu use a i ices o scaling
ha esul in ic i ious o ces in o de o sa is y some
adop ed de ini ion o cons an Tand p. Howe e , i is
possible o show igo ously ha ce ain o hese me hods
gi e p ope equilib ium ensemble a e ages.
Se e al me hods ha e been de eloped o simula e he
canonical ensemble, including he "s ochas ic" me hod
in which collision wi h an imagina y hea -ba h pa icle is
aken in o accoun , he "ex ended-sys em" me hod" in
which adeg ee o eedom is included which ep esen s
he ese oi in con ac wi h he sys em, and he "con-
s ain s" me hod in which he eloci ies a e escaled a
each ime s ep by some ac o in o de o ix he kine ic
ene gy o he sys em. Al hough o he me hods do no
gene a e s a es in he canonical ensemble, hey seem e y
use ul o s udying he changing s a e and sys ems each-
ing equilib ium a adi e en empe a u e. Fo he
cons an -p essu e MD simula ion, in all app oaches i is
ine i able ha he sys em box mus change olume as in
cons an -p essu e Mon e Ca lo (MC) simula ions. The
nomencla u e is he same as o he MD(T, V,N). The
s ochas ic me hod is aMC calcula ion in which he
Me opolis echnique is applied o he andomly
changed olume AV. The ex ended me hod in ol es
coupling he sys em o an ex e nal a iable V, he olume
o he sys em, whe e he coupling mimics he ac ion o a
pis on on a eal sys em. The cons ain s me hod makes
he ins an aneous p essu e acons an o he mo ion. In
wo o he pa icula me hods, one allows he simula ion
box o change shape as well as size, and in he o he he
app op ia e ensemble is no well iden i ied.
In a o me s udy' we p esen ed some MD simula ions
o asys em e y close o he mel ing zone whe e Auc ua-
ions play an impo an ole in he dynamic p ope ies o
he sys em. The simula ions we e pe o med using he
adi ional mic ocanonical ensemble, i.e.,wi h no scaling
p ocedu e, and we ound luc ua ions o empe a u e
and p essu e which spon aneously make he sys em oscil-
la e be ween wo poin s o s a e. We conclude ha simu-
la ions wi h any cons ain s me hod, such as he
cons an - empe a u e and/o cons an -p essu e me hods,
a e inco ec nea he mel ing zone because he e a e
p oblems ela ed o luc ua ions wi h espec o space and
ime and he cons ain s can diminish he luc ua ions o
no allow hem a all. Ne e heless, we now see his con-
clusion o be no comple ely gene al, because he sys em
beha io depends on he me hod used and consequen ly
on he way in which he cons ain s a e applied.
Among he echniques e iewed he e we chose he
combina ion o he ex ended sys em me hod due o Nose
and he cons ain s me hod due o Hoo e , "leading o
Nose-Hoo e (NH) o mula ion. 'The e we e wo
easons o he choice: i s , because he NH scaling p o-
cedu e ep oduces bo h he canonical and he
iso he mal-isoba ic p obabili y densi y in he phase space
43 3514 QC1991 The Ame ican Physical Socie y
43 TEMPERATURE AND PRESSURE CONSTRAINTS NEAR THE. . . 3515
II. SIMULATIONS AND RESULTS
The equa ions o mo ion o he MD(T,p, N) ensemble
in he NH o mula ion a e '
dx; p;
m "D (3)
dp) = ;—
(g —
E)p;, (4)
whe e x;= , /V' a e he educed coo dina es o he
pa icle i, Dis he dimension o he space, , is he o ce,
and gand e=V/DV a e he iso he mal and isoba ic ic-
ion coe icien s which couple he sys em o a ese oi a
T=T„and p=p„, espec i ely.
The dynamical equa ions o he ic ion coe icien s
a e
N2
DNk T—
d ,.&I(DNkT )(5)
o aclassical mechanical sys em, and second, because his
me hod had been used p e iously o show ha he
dynamical beha io o he NH equa ions does no des oy
he long- ange bond co ela ion in he Auid nea he
eezing poin . 'In he ollowing we shall s udy he be-
ha io o he empe a u e, p essu e, and densi y Auc ua-
ions in he MD( T,p, N) ensemble in he wo-dimensional
(2D) sys em nea he eezing poin and he esul s will be
compa ed wi h hose ob ained by MD(E, V, N) nea he
mel ing poin . '
and
=(p —
p,„)V/(kT ),
d (6)
whe e ~z and ~a e he iso he mal and isoba ic elaxa-
ion imes o he sys em espec i ely which gi e he cou-
pling s eng hs. Bo h elaxa ion imes ha e aGaussian
dis ibu ion in he NH heo y and hey a e ela ed o he
mean collision ime. ''
The NH iso he mal-isoba ic simula ions we e pe -
o med o a20 Lenna d-Jones sys em o X=400 pa i-
cles wi h acu oA dis ance o , =2.5o..The s a e poin
was se up a kT/e=0. 7and po =0.83. This poin is a
Auid e y close o he eezing poin ' wi h aco espond-
ing p essu e o pe /E =2 64. .These alues o empe a-
u e and p essu e we e chosen as he ixed ex e nal Tand
p,„o he ese oi in Eqs. (5) and (6). The o he inpu
alues we e an=0. 005 and *= /&N =0'.3.'To in-
eg a e he equa ions o mo ion he Tox ae d algo i hm'
was used wi h a ime s ep h=0.005(mo /E)', and he
sys em was allowed o e ol e in ime abou 120000 h ak-
ing a e ages alues e e y 800 h.
In Fig. 1 he Auc ua ion o empe a u e is plo ed
agains he numbe o ime s eps o he MD(Tp, N) en-
semble. The luc ua ions a e e y small h oughou he
dynamic e olu ion o he sys em, in spi e o he icini y
o he eezing zone whe e he Auc ua ions a e la ge.
This means ha he iso he mal ic ion coe icien wo ks
well a in e changing ene gy be ween sys em and ese -
oi h ough Eq. (5), diminishing he luc ua ions a ound
he ixed empe a u e as much as possible, as was o be
expec ed in p inciple. The inal mean alue o he em-
pe a u e was kT/c, =0.7000+0.0007. These Auc ua ions
0.75
4J
CC
l—
UJ
CL
LU
LU 0.65
50 100 x800
NUMBER 0F TINE STEPS
FIG. 1. Reduced- empe a u e luc ua ions s numbe o ime s eps o a wo-dimensional luid close o he eezing zone.
3516 JUAN J.MORALES, MARIA J.NUEVO, AND LUIS F.RULL 43
2.70,
2.64
2.58
100 x800
NUNBER OF TINE STEPS
FIG. 2. As Fig. 1, bu o he educed p essu e.
a e abou 14 imes smalle han hose o he co espond-
ing MD(E, V,X) simula ion in he mel ing zone. '
Simila beha io is ound o he p essu e Auc ua ions
whe e he dynamical equa ion o he isoba ic ic ion
coe icien Eq. (6) main ains he sys em a cons an p es-
su e while allowing i o change i s olume. These Auc-
ua ions look mo e egula ha he empe a u e luc ua-
ions and g ea e by abou 3 imes. The inal mean alue
o he p essu e was po. / =2. 641+0.002, wi h Iuc ua-
ions 40 imes smalle han he p essu e Auc ua ion o
he MD(E, V, N), as can be seen by compa ison o he
s anda d de ia ion and Fig. 2wi h he esul s gi en in
Re . 18.
The ime e olu ion o he densi y is e y di6'e en , as
Fig. 3shows. The luc ua ion is i egula and la ge due
o he absence o cons ain s. On wo occasions, a e
l—
LLJ
0.87
0,82
0.77
50 j.00 x800
NUNBER OF TINE STEPS
FICx. 3. As Fig. 1, bu o he educed densi y.
43 TEMPERATURE AND PRESSURE CONSTRAINTS NEAR THE. ..3517
abou 44000 and 104000 ime s eps, luc ua ions appea
spon aneously which allow he sys em o each densi ies
wi hin he eezing zone. 'The inal mean alue o he
densi y was po. =0.827+0.010, his s anda d de ia ion
being he same as ha co esponding o he empe a u e
in he MD(E, V, N) ensemble a mel ing. These la ge
spon aneous luc ua ions disappea when he sys em is a
away om he eezing zone: simula ion in he liquid sys-
em a po. =0.7606 showed as anda d de ia ion o
+0.0027, luc ua ions nea ly ou imes smalle han be-
o e.
III. DISCUSSION AND CONCLUSIONS
I we compa e hese esul s MD(T,p, N) nea he eez-
ing poin wi h hose MD(E, V, N) co esponding o he
solid nea he mel ing zone, 'we can see ha he wo sys-
ems beha e in he same way physically. The mic o-
canonical ensemble allows empe a u e and p essu e o
Auc u a e in he same way as he densi y in he
iso he mal-isoba ic ensemble. Nea aphase ansi ion,
bo h ensembles mus unde go spon aneous luc ua ions
ha d i e he sys em back and o h be ween wo poin s
o s a e. The p oblem o he non-appea ance o he Auc-
ua ions a ises when he me hod used o simula e he
MD(T,p, N) is inco ec , e en hough i yields ime a e -
ages o he modynamically in e es ing a iables o
wi hin he s a is ical accu acy o mos expe imen al
da a. A e y in e es ing discussion can be ound in Re s.
7and 16, whe e he alidi y o he s ochas ic me hod
used o he olume change in Re . 7 o in es iga e he
mel ing ansi ion is e u ed in Re . 16, since he e o ob-
ained wi h ha me hod no only changes he co ec lo-
ca ion o he ansi ion bu also mis ea s any densi y
luc ua ions in i s icini y.
We can he e o e conclude ha he cons ain s o
empe a u e and p essu e in he ex ended Nose and Ho-
o e me hod no only do no al e he na u al e olu ion
o he sys em, bu also gi e he co ec dynamical beha -
io including he zone o la ge luc ua ions such as he
eezing ansi ion.
Du ing he p epa a ion o he p esen pape me hods
o implemen ing empe a u e and p essu e con ols a
equilib ium and nonequilib ium in he NH dynamics
ha e been published. 'These me hods inc ease he
e iciency o la ge-seal: simula ions bu he de ailed e-
sul s a e no ye a ailable.
ACKNOWLEDGMENT
The coope a ion o he Uni e si y o Ex emadu a
Compu e Cen e is g a e ully acknowledged.
See, o example, H. Golds ein, Classical Mechanics (Addison-
Wesley, Reading MA, 1980), o M. P. Allen and D. J.Tildes-
ley, Compu e Simula ion o Liquids (Cla endon, Ox o d,
1989}.
F.F.Ab aham, Ad . Phys. 35, 1(1985).
H. C. Ande sen, J.Chem. Phys. 72, 2384 (1980).
4S. Nose, Mol. Phys. 52, 255 (1984).
5D. J. E ans and G. P. Mo is, Compu . Phys. Rep. 1, 297
(1984).
H. J. C. Be endsen, J. P. M. Pos ma, W. F. Van Gus e en, A.
Di Nola, and J.R. Haak, J.Chem. Phys. 81, 3684 (1984).
7S. W. Koch and F.F.Ab aham, Phys. Re . B27, 2964 (1983).
N. Me opolis, A. W. Rosenblu h, M. N. Rosenblu h, A. H.
Telle , and E.Telle , J.Chem. Phys. 21, 1087 (1953).
J.R. Ray, J.Chem. Phys. 79, 5128 (1983).
L. F. Rull, J. J. Mo ales, and F. Cuad os, Phys. Re . B32,
6050 (1985).
W. G. Hoo e , Phys. Re . A31, 1695 (1985).
S.Nose, Mol. Phys. 57, 187 (1986).
'J. J. Mo ales, S. Tox ae d, and L. F. Rull, Phys. Re . A34,
1495 (1986).
~4D. J.E ans and B.L. Holian, J.Chem. Phys. 83, 4069 (1985).
~5L. F. Rull, J. J. Mo ales, and S. Tox ae d, Phys. Re . A38,
4309 (1988}.
S.Tox ae d, Phys. Re . B29, 2821 (1984).
S.Tox ae d, J.Compu . Phys. 47, 444 (1982}.
8J. J. Mo ales, F. Cuad os, and L. F. Rull, J. Chem. Phys. S6,
2960 (1987).
B. L. Holian, A. J. De G oo , W. G. Hoo e , and C. G. Ho-
o e , Phys. Re . A41, 4552 (1990).