Lattice-switch Monte Carlo method: application to soft potentials.
Physical Review E - Statistical, Nonlinear and Soft Matter Physics (2002)
- PubMed: 11909310
Available from
Andrew Jackson's profile on Mendeley.
or
Abstract
The lattice-switch Monte Carlo method, recently introduced and applied in the context of hard spheres, is extended to particles interacting through a soft potential. The method utilizes a transformation that switches between configurations of two different crystalline structures, allowing the phase space of both structures to be explored in a single simulation and the difference between their free energies to be determined directly. We apply the method to determine the fcc-hcp crystalline phase behavior of the classical Lennard-Jones solid.
Available from
Andrew Jackson's profile on Mendeley.
Page 1
Lattice-switch Monte Carlo method: application to soft potentials.
oru
Ed
;
tro
tia
ct
w
ph
S
PHYSICAL REVIEW E, VOLUME 65, 036710~i.e., crystal-to-crystal transitions induced by changes of tem-
perature and pressure! @9#; and harmonic crystal dynamics
and its refinements @10# provide an alternative, approximate,
strategy against which to benchmark the LSMC calculations.
The LJ model is also of rather wider physical interest: it
provides an excellent account of many features of the behav-
ior of rare gas solids ~RGS! @11#; one might hope it would
also account for the observed phase behavior of RGS, at
least under ‘‘normal’’ conditions @12#. But the observed
phase behavior ~the favored RGS crystal structure appears to
be fcc, almost everywhere @13#! has proved surprisingly hard
to elucidate @14#. The difference between the ground-state
energies of hcp and fcc actually favors the hcp structure; the
difference is small and its sign can be changed both by small
intrinsic smallness. We thus have confidence that the implied
phase diagram is truly representative of the classical LJ
model. Finally in Sec. VIII, we review what we have learned
from this study about the LS methodology and its value in
relation to the many other strategies that have been used to
attack the problem of free-energy measurement; and about
LJ systems themselves.
II. MODEL AND GROUND STATE
We consider a system of N particles of spatial coordinates
$
rW
%
, confined within a volume V and subject to periodic
boundary conditions. Each particle is associated with aLattice-switch Monte Carlo meth
A. N. Jackson, A. D. B
Department of Physics and Astronomy, The University of
~Received 6 October 2001
The lattice-switch Monte Carlo method, recently in
extended to particles interacting through a soft poten
between configurations of two different crystalline stru
explored in a single simulation and the difference bet
apply the method to determine the fcc-hcp crystalline
DOI: 10.1103/PhysRevE.65.036710 PAC
I. INTRODUCTION
In recent papers @1,2#, we introduced a new technique for
tackling an old problem. The technique is called lattice-
switch Monte Carlo ~LSMC!; the problem is the determina-
tion of the relative stability of two crystalline structures. The
LSMC method is built on a transformation that maps a con-
figuration of one structure onto a candidate configuration of
the other by ‘‘switching’’ one set of lattice vectors for the
other, while keeping the displacements with respect to the
lattice sites constant. The sampling of the displacement con-
figurations is multicanonically biased @3# to favor paths lead-
ing to gateway arrangements for which the Monte Carlo
switch to the candidate configuration will be accepted. The
configurations of both structures are sampled within a single
simulation, and the difference between their free energies
evaluated from the ratio of their measured probabilities.
In @1,2# we set out the LSMC method in some detail; we
argued that it offers significant advances with respect to ex-
isting strategies @4#, preeminently integration methods @5,6#;
and we explored its operation through extensive studies of
the relative stability of the crystalline ~fcc, hcp! phases of
hard spheres. More recently @7#, we have shown that LSMC
can be set in a somewhat more general framework that al-
lows exploration of freezing ~crystal-liquid phase coexist-
ence!; that work also was set in the context of hard spheres.
In this paper we describe the application of the method to
model crystalline solids with a soft interaction potential of
Lennard-Jones ~LJ! form. The generalization from ‘‘hard’’ to
‘‘soft’’ potentials has a number of consequences. Tempera-
ture and pressure become independent control parameters
@8#. There is the possibility of real structural phase behavior1063-651X/2002/65~3!/036710~12!/$20.00 65 0367d: Application to soft potentials
ce, and G. J. Ackland
inburgh, Edinburgh EH9 3JZ, Scotland, United Kingdom
published 7 March 2002!
duced and applied in the context of hard spheres, is
l. The method utilizes a transformation that switches
ures, allowing the phase space of both structures to be
een their free energies to be determined directly. We
ase behavior of the classical Lennard-Jones solid.
number~s!: 02.70.Rr, 05.10.Ln, 65.20.1w, 64.70.Kb
deformations of the potential @15# and by truncating its range
@16–18#. In such circumstances fluctuation effects are clearly
crucial. Such effects may have either an anharmonic or a
quantum character ~sometimes both! and will necessarily do
so, at high enough or low enough temperatures, respectively.
A satisfactory, general, computational strategy must deal co-
herently with both anharmonic and quantum effects. The
LSMC method can be developed to do so. But in this work
~in keeping with the majority of existing simulation ap-
proaches to this problem! we shall neglect quantum effects:
we study the structural phase behavior of the ~one might say
a, given the sensitivity to the potential! classical LJ model.
The layout of the paper is as follows. Section II defines
the model, and explores the sensitivity of the ground state to
the range at which the potential is truncated. In Sec. III we
summarize the statistical mechanics relevant to the phase sta-
bility problem. Section IV describes the standard approach to
this problem, based on the small-amplitude ~quasiharmonic!
treatment of the fluctuation spectrum, and again explores the
consequences of potential truncation. Section V formulates
the LSMC method in the context of soft potentials; Sec. VII
presents the results it delivers. They show that the LSMC
method captures the results of the harmonic theory in the
regime ~high enough density, low enough temperature!
where one would expect the harmonic approximation to be
satisfactory; and that it extends naturally beyond this regime,
revealing anharmonic corrections that grow with reducing
density, in an intelligible way. Since the method focuses di-
rectly on the difference of the free energies between the
phases it is able to measure it with precision, despite its©2002 The American Physical Society10-1
Ed
;
tro
tia
ct
w
ph
S
PHYSICAL REVIEW E, VOLUME 65, 036710~i.e., crystal-to-crystal transitions induced by changes of tem-
perature and pressure! @9#; and harmonic crystal dynamics
and its refinements @10# provide an alternative, approximate,
strategy against which to benchmark the LSMC calculations.
The LJ model is also of rather wider physical interest: it
provides an excellent account of many features of the behav-
ior of rare gas solids ~RGS! @11#; one might hope it would
also account for the observed phase behavior of RGS, at
least under ‘‘normal’’ conditions @12#. But the observed
phase behavior ~the favored RGS crystal structure appears to
be fcc, almost everywhere @13#! has proved surprisingly hard
to elucidate @14#. The difference between the ground-state
energies of hcp and fcc actually favors the hcp structure; the
difference is small and its sign can be changed both by small
intrinsic smallness. We thus have confidence that the implied
phase diagram is truly representative of the classical LJ
model. Finally in Sec. VIII, we review what we have learned
from this study about the LS methodology and its value in
relation to the many other strategies that have been used to
attack the problem of free-energy measurement; and about
LJ systems themselves.
II. MODEL AND GROUND STATE
We consider a system of N particles of spatial coordinates
$
rW
%
, confined within a volume V and subject to periodic
boundary conditions. Each particle is associated with aLattice-switch Monte Carlo meth
A. N. Jackson, A. D. B
Department of Physics and Astronomy, The University of
~Received 6 October 2001
The lattice-switch Monte Carlo method, recently in
extended to particles interacting through a soft poten
between configurations of two different crystalline stru
explored in a single simulation and the difference bet
apply the method to determine the fcc-hcp crystalline
DOI: 10.1103/PhysRevE.65.036710 PAC
I. INTRODUCTION
In recent papers @1,2#, we introduced a new technique for
tackling an old problem. The technique is called lattice-
switch Monte Carlo ~LSMC!; the problem is the determina-
tion of the relative stability of two crystalline structures. The
LSMC method is built on a transformation that maps a con-
figuration of one structure onto a candidate configuration of
the other by ‘‘switching’’ one set of lattice vectors for the
other, while keeping the displacements with respect to the
lattice sites constant. The sampling of the displacement con-
figurations is multicanonically biased @3# to favor paths lead-
ing to gateway arrangements for which the Monte Carlo
switch to the candidate configuration will be accepted. The
configurations of both structures are sampled within a single
simulation, and the difference between their free energies
evaluated from the ratio of their measured probabilities.
In @1,2# we set out the LSMC method in some detail; we
argued that it offers significant advances with respect to ex-
isting strategies @4#, preeminently integration methods @5,6#;
and we explored its operation through extensive studies of
the relative stability of the crystalline ~fcc, hcp! phases of
hard spheres. More recently @7#, we have shown that LSMC
can be set in a somewhat more general framework that al-
lows exploration of freezing ~crystal-liquid phase coexist-
ence!; that work also was set in the context of hard spheres.
In this paper we describe the application of the method to
model crystalline solids with a soft interaction potential of
Lennard-Jones ~LJ! form. The generalization from ‘‘hard’’ to
‘‘soft’’ potentials has a number of consequences. Tempera-
ture and pressure become independent control parameters
@8#. There is the possibility of real structural phase behavior1063-651X/2002/65~3!/036710~12!/$20.00 65 0367d: Application to soft potentials
ce, and G. J. Ackland
inburgh, Edinburgh EH9 3JZ, Scotland, United Kingdom
published 7 March 2002!
duced and applied in the context of hard spheres, is
l. The method utilizes a transformation that switches
ures, allowing the phase space of both structures to be
een their free energies to be determined directly. We
ase behavior of the classical Lennard-Jones solid.
number~s!: 02.70.Rr, 05.10.Ln, 65.20.1w, 64.70.Kb
deformations of the potential @15# and by truncating its range
@16–18#. In such circumstances fluctuation effects are clearly
crucial. Such effects may have either an anharmonic or a
quantum character ~sometimes both! and will necessarily do
so, at high enough or low enough temperatures, respectively.
A satisfactory, general, computational strategy must deal co-
herently with both anharmonic and quantum effects. The
LSMC method can be developed to do so. But in this work
~in keeping with the majority of existing simulation ap-
proaches to this problem! we shall neglect quantum effects:
we study the structural phase behavior of the ~one might say
a, given the sensitivity to the potential! classical LJ model.
The layout of the paper is as follows. Section II defines
the model, and explores the sensitivity of the ground state to
the range at which the potential is truncated. In Sec. III we
summarize the statistical mechanics relevant to the phase sta-
bility problem. Section IV describes the standard approach to
this problem, based on the small-amplitude ~quasiharmonic!
treatment of the fluctuation spectrum, and again explores the
consequences of potential truncation. Section V formulates
the LSMC method in the context of soft potentials; Sec. VII
presents the results it delivers. They show that the LSMC
method captures the results of the harmonic theory in the
regime ~high enough density, low enough temperature!
where one would expect the harmonic approximation to be
satisfactory; and that it extends naturally beyond this regime,
revealing anharmonic corrections that grow with reducing
density, in an intelligible way. Since the method focuses di-
rectly on the difference of the free energies between the
phases it is able to measure it with precision, despite its©2002 The American Physical Society10-1
Page 2
A. N. JACKSON, A. D. BRUCE, AND G. J. ACKLAND PHYSICAL REVIEW E 65 036710site of a lattice @19# defined by a set of vectors
$
RW
%
a
where a
indexes the lattice type @20#.
The particles interact through some approximation to the
LJ potential
f
~
ri j!54eF S
s
ri j
D
12
2
S
s
ri j
D
6
G
, ~1!
where ri j5urW i2rW ju. The customary approximation is some
form of truncation, which is required, in principle, to avoid
self-interactions in a finite periodic system and is also a prac-
tical expedient to limit computational cost. One may define a
truncated potential in a number of different ways. Thus, one
may truncate the potential at a fixed cutoff rx5ms ~for some
chosen m), defining
fa~ri j!5H
f
~
ri j!2f~rx! ri j<rx
0 ri j.rx
. ~2!
In this case the potential is also ‘‘shifted’’ so as to suppress
the discontinuous change occurring when particle separa-
tions pass through the cutoff value, which they do, both sto-
chastically in the course of the crystal dynamics and system-
atically with changes in control parameters such as density or
pressure.
Alternatively one may choose a truncation scheme in
which pairs of particles are identified at the outset as inter-
acting or not, according to the separation of the lattice sites
with which they are associated. Thus, setting rx5mr0 where
r0 is the equilibrium nearest-neighbor separation, we define
fb~ri j!5H
f
~
ri j! Ri j<rx
0 Ri j.rx
, ~3!
where Ri j[uRW i2RW ju. The fact that the cutoff rx now scales
appropriately with density ensures that the list of interacting
particles prescribed by Eq. ~3! is indeed fixed ~independent
of density! so that a potential shift is redundant in this case.
The significance of the choice of truncation is revealed in
the ~classical! ground-state energy
E05E
~
$
RW
%
!5
(
^
i j
&
f
~
Ri j!, ~4!
where the sum counts ~once! the contribution of each pair of
interacting particles. Figure 1 shows the difference @21# be-
tween the ground-state energies of the two structures ~per
particle!
De
a ,a
˜
0
[
E
~
$
RW
%
a
!2E~
$
RW
%
a
˜
!
N , ~5!
as a function of density, with a[fcc and a˜ [hcp. In addition
to the exact result following from implementing Eq. ~5!
without approximation @22–24#, we show the results that fol-
low if the potential is approximated by the truncated forms
fa @Eq. ~2!# and fb @Eq. ~3!#.03671The exact result shows the well-known T50 phase be-
havior @11#: there is a single-phase boundary ~at r5r0
52.1727s23) separating hcp (r,ro) and fcc (r.ro) re-
gions. By contrast the fixed-cutoff approximation @Eq. ~2!
with rx52.5s# leads to a ground-state energy that varies
wildly with density, as the number of particles within the
interaction range evolves. The implied ~seemingly rich!
phase behavior @17# can be traced directly to the truncation.
The fixed-interaction-list approximation @Eq. ~3! with rx
52.5r0# fares differently, but no better: the evolution is
smooth, but shows no phase boundary at all.
The moral is clear: to capture phase behavior that is any
trustworthy sense representative of the true LJ model one
must handle the ground-state energies without recourse to
truncation. While it is not our aim to give a definitive treat-
ment of the LJ problem ~our treatment is classical, after all!
we shall respect this constraint. Thus, in what follows, we
shall take for the general configurational energy the form
E
~
$
rW
%
!5E01
(
^
i j
&
fc~ri j!, ~6!
with
fc~ri j!5H
f
~
ri j!2f~Ri j! Ri j<rx
0 Ri j.rx
. ~7!
The ground-state energy E0 is treated essentially exactly
@25#. The effects of the truncation are restricted to the fluc-
tuation spectrum; we shall discuss them ~and the choice of
rx) in due course.
FIG. 1. The difference between the ground-state energies per
particle for fcc and hcp structures, as function of density. The solid
line shows the result with no potential truncation. The dotted and
dashed lines show the results based on the fixed-cutoff truncation
scheme @Eq. ~2! with rx52.5s# and the fixed-neighbor-list trunca-
tion scheme @Eq. ~3! with rx52.5r0# respectively. Note the order-
of-magnitude difference between the scales on the left and right
hand axes that are used for ‘‘exact’’ and ‘‘truncated’’ results respec-
tively.0-2
$
RW
%
a
where a
indexes the lattice type @20#.
The particles interact through some approximation to the
LJ potential
f
~
ri j!54eF S
s
ri j
D
12
2
S
s
ri j
D
6
G
, ~1!
where ri j5urW i2rW ju. The customary approximation is some
form of truncation, which is required, in principle, to avoid
self-interactions in a finite periodic system and is also a prac-
tical expedient to limit computational cost. One may define a
truncated potential in a number of different ways. Thus, one
may truncate the potential at a fixed cutoff rx5ms ~for some
chosen m), defining
fa~ri j!5H
f
~
ri j!2f~rx! ri j<rx
0 ri j.rx
. ~2!
In this case the potential is also ‘‘shifted’’ so as to suppress
the discontinuous change occurring when particle separa-
tions pass through the cutoff value, which they do, both sto-
chastically in the course of the crystal dynamics and system-
atically with changes in control parameters such as density or
pressure.
Alternatively one may choose a truncation scheme in
which pairs of particles are identified at the outset as inter-
acting or not, according to the separation of the lattice sites
with which they are associated. Thus, setting rx5mr0 where
r0 is the equilibrium nearest-neighbor separation, we define
fb~ri j!5H
f
~
ri j! Ri j<rx
0 Ri j.rx
, ~3!
where Ri j[uRW i2RW ju. The fact that the cutoff rx now scales
appropriately with density ensures that the list of interacting
particles prescribed by Eq. ~3! is indeed fixed ~independent
of density! so that a potential shift is redundant in this case.
The significance of the choice of truncation is revealed in
the ~classical! ground-state energy
E05E
~
$
RW
%
!5
(
^
i j
&
f
~
Ri j!, ~4!
where the sum counts ~once! the contribution of each pair of
interacting particles. Figure 1 shows the difference @21# be-
tween the ground-state energies of the two structures ~per
particle!
De
a ,a
˜
0
[
E
~
$
RW
%
a
!2E~
$
RW
%
a
˜
!
N , ~5!
as a function of density, with a[fcc and a˜ [hcp. In addition
to the exact result following from implementing Eq. ~5!
without approximation @22–24#, we show the results that fol-
low if the potential is approximated by the truncated forms
fa @Eq. ~2!# and fb @Eq. ~3!#.03671The exact result shows the well-known T50 phase be-
havior @11#: there is a single-phase boundary ~at r5r0
52.1727s23) separating hcp (r,ro) and fcc (r.ro) re-
gions. By contrast the fixed-cutoff approximation @Eq. ~2!
with rx52.5s# leads to a ground-state energy that varies
wildly with density, as the number of particles within the
interaction range evolves. The implied ~seemingly rich!
phase behavior @17# can be traced directly to the truncation.
The fixed-interaction-list approximation @Eq. ~3! with rx
52.5r0# fares differently, but no better: the evolution is
smooth, but shows no phase boundary at all.
The moral is clear: to capture phase behavior that is any
trustworthy sense representative of the true LJ model one
must handle the ground-state energies without recourse to
truncation. While it is not our aim to give a definitive treat-
ment of the LJ problem ~our treatment is classical, after all!
we shall respect this constraint. Thus, in what follows, we
shall take for the general configurational energy the form
E
~
$
rW
%
!5E01
(
^
i j
&
fc~ri j!, ~6!
with
fc~ri j!5H
f
~
ri j!2f~Ri j! Ri j<rx
0 Ri j.rx
. ~7!
The ground-state energy E0 is treated essentially exactly
@25#. The effects of the truncation are restricted to the fluc-
tuation spectrum; we shall discuss them ~and the choice of
rx) in due course.
FIG. 1. The difference between the ground-state energies per
particle for fcc and hcp structures, as function of density. The solid
line shows the result with no potential truncation. The dotted and
dashed lines show the results based on the fixed-cutoff truncation
scheme @Eq. ~2! with rx52.5s# and the fixed-neighbor-list trunca-
tion scheme @Eq. ~3! with rx52.5r0# respectively. Note the order-
of-magnitude difference between the scales on the left and right
hand axes that are used for ‘‘exact’’ and ‘‘truncated’’ results respec-
tively.0-2
Page 3
LATTICE-SWITCH MONTE CARLO METHOD: . . . PHYSICAL REVIEW E 65 036710III. STATISTICAL MECHANICS
The fluctuation spectrum ~the ‘‘significant’’ configurations
$
rW
%
) of a crystalline phase a is naturally described by the
displacements with respect to the lattice sites
$
RW
%
a
. Thus we
make the familiar decomposition
rW i5RW i
a
1uW i . ~8!
The partition function ~configurational weight! of structure
a , at temperature T, may then be written as
Z
~
N ,T ,V ,a!5
)
i
F
E
[a]
duW i
G
e2bE($r
W
%
)
. ~9!
Here * [a] signifies integration subject to an appropriate con-
figurational constraint that picks out ~from the entire con-
figuration space of N particles in volume V) those configu-
rations that ‘‘belong’’ to crystalline phase a @26#. The
partition function defines the Helmholtz free-energy density
through
f
~
N ,T ,V ,a!52
1
Nb ln Z~N ,T ,V ,a!. ~10!
The relative stability of the two phases reflects the difference
between their free energies, and thus the ratio of their parti-
tion functions @27#
D f
aa
˜
[ f
~
N ,T ,V ,a!2 f ~N ,T ,V ,a˜ !52 1Nb ln Raa˜ ~N ,T ,V !,
~11!
where
R
aa
˜
~
N ,T ,V ![
Z
~
N ,T ,V ,a!
Z
~
N ,T ,V ,a˜ !
. ~12!
If the difference between the densities of the two phases is
small, the phase boundary follows immediately as the locus
of points for which D f
aa
˜ vanishes ~the two phases have
equal statistical weights!. In general, however, the problem is
more conveniently formulated in the NTP ensemble. The
relevant partition functions are then
Z
~
N ,T ,P ,a!5
E
dVZ
~
N ,T ,V ,a!e2bPV. ~13!
The difference between the Gibbs free-energy densities fol-
lows as
Dg
aa
˜
[g
~
N ,T ,P ,a!2g~N ,T ,P ,a˜ !
52
1
Nb ln Raa˜ ~N ,T ,P !, ~14!
where now
R
aa
˜
~
N ,T ,P ![
Z
~
N ,T ,P ,a!
Z
~
N ,T ,P ,a˜ !
. ~15!03671The phase boundary is defined by the condition that
R
aa
˜
~
N ,T ,P !51, ~16!
in the thermodynamic (N→‘) limit.
IV. HARMONIC APPROXIMATION
The foundation stone of analytical approaches to the prob-
lem of free-energy determination for crystalline solids is the
harmonic approximation. The methodology is long estab-
lished @10#; many studies of LJ systems within the harmonic
approximation ~and its refinements! exist @28#. Here we use
the harmonic framework to provide a benchmark of the
LSMC work to be presented in following sections. We de-
scribe it in outline only.
In the harmonic approximation we expand the configura-
tional energy @Eq. ~6!# to second order in the displacements
$
uW
%
@29#
E
~
$
rW
%
!5E~
$
RW
%
!1
1
2 (
^
i j
&
,m ,n
Ki j
mnui
mu j
n
1 . ~17!
Here ui
m and u j
n denote cartesian components (m , n) of the
displacements uW of the particles associated with sites i , j .
The ‘‘dynamical’’ matrix K is defined by @30#
Ki j
mn
5
H
F
fc8~Ri j!
Ri j
2fc9~Ri j!Gnˆ i j
mnˆ i j
n
2
fc8~Ri j!
Ri j
d
mn
~
jÞi !
2
(
kÞi
Kik
mn
~
j5i !
~18!
with
nˆ i j[
RW i j
Ri j
. ~19!
As a result of the truncation scheme adopted in Eq. ~6! the
ground-state energy in Eq. ~17! is essentially exact ~in the
limit of large enough N), while the fluctuation term contains
contributions only from interactions falling within the
~density-dependent! cutoff.
Within this approximation the NVT partition function
@Eq. ~9!# can be expressed in the form
Z
~
N ,T ,V ,a!5@2p/b#3N/2exp@2bE~
$
RW
%
a
!#@
det K
~
a!#
21/2
,
~20!
where K(a) is the 3N33N matrix whose elements are de-
fined by Eq. ~18! ~applied to phase a). Then the fluctuation
contribution to the free-energy density of phase a in the
harmonic approximation is measured @31# by
f
a
h
[
1
2Nb ln@det K~a!#5
1
2Nb (j51
3(N21)
ln l j
a
, ~21!0-3
The fluctuation spectrum ~the ‘‘significant’’ configurations
$
rW
%
) of a crystalline phase a is naturally described by the
displacements with respect to the lattice sites
$
RW
%
a
. Thus we
make the familiar decomposition
rW i5RW i
a
1uW i . ~8!
The partition function ~configurational weight! of structure
a , at temperature T, may then be written as
Z
~
N ,T ,V ,a!5
)
i
F
E
[a]
duW i
G
e2bE($r
W
%
)
. ~9!
Here * [a] signifies integration subject to an appropriate con-
figurational constraint that picks out ~from the entire con-
figuration space of N particles in volume V) those configu-
rations that ‘‘belong’’ to crystalline phase a @26#. The
partition function defines the Helmholtz free-energy density
through
f
~
N ,T ,V ,a!52
1
Nb ln Z~N ,T ,V ,a!. ~10!
The relative stability of the two phases reflects the difference
between their free energies, and thus the ratio of their parti-
tion functions @27#
D f
aa
˜
[ f
~
N ,T ,V ,a!2 f ~N ,T ,V ,a˜ !52 1Nb ln Raa˜ ~N ,T ,V !,
~11!
where
R
aa
˜
~
N ,T ,V ![
Z
~
N ,T ,V ,a!
Z
~
N ,T ,V ,a˜ !
. ~12!
If the difference between the densities of the two phases is
small, the phase boundary follows immediately as the locus
of points for which D f
aa
˜ vanishes ~the two phases have
equal statistical weights!. In general, however, the problem is
more conveniently formulated in the NTP ensemble. The
relevant partition functions are then
Z
~
N ,T ,P ,a!5
E
dVZ
~
N ,T ,V ,a!e2bPV. ~13!
The difference between the Gibbs free-energy densities fol-
lows as
Dg
aa
˜
[g
~
N ,T ,P ,a!2g~N ,T ,P ,a˜ !
52
1
Nb ln Raa˜ ~N ,T ,P !, ~14!
where now
R
aa
˜
~
N ,T ,P ![
Z
~
N ,T ,P ,a!
Z
~
N ,T ,P ,a˜ !
. ~15!03671The phase boundary is defined by the condition that
R
aa
˜
~
N ,T ,P !51, ~16!
in the thermodynamic (N→‘) limit.
IV. HARMONIC APPROXIMATION
The foundation stone of analytical approaches to the prob-
lem of free-energy determination for crystalline solids is the
harmonic approximation. The methodology is long estab-
lished @10#; many studies of LJ systems within the harmonic
approximation ~and its refinements! exist @28#. Here we use
the harmonic framework to provide a benchmark of the
LSMC work to be presented in following sections. We de-
scribe it in outline only.
In the harmonic approximation we expand the configura-
tional energy @Eq. ~6!# to second order in the displacements
$
uW
%
@29#
E
~
$
rW
%
!5E~
$
RW
%
!1
1
2 (
^
i j
&
,m ,n
Ki j
mnui
mu j
n
1 . ~17!
Here ui
m and u j
n denote cartesian components (m , n) of the
displacements uW of the particles associated with sites i , j .
The ‘‘dynamical’’ matrix K is defined by @30#
Ki j
mn
5
H
F
fc8~Ri j!
Ri j
2fc9~Ri j!Gnˆ i j
mnˆ i j
n
2
fc8~Ri j!
Ri j
d
mn
~
jÞi !
2
(
kÞi
Kik
mn
~
j5i !
~18!
with
nˆ i j[
RW i j
Ri j
. ~19!
As a result of the truncation scheme adopted in Eq. ~6! the
ground-state energy in Eq. ~17! is essentially exact ~in the
limit of large enough N), while the fluctuation term contains
contributions only from interactions falling within the
~density-dependent! cutoff.
Within this approximation the NVT partition function
@Eq. ~9!# can be expressed in the form
Z
~
N ,T ,V ,a!5@2p/b#3N/2exp@2bE~
$
RW
%
a
!#@
det K
~
a!#
21/2
,
~20!
where K(a) is the 3N33N matrix whose elements are de-
fined by Eq. ~18! ~applied to phase a). Then the fluctuation
contribution to the free-energy density of phase a in the
harmonic approximation is measured @31# by
f
a
h
[
1
2Nb ln@det K~a!#5
1
2Nb (j51
3(N21)
ln l j
a
, ~21!0-3
Page 4
A. N. JACKSON, A. D. BRUCE, AND G. J. ACKLAND PHYSICAL REVIEW E 65 036710where l j
a
@
j51, . . . ,3(N21)
#
denote the nonzero eigen-
values of the matrix K(a). The difference between the har-
monic contributions to the free-energy densities of the two
phases follows as
D f
aa
˜
h
5 f
a
h
2 f
a
˜
h
5
1
2Nb lnF
det K
~
a!
det K
~
a
˜
!
G
. ~22!
These expressions can be evaluated numerically with rela-
tive ease ~for modest N values!, allowing investigation of the
influence of the potential truncation. Figure 2 shows f h, as
defined by Eq. ~21! @32#, for the two structures in systems of
N5216 particles, at a density rs35A2, for different choices
of cutoff. We have chosen @33# to parameterize the cutoff by
the inverse of the number nx of interacting neighbors ~the
number of particles falling within the cutoff range!. The free-
energy ‘‘scale’’ on which to assess the significance of these
results is set by the difference between the values for the two
structures. On this scale the values of f h for the two struc-
tures each evolve significantly between the two extremes.
But the difference between these free energies varies rela-
tively little: the values for the two extremes differ by only
2%. Figure 3 shows that this statement remains true over a
wide range of densities, and for a substantially larger system
size N. The conclusion we draw is that, in contrast to the
ground-state energy ~Sec. II, Fig. 1!, the fluctuation contri-
bution to the free-energy difference of interest in the har-
monic approximation can be accurately captured by a model
incorporating only relatively short-range interactions. It
might seem reasonable to suppose that this state of affairs
extends beyond the harmonic limit and we have indeed pro-
ceeded on this basis: for the most part our LSMC studies
have focused on the truncated model defined in Eq. ~7!, with
rx51.5r0 and nx518. But the results presented below in-
clude a retrospective check on the consistency of this as-
sumption.
Figure 3 also shows that the dependence of the free-
energy difference upon system size N ~finite-size effects! is
FIG. 2. The harmonic contributions to the free-energy density
@Eq. ~21!; note @32## for fcc and hcp structures with N5216, and
rs
3
5
A2. The results are expressed as functions of the number of
interacting neighbors nx @33#.03671significant on the relevant scale: the harmonic free-energy
difference changes by some 25% between N5635216 and
N512351728. Our direct matrix-diagonalization treatment
is not readily extended to larger N values. But our results for
the harmonic free-energy difference between the two struc-
tures for N5123 lie within 3% of the harmonic thermody-
namic limit reported in @34#; and we find that a 1/N extrapo-
lation of our N563 and N5123 results is consistent with
that limit. Accordingly our LSMC studies have also been
largely restricted to these two system sizes.
V. LATTICE-SWITCH MC METHOD
The problem of phase behavior is succinctly expressed in
Eqs. ~12! or ~15!: the relative stability of ~free-energy differ-
ence between! two phases is determined by the ratio of the
associated partition functions. It is helpful ~both conceptually
and practically! to recognize that ratio for what it is: the ratio
of the equilibrium probabilities with which a system, free to
explore both phases, will visit each one. The LSMC method
provides the strategy needed to captitalize on this identifica-
tion. We have described the method fully elsewhere @2#, in
the context of a system of hard spheres. The description here
will, therefore, be brief and will focus on the adaptations
needed to accommodate soft interactions.
In principle, LSMC represents a simple and intuitive ex-
tension of the MC framework @6# conventionally used for
exploration of a single ~structural! phase. The ‘‘single-phase’’
MC method updates the particle displacements for fixed lat-
tice vectors; LSMC augments this procedure by updates of
the lattice vectors with the displacements held fixed. By ‘‘up-
date’’ of the lattice vectors we mean replacing the entire set
of lattice vectors describing the current structure with a set
that describes the other structure. The current particle con-
figuration, which ‘‘belongs’’ to the phase space associated
with the current structure, is thus mapped onto a particle
configuration that ‘‘belongs’’ to the phase space of the other
structure.
Without further refinement this simple procedure will not
FIG. 3. The difference between the harmonic contributions to
the free-energy density @Eq. ~21!# of the two structures for system
sizes N5216 and N51728 as function of density. The solid and
dashed lines correspond to different fixed-neighbor-list truncation
schemes.0-4
a
@
j51, . . . ,3(N21)
#
denote the nonzero eigen-
values of the matrix K(a). The difference between the har-
monic contributions to the free-energy densities of the two
phases follows as
D f
aa
˜
h
5 f
a
h
2 f
a
˜
h
5
1
2Nb lnF
det K
~
a!
det K
~
a
˜
!
G
. ~22!
These expressions can be evaluated numerically with rela-
tive ease ~for modest N values!, allowing investigation of the
influence of the potential truncation. Figure 2 shows f h, as
defined by Eq. ~21! @32#, for the two structures in systems of
N5216 particles, at a density rs35A2, for different choices
of cutoff. We have chosen @33# to parameterize the cutoff by
the inverse of the number nx of interacting neighbors ~the
number of particles falling within the cutoff range!. The free-
energy ‘‘scale’’ on which to assess the significance of these
results is set by the difference between the values for the two
structures. On this scale the values of f h for the two struc-
tures each evolve significantly between the two extremes.
But the difference between these free energies varies rela-
tively little: the values for the two extremes differ by only
2%. Figure 3 shows that this statement remains true over a
wide range of densities, and for a substantially larger system
size N. The conclusion we draw is that, in contrast to the
ground-state energy ~Sec. II, Fig. 1!, the fluctuation contri-
bution to the free-energy difference of interest in the har-
monic approximation can be accurately captured by a model
incorporating only relatively short-range interactions. It
might seem reasonable to suppose that this state of affairs
extends beyond the harmonic limit and we have indeed pro-
ceeded on this basis: for the most part our LSMC studies
have focused on the truncated model defined in Eq. ~7!, with
rx51.5r0 and nx518. But the results presented below in-
clude a retrospective check on the consistency of this as-
sumption.
Figure 3 also shows that the dependence of the free-
energy difference upon system size N ~finite-size effects! is
FIG. 2. The harmonic contributions to the free-energy density
@Eq. ~21!; note @32## for fcc and hcp structures with N5216, and
rs
3
5
A2. The results are expressed as functions of the number of
interacting neighbors nx @33#.03671significant on the relevant scale: the harmonic free-energy
difference changes by some 25% between N5635216 and
N512351728. Our direct matrix-diagonalization treatment
is not readily extended to larger N values. But our results for
the harmonic free-energy difference between the two struc-
tures for N5123 lie within 3% of the harmonic thermody-
namic limit reported in @34#; and we find that a 1/N extrapo-
lation of our N563 and N5123 results is consistent with
that limit. Accordingly our LSMC studies have also been
largely restricted to these two system sizes.
V. LATTICE-SWITCH MC METHOD
The problem of phase behavior is succinctly expressed in
Eqs. ~12! or ~15!: the relative stability of ~free-energy differ-
ence between! two phases is determined by the ratio of the
associated partition functions. It is helpful ~both conceptually
and practically! to recognize that ratio for what it is: the ratio
of the equilibrium probabilities with which a system, free to
explore both phases, will visit each one. The LSMC method
provides the strategy needed to captitalize on this identifica-
tion. We have described the method fully elsewhere @2#, in
the context of a system of hard spheres. The description here
will, therefore, be brief and will focus on the adaptations
needed to accommodate soft interactions.
In principle, LSMC represents a simple and intuitive ex-
tension of the MC framework @6# conventionally used for
exploration of a single ~structural! phase. The ‘‘single-phase’’
MC method updates the particle displacements for fixed lat-
tice vectors; LSMC augments this procedure by updates of
the lattice vectors with the displacements held fixed. By ‘‘up-
date’’ of the lattice vectors we mean replacing the entire set
of lattice vectors describing the current structure with a set
that describes the other structure. The current particle con-
figuration, which ‘‘belongs’’ to the phase space associated
with the current structure, is thus mapped onto a particle
configuration that ‘‘belongs’’ to the phase space of the other
structure.
Without further refinement this simple procedure will not
FIG. 3. The difference between the harmonic contributions to
the free-energy density @Eq. ~21!# of the two structures for system
sizes N5216 and N51728 as function of density. The solid and
dashed lines correspond to different fixed-neighbor-list truncation
schemes.0-4
Page 5
LATTICE-SWITCH MONTE CARLO METHOD: . . . PHYSICAL REVIEW E 65 036710in general work: an equilibrium configuration of one struc-
ture will typically have as its conjugate ~i.e., will map into! a
relatively high-energy configuration of the other structure.
The lattice switch MC step will typically be rejected. To
make it work the sampling of the configurations of each
structure must be multicanonically biased @3# to enhance the
probabilities of those states ~we call them gateway states!
from which a switch can be attempted with a reasonable
likelihood of success. The bias needs to reflect the difference
between the excitation energies @35# of conjugate pairs. To
that end we define an order parameter @36#
M[M
~
$
uW
%
!5E~
$
uW
%
,
$
RW
%hcp!2E~$uW %,$RW %fcc!, ~23!
where
E
~
$
uW
%
,
$
RW
%
a
!5b@
E
~
$
rW
%
!2E~
$
RW
%
a
!#
, ~24!
measures the excitation energy of configuration
$
rW
%
[
$
uW
%
,a
~the amount by which it exceeds the ground-state energy of
structure a) in units of kT5b21. The energy difference de-
fined in Eq. ~23! is the analog of the ‘‘overlap order param-
eter’’ utilized in our studies of hard spheres @1,2#.
Given the relationship between the energy of a configura-
tion and the energy of its conjugate one sees that displace-
ment patterns
$
uW
%
drawn from equilibrium sampling of the
fcc ~hcp! configuration space will give predominantly posi-
tive ~negative! values of
$
M
%
. The gateway states lie in be-
tween @37#; in the thermodynamic limit they have vanish-
ingly small equilibrium probability. We thus require to
engineer a multicanonical sampling algorithm to enhance the
probability along a notional line in M space, extending from
the ‘‘equilibrium’’ M values ~whose magnitude reflects the
energy cost of a LS acting on a typical configuration!
through to the small-M gateway configurations. This aim is
realized by augmenting the system energy function Eq. ~6!
appropriately,
bE
~
$
rW
%
!→bE~
$
rW
%
!1h@
M
#
[E˜
~
$
rW
%
!. ~25!
Here h
@
M
#
is a multicanonical weight function @3# that
must be chosen so as to allow the system to access the gate-
way configurations. For this purpose it is sufficient ~though
not necessarily optimal! to construct the weights such that
the multicanonical distribution P(MuN ,T ,
$
h
%
. . . ) @38# is
roughly flat. The desired ratio of partition functions @Eq. ~12!
or ~15!#, is identified with the ratio of the aggregate canoni-
cal probabilities of the two phases
R
aa
˜
~
N ,T , . . . ![
P
~
auN ,T , . . . !
P
~
a
˜
uN ,T , . . . !
5
(
$
M
%
P
~
M,auN ,T , . . . !
(
$
M
%
P
~
M,a˜ uN ,T , . . . !
~26!
and can be obtained from the corresponding ~measured! mul-
ticanonical distributions by reweighting ~folding out the
weights!03671(
$
M
%
P
~
M,auN ,T , . . . !
(
$
M
%
P
~
M,a˜ uN ,T , . . . !
5
(
$
M
%
eh[M]P
~
M,auN ,T ,
$
h
%
, . . . !
(
$
M
%
eh[M]P
~
M,a˜ uN ,T ,
$
h
%
, . . . !
. ~27!
The last step removes the explicit bias resulting from the
weights, leaving only the desired legacy @39#—that both
phases are visited with the canonical probabilities. Notice
that, in contrast to the hard-sphere problem, where M serves
as a unique identifier of the phase, the relevant sums in Eqs.
~26! and ~27! involve the phase labels explicitly.
VI. LATTICE-SWITCH MC IMPLEMENTATION DETAILS
Our MC procedure involves stochastic sampling of three
classes of variable: the displacement vectors
$
uW
%
; the volume
V; and the set of lattice vectors
$
RW
%
a
indexed, collectively, by
the lattice label a . The sampling of
$
uW
%
and V follows stan-
dard prescriptions @6#, with the effective configurational en-
ergy ~controlling metropolis acceptance probabilities! given
by Eq. ~25!. Stochastic updating of the lattice label a ~the
lattice switch! is discussed extensively in @2#. We summarize
the essentials.
There are many ways of defining the switch transforma-
tion ~lattice-to-lattice mapping!; for the fcc-hcp problem the
prescription illustrated in Fig. 4 is simple and efficient.
In the cases in which the two phases related by the switch
have significantly different volumes the switch operation can
accommodate an appropriate dilation @40#. In this case we
did not find that necessary: thus, as implemented here, the
switch transformation preserves both V and
$
uW
%
@41,42#.
With the prescription adopted in Eq. ~23! the switch also
leaves M ~and thus h
@
M
#
) unchanged @43#. The switch
FIG. 4. The LS transformation applied to the ground-state con-
figuration. The diagram shows six close-packed (x-y) layers. ~The
additional bracketed layer at the bottom is the periodic image of the
layer at the top.! The circles show the boundaries of particles lo-
cated at the sites of the two close-packed structures. In this realiza-
tion of the fcc→hcp lattice switch, the top pair of planes are left
unaltered, while the other pairs of planes are relocated by transla-
tions, specified by the black and white arrows.0-5
ture will typically have as its conjugate ~i.e., will map into! a
relatively high-energy configuration of the other structure.
The lattice switch MC step will typically be rejected. To
make it work the sampling of the configurations of each
structure must be multicanonically biased @3# to enhance the
probabilities of those states ~we call them gateway states!
from which a switch can be attempted with a reasonable
likelihood of success. The bias needs to reflect the difference
between the excitation energies @35# of conjugate pairs. To
that end we define an order parameter @36#
M[M
~
$
uW
%
!5E~
$
uW
%
,
$
RW
%hcp!2E~$uW %,$RW %fcc!, ~23!
where
E
~
$
uW
%
,
$
RW
%
a
!5b@
E
~
$
rW
%
!2E~
$
RW
%
a
!#
, ~24!
measures the excitation energy of configuration
$
rW
%
[
$
uW
%
,a
~the amount by which it exceeds the ground-state energy of
structure a) in units of kT5b21. The energy difference de-
fined in Eq. ~23! is the analog of the ‘‘overlap order param-
eter’’ utilized in our studies of hard spheres @1,2#.
Given the relationship between the energy of a configura-
tion and the energy of its conjugate one sees that displace-
ment patterns
$
uW
%
drawn from equilibrium sampling of the
fcc ~hcp! configuration space will give predominantly posi-
tive ~negative! values of
$
M
%
. The gateway states lie in be-
tween @37#; in the thermodynamic limit they have vanish-
ingly small equilibrium probability. We thus require to
engineer a multicanonical sampling algorithm to enhance the
probability along a notional line in M space, extending from
the ‘‘equilibrium’’ M values ~whose magnitude reflects the
energy cost of a LS acting on a typical configuration!
through to the small-M gateway configurations. This aim is
realized by augmenting the system energy function Eq. ~6!
appropriately,
bE
~
$
rW
%
!→bE~
$
rW
%
!1h@
M
#
[E˜
~
$
rW
%
!. ~25!
Here h
@
M
#
is a multicanonical weight function @3# that
must be chosen so as to allow the system to access the gate-
way configurations. For this purpose it is sufficient ~though
not necessarily optimal! to construct the weights such that
the multicanonical distribution P(MuN ,T ,
$
h
%
. . . ) @38# is
roughly flat. The desired ratio of partition functions @Eq. ~12!
or ~15!#, is identified with the ratio of the aggregate canoni-
cal probabilities of the two phases
R
aa
˜
~
N ,T , . . . ![
P
~
auN ,T , . . . !
P
~
a
˜
uN ,T , . . . !
5
(
$
M
%
P
~
M,auN ,T , . . . !
(
$
M
%
P
~
M,a˜ uN ,T , . . . !
~26!
and can be obtained from the corresponding ~measured! mul-
ticanonical distributions by reweighting ~folding out the
weights!03671(
$
M
%
P
~
M,auN ,T , . . . !
(
$
M
%
P
~
M,a˜ uN ,T , . . . !
5
(
$
M
%
eh[M]P
~
M,auN ,T ,
$
h
%
, . . . !
(
$
M
%
eh[M]P
~
M,a˜ uN ,T ,
$
h
%
, . . . !
. ~27!
The last step removes the explicit bias resulting from the
weights, leaving only the desired legacy @39#—that both
phases are visited with the canonical probabilities. Notice
that, in contrast to the hard-sphere problem, where M serves
as a unique identifier of the phase, the relevant sums in Eqs.
~26! and ~27! involve the phase labels explicitly.
VI. LATTICE-SWITCH MC IMPLEMENTATION DETAILS
Our MC procedure involves stochastic sampling of three
classes of variable: the displacement vectors
$
uW
%
; the volume
V; and the set of lattice vectors
$
RW
%
a
indexed, collectively, by
the lattice label a . The sampling of
$
uW
%
and V follows stan-
dard prescriptions @6#, with the effective configurational en-
ergy ~controlling metropolis acceptance probabilities! given
by Eq. ~25!. Stochastic updating of the lattice label a ~the
lattice switch! is discussed extensively in @2#. We summarize
the essentials.
There are many ways of defining the switch transforma-
tion ~lattice-to-lattice mapping!; for the fcc-hcp problem the
prescription illustrated in Fig. 4 is simple and efficient.
In the cases in which the two phases related by the switch
have significantly different volumes the switch operation can
accommodate an appropriate dilation @40#. In this case we
did not find that necessary: thus, as implemented here, the
switch transformation preserves both V and
$
uW
%
@41,42#.
With the prescription adopted in Eq. ~23! the switch also
leaves M ~and thus h
@
M
#
) unchanged @43#. The switch
FIG. 4. The LS transformation applied to the ground-state con-
figuration. The diagram shows six close-packed (x-y) layers. ~The
additional bracketed layer at the bottom is the periodic image of the
layer at the top.! The circles show the boundaries of particles lo-
cated at the sites of the two close-packed structures. In this realiza-
tion of the fcc→hcp lattice switch, the top pair of planes are left
unaltered, while the other pairs of planes are relocated by transla-
tions, specified by the black and white arrows.0-5
Page 6
A. N. JACKSON, A. D. BRUCE, AND G. J. ACKLAND PHYSICAL REVIEW E 65 036710constructed such that the associated weighted distribution
P(Muh
@
M
#
) is effectively ‘‘flat’’ ~the multicanonical ideal
@3#!. To the extent that it is flat one can regard the weight
function as a logarithmic measure of the unweighted ~equi-
librium! distribution P(M). The peaks in the weight func-
tions thus locate the ‘‘equilibrium’’ values of M for each
phase, the peaks at positive ~negative! M corresponding to
fcc ~hcp! phases respectively. The heights of these peaks
measure, on a logarithmic scale, the size of the free-energy
barrier to be negotiated for an escape from that phase, along
the route offered by LS @46# .
Figure 5~a! compares the results for the soft LJ potential
of primary interest here, with the results for a system of the
same number of particles interacting as hard spheres @2#. The
thermodynamic parameters were chosen so that the two sys-
FIG. 6. The multicanonical weights h ~solid line! and the prob-
ability of acceptance of a switch pacc as a function of M in the
region near M50. The system parameters are N563, p50, and
kT/e50.285.acceptance probability is then simply
Pa~a→a˜ !5min@1,exp~2DE˜ !# , ~28!
where @44#
DE˜ [b
@
E
~
$
uW
%
,
$
RW
%
a
˜
!2E~
$
uW
%
,
$
RW
%
a
!#
, ~29!
is the energy cost of the switch.
The practical core of the whole procedure lies in the de-
termination of the multicanonical weight function h
@
M
#
. In
generic terms the task is simple: the weight function is pa-
rameterized by a discrete set of weights each associated with
an interval of M space of width m0 ~that we set to a constant
value, m051) and together extending across the relevant
range of M space ~the values appropriate for the equilibrium
structures!; the weights are refined iteratively until the sam-
pling distribution is sufficiently close to the multicanonical
ideal. Details of the techniques we used are described else-
where @45#.
VII. LATTICE-SWITCH MC RESULTS
A. Features of the LS operation
We begin by presenting a selection of results chosen to
illustrate features of the LS operation itself; results for the
quantities of more immediate physical interest appear in the
subsections that follow.
Figure 5 provide examples of weight functions h
@
M
#03671tems have essentially the same densities. It is clear that
‘‘softening’’ the potential greatly reduces the barriers against
switching @47#.
Figure 5~b! shows the weights for the ~zero-pressure! LJ
system for a range of temperatures. The ‘‘barrier heights’’
grow with increasing temperature reflecting the growing en-
ergy mismatch between typical configurations and their con-
jugates.
The value of M for which h
@
M
#
assumes its minimum
~zero @46#! value—corresponding to the macrostates of low-
est probability in the equilibrium ensemble—identifies the
gateway macrostates. This identification is apparent opera-
tionally in Fig. 6 that shows that the probability of a switch
being accepted is maximal in the region where h
@
M
#
has its
minimum. Though this identification is not self-evidently
true, it is not surprising: the defining constraint on the ‘‘gate-
way’’ macrostates ~namely that they should allow a switch! is
sufficiently tight to ensure that they have the lowest probabil-
ity of all the macrostates that it is necessary to sample.
Close scrutiny of Fig. 5~b! shows that the M value MG ,
locating the gateway states is nonzero and shifts with tem-
perature. To understand this one should recall that the gate-
FIG. 5. Examples of ‘‘best-estimate’’ weight
functions for systems of N563 particles. ~a! The
Lennard-Jones ~LJ! weight function ~for p
50,kT/e50.285) is compared with the weight
function for hard spheres ~HS; see @2#! at a pres-
sure chosen such that the densities of the two
systems are comparable. ~b! Weight functions for
the p50 LJ system at three different tempera-
tures.0-6
P(Muh
@
M
#
) is effectively ‘‘flat’’ ~the multicanonical ideal
@3#!. To the extent that it is flat one can regard the weight
function as a logarithmic measure of the unweighted ~equi-
librium! distribution P(M). The peaks in the weight func-
tions thus locate the ‘‘equilibrium’’ values of M for each
phase, the peaks at positive ~negative! M corresponding to
fcc ~hcp! phases respectively. The heights of these peaks
measure, on a logarithmic scale, the size of the free-energy
barrier to be negotiated for an escape from that phase, along
the route offered by LS @46# .
Figure 5~a! compares the results for the soft LJ potential
of primary interest here, with the results for a system of the
same number of particles interacting as hard spheres @2#. The
thermodynamic parameters were chosen so that the two sys-
FIG. 6. The multicanonical weights h ~solid line! and the prob-
ability of acceptance of a switch pacc as a function of M in the
region near M50. The system parameters are N563, p50, and
kT/e50.285.acceptance probability is then simply
Pa~a→a˜ !5min@1,exp~2DE˜ !# , ~28!
where @44#
DE˜ [b
@
E
~
$
uW
%
,
$
RW
%
a
˜
!2E~
$
uW
%
,
$
RW
%
a
!#
, ~29!
is the energy cost of the switch.
The practical core of the whole procedure lies in the de-
termination of the multicanonical weight function h
@
M
#
. In
generic terms the task is simple: the weight function is pa-
rameterized by a discrete set of weights each associated with
an interval of M space of width m0 ~that we set to a constant
value, m051) and together extending across the relevant
range of M space ~the values appropriate for the equilibrium
structures!; the weights are refined iteratively until the sam-
pling distribution is sufficiently close to the multicanonical
ideal. Details of the techniques we used are described else-
where @45#.
VII. LATTICE-SWITCH MC RESULTS
A. Features of the LS operation
We begin by presenting a selection of results chosen to
illustrate features of the LS operation itself; results for the
quantities of more immediate physical interest appear in the
subsections that follow.
Figure 5 provide examples of weight functions h
@
M
#03671tems have essentially the same densities. It is clear that
‘‘softening’’ the potential greatly reduces the barriers against
switching @47#.
Figure 5~b! shows the weights for the ~zero-pressure! LJ
system for a range of temperatures. The ‘‘barrier heights’’
grow with increasing temperature reflecting the growing en-
ergy mismatch between typical configurations and their con-
jugates.
The value of M for which h
@
M
#
assumes its minimum
~zero @46#! value—corresponding to the macrostates of low-
est probability in the equilibrium ensemble—identifies the
gateway macrostates. This identification is apparent opera-
tionally in Fig. 6 that shows that the probability of a switch
being accepted is maximal in the region where h
@
M
#
has its
minimum. Though this identification is not self-evidently
true, it is not surprising: the defining constraint on the ‘‘gate-
way’’ macrostates ~namely that they should allow a switch! is
sufficiently tight to ensure that they have the lowest probabil-
ity of all the macrostates that it is necessary to sample.
Close scrutiny of Fig. 5~b! shows that the M value MG ,
locating the gateway states is nonzero and shifts with tem-
perature. To understand this one should recall that the gate-
FIG. 5. Examples of ‘‘best-estimate’’ weight
functions for systems of N563 particles. ~a! The
Lennard-Jones ~LJ! weight function ~for p
50,kT/e50.285) is compared with the weight
function for hard spheres ~HS; see @2#! at a pres-
sure chosen such that the densities of the two
systems are comparable. ~b! Weight functions for
the p50 LJ system at three different tempera-
tures.0-6
Page 7
LATTICE-SWITCH MONTE CARLO METHOD: . . . PHYSICAL REVIEW E 65 036710distributions for the appropriate state point and applying Eq.
~11!. Evidently, the two methods agree at low enough tem-
peratures, where the harmonic theory is trustworthy; this pro-
vides a reassuring consistency check. As we have already
noted, the results of Ref. @34# provide consistency checks on
the harmonic studies themselves.
The differences between the results of the two methods
should, then simply reflect the anharmonic effects that are
captured by the LS studies. The differences grow with in-
FIG. 8. The difference @21# between the Helmholtz free-energy
densities of the two structures for a system of N5123 particles at a
density of ~a! rs351.092 close to the zero-pressure density and ~b!
rs
3
51.538 the density at which De fcc, hcp0 is maximal. In each case
the solid line represents the results of a harmonic calculation, which
follows Eq. ~30! with intercepts and gradients given by
De
0/e50.000 869 8
@
0.001 152
#
and D f h/kT520.002 574
@
20.002 949
#
for cases ~a! @~b!# respectively. The data points are
the results of LSMC studies utilizing Eq. ~11!.way states are such that the switch energy cost DE˜ @Eq. ~29!#
is close to zero @48#. Recalling the definition of M we see
that the condition DE˜ 50 implies M5MG5NbDe fcc, hcp0 :
the behavior of MG can thus be traced to the difference in
ground-state energies.
Figure 7 shows the equilibrium M distributions for the
zero-pressure LJ system at the three different temperatures
featuring in Fig. 5~b!: in each case we sample the multica-
nonical distributions defined by the appropriate set of
weights ~allowing both phases to be sampled! and then re-
weight to obtain the equilibrium distribution. Equation ~15!
shows that the relative stability of the two phases is directly
reflected in the ratio of the integrated areas of the peaks in
the equilibrium M distribution. We thus see immediately
and transparently from Fig. 7 that, while the hcp structure is
favored at low enough temperature ~as implied by the differ-
ence between the ground-state energies at this pressure, Fig.
1!, the favored phase at high enough temperatures is fcc.
B. Benchmarking: LSMC and the harmonic approximation
A comparison between the results based on harmonic
crystal dynamics ~Sec. IV! and LSMC is worthwhile for two
reasons. In the regime of low enough T, the harmonic theory
is necessarily correct: comparison with LSMC results then
provides a check on the overall consistency of the LS frame-
work. In the regime of higher T, we must expect the har-
monic theory to be growingly inaccurate: comparison with
LSMC allows us to explore the extent of these inaccuracies.
Figure 8 shows the results of LSMC studies in the NVT
ensemble, at two different densities, for a range of tempera-
tures. In each case the solid line is the free-energy difference
in the harmonic approximation, which combines the ground-
state-energy difference @Eq. ~5!# and the harmonic contribu-
tion to the free energy @Eq. ~22!#:
D f
e
5
De
0
e
1
D f h
kT
kT
e
~
harmonic!. ~30!
The LS values follow from simulations determining the M03671creasing temperature and decrease with increasing density;
this behavior is consistent with simple qualitative expecta-
tions @49#.
At both densities the deviations from harmonic behavior
are such as to favor hcp: we see no change in the sign of the
‘‘anharmonic’’ contributions, with increasing density, sug-
gested by the perturbation-theory calculations in @18#.
At kT/e50.4 and rs351.092, we find D f /e
520.000 074(16), which is virtually indistinguishable from
the predictions of the harmonic theory. In contrast Ref. @50#
~which uses an integration method to determine the free en-
ergies of each phase separately! reports D f /e520.066(4) at
kT/e50.4 and rs351; this seems irreconcilable with our
results.
At the temperatures of most immediate physical interest
~where D f 50) the discrepancies between the two methods
are small, even for the zero-pressure density. We may thus
expect that the harmonic theory of the fcc-hcp phase bound-
ary will not prove far wrong; this is indeed the case.
C. The fcc-hcp phase boundary
We have established that LSMC provides us with a way
of determining the difference between the free energies of
the two phases at some point in the space of thermodynamic
coordinates. We must now turn to the practical task of inter-
est here: the determination of the phase behavior itself. This
task may be divided into two—extrapolation and tracking.
First ~‘‘extrapolation’’! we need to use the data accumulated
FIG. 7. Temperature dependence of the prob-
ability distribution P(M) for a system of N
563 Lennard-Jones particles at zero pressure.0-7
~11!. Evidently, the two methods agree at low enough tem-
peratures, where the harmonic theory is trustworthy; this pro-
vides a reassuring consistency check. As we have already
noted, the results of Ref. @34# provide consistency checks on
the harmonic studies themselves.
The differences between the results of the two methods
should, then simply reflect the anharmonic effects that are
captured by the LS studies. The differences grow with in-
FIG. 8. The difference @21# between the Helmholtz free-energy
densities of the two structures for a system of N5123 particles at a
density of ~a! rs351.092 close to the zero-pressure density and ~b!
rs
3
51.538 the density at which De fcc, hcp0 is maximal. In each case
the solid line represents the results of a harmonic calculation, which
follows Eq. ~30! with intercepts and gradients given by
De
0/e50.000 869 8
@
0.001 152
#
and D f h/kT520.002 574
@
20.002 949
#
for cases ~a! @~b!# respectively. The data points are
the results of LSMC studies utilizing Eq. ~11!.way states are such that the switch energy cost DE˜ @Eq. ~29!#
is close to zero @48#. Recalling the definition of M we see
that the condition DE˜ 50 implies M5MG5NbDe fcc, hcp0 :
the behavior of MG can thus be traced to the difference in
ground-state energies.
Figure 7 shows the equilibrium M distributions for the
zero-pressure LJ system at the three different temperatures
featuring in Fig. 5~b!: in each case we sample the multica-
nonical distributions defined by the appropriate set of
weights ~allowing both phases to be sampled! and then re-
weight to obtain the equilibrium distribution. Equation ~15!
shows that the relative stability of the two phases is directly
reflected in the ratio of the integrated areas of the peaks in
the equilibrium M distribution. We thus see immediately
and transparently from Fig. 7 that, while the hcp structure is
favored at low enough temperature ~as implied by the differ-
ence between the ground-state energies at this pressure, Fig.
1!, the favored phase at high enough temperatures is fcc.
B. Benchmarking: LSMC and the harmonic approximation
A comparison between the results based on harmonic
crystal dynamics ~Sec. IV! and LSMC is worthwhile for two
reasons. In the regime of low enough T, the harmonic theory
is necessarily correct: comparison with LSMC results then
provides a check on the overall consistency of the LS frame-
work. In the regime of higher T, we must expect the har-
monic theory to be growingly inaccurate: comparison with
LSMC allows us to explore the extent of these inaccuracies.
Figure 8 shows the results of LSMC studies in the NVT
ensemble, at two different densities, for a range of tempera-
tures. In each case the solid line is the free-energy difference
in the harmonic approximation, which combines the ground-
state-energy difference @Eq. ~5!# and the harmonic contribu-
tion to the free energy @Eq. ~22!#:
D f
e
5
De
0
e
1
D f h
kT
kT
e
~
harmonic!. ~30!
The LS values follow from simulations determining the M03671creasing temperature and decrease with increasing density;
this behavior is consistent with simple qualitative expecta-
tions @49#.
At both densities the deviations from harmonic behavior
are such as to favor hcp: we see no change in the sign of the
‘‘anharmonic’’ contributions, with increasing density, sug-
gested by the perturbation-theory calculations in @18#.
At kT/e50.4 and rs351.092, we find D f /e
520.000 074(16), which is virtually indistinguishable from
the predictions of the harmonic theory. In contrast Ref. @50#
~which uses an integration method to determine the free en-
ergies of each phase separately! reports D f /e520.066(4) at
kT/e50.4 and rs351; this seems irreconcilable with our
results.
At the temperatures of most immediate physical interest
~where D f 50) the discrepancies between the two methods
are small, even for the zero-pressure density. We may thus
expect that the harmonic theory of the fcc-hcp phase bound-
ary will not prove far wrong; this is indeed the case.
C. The fcc-hcp phase boundary
We have established that LSMC provides us with a way
of determining the difference between the free energies of
the two phases at some point in the space of thermodynamic
coordinates. We must now turn to the practical task of inter-
est here: the determination of the phase behavior itself. This
task may be divided into two—extrapolation and tracking.
First ~‘‘extrapolation’’! we need to use the data accumulated
FIG. 7. Temperature dependence of the prob-
ability distribution P(M) for a system of N
563 Lennard-Jones particles at zero pressure.0-7
Page 8
A. N. JACKSON, A. D. BRUCE, AND G. J. ACKLAND PHYSICAL REVIEW E 65 036710at our initially chosen state point to infer the location of some
point on the phase boundary ~and refine it as necessary!.
Second ~‘‘tracking’’! we need to map out the phase boundary
emanating from that point. There are many ways of address-
ing these tasks; we shall focus on the strategies we have
found most useful, while noting other approaches that proved
less satisfactory here.
Both extrapolation and tracking stages make use of the
fact that the derivatives of the free-energy difference can be
expressed in terms of the differences between single-phase
averages of appropriate observables
]
@
bDg
#
]b
5
1
N ^DE1pDV&5Dh ~31a!
and
]
@
bDg
#
]p 5
1
N ^DV&5Dv , ~31b!
where Dh and Dv represent the differences between the en-
thalpy per particle and mean volume per particle, for the two
phases. Since the densities of fcc and hcp phases are very
similar at coexistence @51# the p derivative @Eq. ~31b!# is
hard to determine accurately. But the b derivative @Eq. ~31a!#
is readily measurable to the required precision, using the
appropriate single-phase averages already available in the LS
simulations used to determine Dg itself. Given Dg and its b
derivative a simple linear extrapolation ~LE! scheme pro-
vides an estimate of the phase boundary temperature for the
chosen pressure. Figure 9 provides an example. The extrapo-
lation works well—a reflection of the fact that the behavior
is close to harmonic over the region involved. Figure 9 also
shows the results generated by a single histogram extrapola-
tion ~HE! @52# in which the M distribution measured
through the histogram accumulated at the initial state point is
used to estimate the distribution at other temperatures. In this
case HE fares much worse than LE. The reason is that the
HE predictions for some extrapolated temperature b explic-
FIG. 9. Comparison of linear extrapolation ~LE! and histogram
extrapolation ~HE! methods for determining the p50 hcp-fcc phase
boundary temperature T0(p50) on the basis of simulation data
accumulated at kT/e50.1 and p50 for a system of N563
Lennard-Jones particles. Our best estimate of T0(p50) is also
shown.03671itly utilize the statistical information that the simulations
~performed at b0) give about the dominant configurations at
b; this information typically resides in the poorly sampled
wings of the measured distribution. By contrast LE implicitly
models the measured distribution ~in this case! as two Gaus-
sians, each centered on the average for the phase, whose
characteristics are then determined by the well-sampled re-
gions of the measured distribution. LE does well when ~as
here! the near-Gaussian assumption is good.
Once a point on the phase boundary is established the rest
of the boundary can, in principle, be mapped by numerical
integration of the phase boundary derivative that follows
from Eqs. ~31a! and ~31b!. This strategy—Gibbs-Duhem
~GD! integration—has been widely and effectively used @53#.
In this context, however, it is undermined by the problem we
have already noted with the measurement of the volume dif-
ference at coexistence. The resulting inaccuracies are com-
pounded by the generic draw back of GD—the absence of
any self-correction mechanism. Figure 10 shows the results
of a low-order predict-evaluate-correct ~PEC! GD integra-
tion: on its own GD is clearly not suited to the present prob-
lem. Instead we chose a strategy which, while making use of
GD to the extent it is trustworthy, employs LS to correct it.
Starting at some coexistence point p1 ,b1 say we choose a
different pressure p25p11Dp . We use the GD derivative to
predict the associated transition temperature b2
e
. If the dif-
ference between b2
e and b1 is not statistically significant
~given the errors associated with the GD derivative! we reset
b2
e
5b1. The multicanonical weights used at b1 ,p1 provide
estimates of the weights at b2
e
,p2 that can be easily refined.
The process is completed by a b extrapolation to the phase
boundary ~as discussed above!. Sustaining switching be-
tween phases as the phase boundary is traced out delivers a
running consistency check. In this case the computational
cost is actually less than it would be to determine the GD
derivative to the accuracy required to make GD itself viable.
Figure 10 shows some representative results.
Using these techniques we have mapped out the hcp-fcc
phase boundary from the zero-pressure limit, through to the
maximum pressure at which hcp is classically stable at T
50. The results are gathered together in Fig. 11. In the case
FIG. 10. Comparison of lattice-switch ~LS! and Gibbs-Duhem
~GD! methods for mapping out a high-pressure section of the phase
boundary for a system of N563 Lennard-Jones particles. The re-
sults of harmonic theory are also shown.0-8
point on the phase boundary ~and refine it as necessary!.
Second ~‘‘tracking’’! we need to map out the phase boundary
emanating from that point. There are many ways of address-
ing these tasks; we shall focus on the strategies we have
found most useful, while noting other approaches that proved
less satisfactory here.
Both extrapolation and tracking stages make use of the
fact that the derivatives of the free-energy difference can be
expressed in terms of the differences between single-phase
averages of appropriate observables
]
@
bDg
#
]b
5
1
N ^DE1pDV&5Dh ~31a!
and
]
@
bDg
#
]p 5
1
N ^DV&5Dv , ~31b!
where Dh and Dv represent the differences between the en-
thalpy per particle and mean volume per particle, for the two
phases. Since the densities of fcc and hcp phases are very
similar at coexistence @51# the p derivative @Eq. ~31b!# is
hard to determine accurately. But the b derivative @Eq. ~31a!#
is readily measurable to the required precision, using the
appropriate single-phase averages already available in the LS
simulations used to determine Dg itself. Given Dg and its b
derivative a simple linear extrapolation ~LE! scheme pro-
vides an estimate of the phase boundary temperature for the
chosen pressure. Figure 9 provides an example. The extrapo-
lation works well—a reflection of the fact that the behavior
is close to harmonic over the region involved. Figure 9 also
shows the results generated by a single histogram extrapola-
tion ~HE! @52# in which the M distribution measured
through the histogram accumulated at the initial state point is
used to estimate the distribution at other temperatures. In this
case HE fares much worse than LE. The reason is that the
HE predictions for some extrapolated temperature b explic-
FIG. 9. Comparison of linear extrapolation ~LE! and histogram
extrapolation ~HE! methods for determining the p50 hcp-fcc phase
boundary temperature T0(p50) on the basis of simulation data
accumulated at kT/e50.1 and p50 for a system of N563
Lennard-Jones particles. Our best estimate of T0(p50) is also
shown.03671itly utilize the statistical information that the simulations
~performed at b0) give about the dominant configurations at
b; this information typically resides in the poorly sampled
wings of the measured distribution. By contrast LE implicitly
models the measured distribution ~in this case! as two Gaus-
sians, each centered on the average for the phase, whose
characteristics are then determined by the well-sampled re-
gions of the measured distribution. LE does well when ~as
here! the near-Gaussian assumption is good.
Once a point on the phase boundary is established the rest
of the boundary can, in principle, be mapped by numerical
integration of the phase boundary derivative that follows
from Eqs. ~31a! and ~31b!. This strategy—Gibbs-Duhem
~GD! integration—has been widely and effectively used @53#.
In this context, however, it is undermined by the problem we
have already noted with the measurement of the volume dif-
ference at coexistence. The resulting inaccuracies are com-
pounded by the generic draw back of GD—the absence of
any self-correction mechanism. Figure 10 shows the results
of a low-order predict-evaluate-correct ~PEC! GD integra-
tion: on its own GD is clearly not suited to the present prob-
lem. Instead we chose a strategy which, while making use of
GD to the extent it is trustworthy, employs LS to correct it.
Starting at some coexistence point p1 ,b1 say we choose a
different pressure p25p11Dp . We use the GD derivative to
predict the associated transition temperature b2
e
. If the dif-
ference between b2
e and b1 is not statistically significant
~given the errors associated with the GD derivative! we reset
b2
e
5b1. The multicanonical weights used at b1 ,p1 provide
estimates of the weights at b2
e
,p2 that can be easily refined.
The process is completed by a b extrapolation to the phase
boundary ~as discussed above!. Sustaining switching be-
tween phases as the phase boundary is traced out delivers a
running consistency check. In this case the computational
cost is actually less than it would be to determine the GD
derivative to the accuracy required to make GD itself viable.
Figure 10 shows some representative results.
Using these techniques we have mapped out the hcp-fcc
phase boundary from the zero-pressure limit, through to the
maximum pressure at which hcp is classically stable at T
50. The results are gathered together in Fig. 11. In the case
FIG. 10. Comparison of lattice-switch ~LS! and Gibbs-Duhem
~GD! methods for mapping out a high-pressure section of the phase
boundary for a system of N563 Lennard-Jones particles. The re-
sults of harmonic theory are also shown.0-8
Page 9
LATTICE-SWITCH MONTE CARLO METHOD: . . . PHYSICAL REVIEW E 65 036710of the smaller (N563) system we gathered LS data along
the entire phase boundary, demonstrating agreement with the
results of harmonic theory except in a low-pressure region
where the break-down of the harmonic approximation begins
to be apparent. In the case of the larger (N5123) system we
focused our LS studies on this low pressure region, mapping
out the phase boundary from p50 through to the pressure at
which it becomes indistinguishable from the predictions of
the harmonic theory. The dashed line through the LS points
is based on a simple phenomenology of the growth of anhar-
monic effects with decreasing density @49,54#. While the de-
viations from the harmonic (NVT) theory are certainly re-
solvable, they are always small. And indeed the same general
structure and scale of the phase diagram is apparent in pio-
neering ~if semiquantitative! harmonic lattice dynamics stud-
ies of 50 years ago @24#. It is harder to make useful connec-
tions with more recent studies. Jackson and Swol @55# and
Choi et al. @17# both appeal to perturbation theories using a
truncated potential giving a complex phase behavior. In @18#
the same method is applied to a system without truncation;
but the perturbation theory presupposes hard-sphere results
now known to be incorrect, and is in any case restricted to
high temperatures (kT/e*1).
VIII. DISCUSSION
We divide this concluding section into three. First we
shall take stock of where we stand on the generic problem
FIG. 11. A variety of approximations to the classical Lennard-
Jones phase diagram. In all cases the difference between the
ground-state energies of the two structures is essentially exact; the
fluctuation contribution to the energy is truncated according to Eq.
~7!, with rx51.5r0. The dashed and solid lines are the results of
harmonic calculations ~for the two system sizes!. The dash-dotted
line is a phenomenological parameterization of the anharmonic ef-
fects @54#. The scale at the top of the figure shows the pressures at
selected points on the ~LS N5123, NPT) coexistence curve. Tie-
line structure is unresolvable on the scale of the figure.03671addressed here: the development of a general computational
technique for the study of the phase behavior of crystalline
matter. Then we shall summarize what we have learned
about the particular problem we have chosen to address: the
equilibrium phase behavior of the LJ model. Finally we shall
consider the implications of the LJ model predictions for
studies of the crystallization process itself.
1. The technique
To assess the LS technique we must now set it in context.
The task of determining the phase behavior of solids has
traditionally been seen as requiring separate calculations of
the absolute free energies of the competing phases, linking
each to some reference system. This approach was pioneered
by Frenkel and Ladd @5# and has been implemented in a
variety of forms, differing according to the choice of refer-
ence system ~Einstein solid, harmonic solid, ideal gas, . . . !
and how the linking path is negotiated ~integration, extended
sampling! @50,56–59#. The problem facing this strategy is
that the quantity of physical interest—the free-energy differ-
ence between the phases—is typically small ~in the vicinity
of a phase boundary necessarily so! compared to the absolute
free energies of the separate phases; one is left with the task
of determining the difference between two ‘‘relatively large’’
numbers @60# that are themselves of little or no physical in-
terest. It is possible to make this strategy deliver results of
comparable accuracy to that of LS, for comparable compu-
tational resource @61#. But the extent of the signal-to-noise
problem provides strong motivation for developing a more
direct strategy.
The essence of a ‘‘direct strategy’’ is that it should focus
on the free-energy difference between ~relative likelihood of!
the two phases: this requires a ‘‘path’’ connecting the disjoint
configuration spaces of the two phases. The pioneering study
here is that of Rahman and Jacucci @62#. Subsequent studies
in this vein can be differentiated according to whether ~in the
terminology of @62#! the path is ‘‘physical’’ or ‘‘nonphysical’’
and according to how the path is traversed @63–65#. Refer-
ence @65#, devoted to the hcp-fcc phase behavior of hard
spheres, provides the most recent instance of the use of a
‘‘physical’’ path between the two structures @66#, explored
with multicanonical sampling @3#. In that case the ‘‘physical
path’’ was found to encounter some stubbornly physical bar-
riers that ~the authors report! made it less efficient than a LS
strategy that they also explored.
The LS technique belongs to the class of ‘‘direct’’ meth-
ods that exploit a ‘‘nonphysical’’ path. One can see it @67# as
a reformulation of the strategy used by Moody et al. @63# in
early ~though numerically inconclusive! studies of polymor-
phic crystalline behavior. The principal additional feature of
LS is that we have used extended ~‘‘multicanonical’’! sam-
pling instead of multistage sampling; that change, combined
with the availability of massively more computational re-
source produces a technique that—we have seen here and in
@2#—delivers the precision one requires, with readily quanti-
fiable uncertainties, in a physically transparent fashion.
There are two noteworthy ways in which the LS frame-
work could be fruitfully developed. First, as noted in @2#, one
may fold into the switch operation any desired transforma-0-9
the entire phase boundary, demonstrating agreement with the
results of harmonic theory except in a low-pressure region
where the break-down of the harmonic approximation begins
to be apparent. In the case of the larger (N5123) system we
focused our LS studies on this low pressure region, mapping
out the phase boundary from p50 through to the pressure at
which it becomes indistinguishable from the predictions of
the harmonic theory. The dashed line through the LS points
is based on a simple phenomenology of the growth of anhar-
monic effects with decreasing density @49,54#. While the de-
viations from the harmonic (NVT) theory are certainly re-
solvable, they are always small. And indeed the same general
structure and scale of the phase diagram is apparent in pio-
neering ~if semiquantitative! harmonic lattice dynamics stud-
ies of 50 years ago @24#. It is harder to make useful connec-
tions with more recent studies. Jackson and Swol @55# and
Choi et al. @17# both appeal to perturbation theories using a
truncated potential giving a complex phase behavior. In @18#
the same method is applied to a system without truncation;
but the perturbation theory presupposes hard-sphere results
now known to be incorrect, and is in any case restricted to
high temperatures (kT/e*1).
VIII. DISCUSSION
We divide this concluding section into three. First we
shall take stock of where we stand on the generic problem
FIG. 11. A variety of approximations to the classical Lennard-
Jones phase diagram. In all cases the difference between the
ground-state energies of the two structures is essentially exact; the
fluctuation contribution to the energy is truncated according to Eq.
~7!, with rx51.5r0. The dashed and solid lines are the results of
harmonic calculations ~for the two system sizes!. The dash-dotted
line is a phenomenological parameterization of the anharmonic ef-
fects @54#. The scale at the top of the figure shows the pressures at
selected points on the ~LS N5123, NPT) coexistence curve. Tie-
line structure is unresolvable on the scale of the figure.03671addressed here: the development of a general computational
technique for the study of the phase behavior of crystalline
matter. Then we shall summarize what we have learned
about the particular problem we have chosen to address: the
equilibrium phase behavior of the LJ model. Finally we shall
consider the implications of the LJ model predictions for
studies of the crystallization process itself.
1. The technique
To assess the LS technique we must now set it in context.
The task of determining the phase behavior of solids has
traditionally been seen as requiring separate calculations of
the absolute free energies of the competing phases, linking
each to some reference system. This approach was pioneered
by Frenkel and Ladd @5# and has been implemented in a
variety of forms, differing according to the choice of refer-
ence system ~Einstein solid, harmonic solid, ideal gas, . . . !
and how the linking path is negotiated ~integration, extended
sampling! @50,56–59#. The problem facing this strategy is
that the quantity of physical interest—the free-energy differ-
ence between the phases—is typically small ~in the vicinity
of a phase boundary necessarily so! compared to the absolute
free energies of the separate phases; one is left with the task
of determining the difference between two ‘‘relatively large’’
numbers @60# that are themselves of little or no physical in-
terest. It is possible to make this strategy deliver results of
comparable accuracy to that of LS, for comparable compu-
tational resource @61#. But the extent of the signal-to-noise
problem provides strong motivation for developing a more
direct strategy.
The essence of a ‘‘direct strategy’’ is that it should focus
on the free-energy difference between ~relative likelihood of!
the two phases: this requires a ‘‘path’’ connecting the disjoint
configuration spaces of the two phases. The pioneering study
here is that of Rahman and Jacucci @62#. Subsequent studies
in this vein can be differentiated according to whether ~in the
terminology of @62#! the path is ‘‘physical’’ or ‘‘nonphysical’’
and according to how the path is traversed @63–65#. Refer-
ence @65#, devoted to the hcp-fcc phase behavior of hard
spheres, provides the most recent instance of the use of a
‘‘physical’’ path between the two structures @66#, explored
with multicanonical sampling @3#. In that case the ‘‘physical
path’’ was found to encounter some stubbornly physical bar-
riers that ~the authors report! made it less efficient than a LS
strategy that they also explored.
The LS technique belongs to the class of ‘‘direct’’ meth-
ods that exploit a ‘‘nonphysical’’ path. One can see it @67# as
a reformulation of the strategy used by Moody et al. @63# in
early ~though numerically inconclusive! studies of polymor-
phic crystalline behavior. The principal additional feature of
LS is that we have used extended ~‘‘multicanonical’’! sam-
pling instead of multistage sampling; that change, combined
with the availability of massively more computational re-
source produces a technique that—we have seen here and in
@2#—delivers the precision one requires, with readily quanti-
fiable uncertainties, in a physically transparent fashion.
There are two noteworthy ways in which the LS frame-
work could be fruitfully developed. First, as noted in @2#, one
may fold into the switch operation any desired transforma-0-9
Page 10
chose for this exploratory study. From some standpoints that
choice proved ill advised; while LJ is regarded as the generic
soft potential, its phase behavior is acutely sensitive to de-
in LJ systems. There have been many studies of this kind.
Much of the focus has been on the initial stages of the
process—in particular, the suggestion @72# that the route to
A. N. JACKSON, A. D. BRUCE, AND G. J. ACKLAND PHYSICAL REVIEW E 65 036710tails ~notably the range of interactions permitted! that are
largely arbitrary adjuncts of the phenomenology. The diver-
sity of LJ phase diagrams featuring in the literature is testi-
mony to this sensitivity. Should one wish to take the LJ
potential at face value and establish a definitive phase dia-
gram one must certainly investigate the effects of potential-
range truncation and ~other! finite-size effects systematically.
We have attempted to do so here: we have dealt exactly with
the ground-state energies; we have seen that our truncation
scheme for the excitation energy captures the thermodynami-
cally limiting behavior well in the harmonic regime ~Figs 2,
3!; and we find no evidence to doubt it when one moves
beyond the harmonic region @70#. We have some confidence
then that Fig. 11 does indeed capture the classical LJ crys-
talline phase behavior.
@1# A. D. Bruce, N. B. Wilding, and G. J. Ackland, Phys. Rev.
Lett. 79, 3002 ~1997!.
@2# A. D. Bruce, A. N. Jackson, G. J. Ackland, and N. B. Wilding,
Phys. Rev. E 61, 906 ~2000!.
@3# B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 ~1992!; B.
A. Berg, J. Stat. Phys. 82, 323 ~1996!.
@4# M. P. Allen, in Monte Carlo and Molecular Dynamics of Con-
densed Matter Systems, edited by K. Binder and G. Ciccotti
~Italian Physical Society, Bologna, 1996!, p. 255.
@5# D. Frenkel and A. J. C. Ladd, J. Chem. Phys. 81, 3188 ~1984!.
@6# D. Frenkel and B. Smit, Understanding Molecular Simulation
~Academic, New York, 1996!.
@7# N. B. Wilding and A. D. Bruce, Phys. Rev. Lett. 85, 5138
~2000!.
@8# In the case of hard spheres the temperature enters only trivially
through the scale of the pressure.
@9# D. A. Young, Phase Diagrams of the Elements ~University of
California Press, Berkeley, 1991!.
@10# T. H. K. Barron and M. L. Klein, in Dynamical Properties of036710the thermodynamically favored phase traverses a region with
a bcc structure. This phenomenon ~for which there is simu-
lation support @73#! is a useful reminder of the importance of
the kinetic pathways by which the equilibrium state is
reached, in the long term. The evidence from the most de-
tailed study known to us @74# is that a fcc structure emerges
in that ‘‘long term’’ @75#. But this study utilized a LJ poten-
tial truncated at rx52.3s; and we have seen that such details
matter as regards the ‘‘true’’ equilibrium behavior. In any
case, given the LJ knife edge, it seems likely that the par-
ticular long term behavior that does emerge may well have
more to do with the kinetics than the thermodynamics. That
interplay is probably better explored in a system other than
LJ. Recent work has begun to address this interesting area
@76#.
Solids, edited by G. K. Horton and A. A. Maradudin ~North-
Holland, Amsterdam, 1974!, p. 391.
@11# J. A. Barker, in Rare Gas Solids, edited by M. L. Klein and J.
A. Venables ~Academic Press, London, 1976!, Vol. 1, p. 212.
@12# One should not expect it to work at sufficiently extreme pres-
sures and it does not: see for example, R. J. Hemley and N. W.
Ashcroft, Phys. Today 51 ~8!, 26 ~1998!.
@13# Helium, of course, is different, see e.g., H. R. Glyde, Rare Gas
Solids ~Ref. @11#!, p 382.
@14# See e.g., K. F. Niebel and J. A. Venables, in Rare Gas Solids
~Ref. @11#!, p. 558.
@15# B. J. Alder and P. H. Paulson, J. Chem. Phys. 43, 4172 ~1965!.
@16# Compare Refs. @17,18#; see also Fig. 1.
@17# Y. Choi, T. Ree, and F. H. Ree, J. Chem. Phys. 95, 7548
~1991!.
@18# Y. Choi, T. Ree, and F. H. Ree, J. Chem. Phys. 99, 9917
~1993!.
@19# We use the term ‘‘lattice’’ rather loosely. The set of vectors
$
RW
%tion of the displacement coordinates, provided one takes ap-
propriate account ~through the Jacobian! of the resulting
warping of the configuration space @68#. In particular one can
define a switch operation that preserves collective ~normal
mode! coordinates rather than local coordinates. In the har-
monic region such an operation will work without the need
for any weighting; the barrier ~to be surmounted by multica-
nonical methods! is reduced to that associated with the an-
harmonic contributions. This approach is currently under in-
vestigation.
Second, the current framework is strictly classical. It is
clearly desirable to extend the framework to include quan-
tum effects; we plan to do so within a path integral formal-
ism @69#.
2. The system
We turn now to the particular system ~the LJ model! we3. The implications
If the calculations reported here are to represent more
than an exercise in equilibrium statistical mechanics then
their results need to be related to studies of the crystallization
process itself. It is not easy to do so in any satisfying way.
If one tries to make contact with observations on RGS one
is faced with a whole range of effects that are potentially
important. We have already noted the need to extend the
framework to include quantum effects—clearly important for
the lighter RGS and essential in understanding the behavior
of Helium. But the link with real systems is also complicated
by two other manifestations of the fcc-hcp knife edge on
which the LJ model sits. That knife edge is presumably re-
flected in the sensitivity of the observed structure to impuri-
ties @71#; it also means that effects that are missing from the
two-body central force LJ phenomenology ~and that are not
crucial to other aspects of the RGS behavior! need to be
addressed @14#. As a more limited objective one may try to
make contact with simulations of the crystallization process-10
choice proved ill advised; while LJ is regarded as the generic
soft potential, its phase behavior is acutely sensitive to de-
in LJ systems. There have been many studies of this kind.
Much of the focus has been on the initial stages of the
process—in particular, the suggestion @72# that the route to
A. N. JACKSON, A. D. BRUCE, AND G. J. ACKLAND PHYSICAL REVIEW E 65 036710tails ~notably the range of interactions permitted! that are
largely arbitrary adjuncts of the phenomenology. The diver-
sity of LJ phase diagrams featuring in the literature is testi-
mony to this sensitivity. Should one wish to take the LJ
potential at face value and establish a definitive phase dia-
gram one must certainly investigate the effects of potential-
range truncation and ~other! finite-size effects systematically.
We have attempted to do so here: we have dealt exactly with
the ground-state energies; we have seen that our truncation
scheme for the excitation energy captures the thermodynami-
cally limiting behavior well in the harmonic regime ~Figs 2,
3!; and we find no evidence to doubt it when one moves
beyond the harmonic region @70#. We have some confidence
then that Fig. 11 does indeed capture the classical LJ crys-
talline phase behavior.
@1# A. D. Bruce, N. B. Wilding, and G. J. Ackland, Phys. Rev.
Lett. 79, 3002 ~1997!.
@2# A. D. Bruce, A. N. Jackson, G. J. Ackland, and N. B. Wilding,
Phys. Rev. E 61, 906 ~2000!.
@3# B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 ~1992!; B.
A. Berg, J. Stat. Phys. 82, 323 ~1996!.
@4# M. P. Allen, in Monte Carlo and Molecular Dynamics of Con-
densed Matter Systems, edited by K. Binder and G. Ciccotti
~Italian Physical Society, Bologna, 1996!, p. 255.
@5# D. Frenkel and A. J. C. Ladd, J. Chem. Phys. 81, 3188 ~1984!.
@6# D. Frenkel and B. Smit, Understanding Molecular Simulation
~Academic, New York, 1996!.
@7# N. B. Wilding and A. D. Bruce, Phys. Rev. Lett. 85, 5138
~2000!.
@8# In the case of hard spheres the temperature enters only trivially
through the scale of the pressure.
@9# D. A. Young, Phase Diagrams of the Elements ~University of
California Press, Berkeley, 1991!.
@10# T. H. K. Barron and M. L. Klein, in Dynamical Properties of036710the thermodynamically favored phase traverses a region with
a bcc structure. This phenomenon ~for which there is simu-
lation support @73#! is a useful reminder of the importance of
the kinetic pathways by which the equilibrium state is
reached, in the long term. The evidence from the most de-
tailed study known to us @74# is that a fcc structure emerges
in that ‘‘long term’’ @75#. But this study utilized a LJ poten-
tial truncated at rx52.3s; and we have seen that such details
matter as regards the ‘‘true’’ equilibrium behavior. In any
case, given the LJ knife edge, it seems likely that the par-
ticular long term behavior that does emerge may well have
more to do with the kinetics than the thermodynamics. That
interplay is probably better explored in a system other than
LJ. Recent work has begun to address this interesting area
@76#.
Solids, edited by G. K. Horton and A. A. Maradudin ~North-
Holland, Amsterdam, 1974!, p. 391.
@11# J. A. Barker, in Rare Gas Solids, edited by M. L. Klein and J.
A. Venables ~Academic Press, London, 1976!, Vol. 1, p. 212.
@12# One should not expect it to work at sufficiently extreme pres-
sures and it does not: see for example, R. J. Hemley and N. W.
Ashcroft, Phys. Today 51 ~8!, 26 ~1998!.
@13# Helium, of course, is different, see e.g., H. R. Glyde, Rare Gas
Solids ~Ref. @11#!, p 382.
@14# See e.g., K. F. Niebel and J. A. Venables, in Rare Gas Solids
~Ref. @11#!, p. 558.
@15# B. J. Alder and P. H. Paulson, J. Chem. Phys. 43, 4172 ~1965!.
@16# Compare Refs. @17,18#; see also Fig. 1.
@17# Y. Choi, T. Ree, and F. H. Ree, J. Chem. Phys. 95, 7548
~1991!.
@18# Y. Choi, T. Ree, and F. H. Ree, J. Chem. Phys. 99, 9917
~1993!.
@19# We use the term ‘‘lattice’’ rather loosely. The set of vectors
$
RW
%tion of the displacement coordinates, provided one takes ap-
propriate account ~through the Jacobian! of the resulting
warping of the configuration space @68#. In particular one can
define a switch operation that preserves collective ~normal
mode! coordinates rather than local coordinates. In the har-
monic region such an operation will work without the need
for any weighting; the barrier ~to be surmounted by multica-
nonical methods! is reduced to that associated with the an-
harmonic contributions. This approach is currently under in-
vestigation.
Second, the current framework is strictly classical. It is
clearly desirable to extend the framework to include quan-
tum effects; we plan to do so within a path integral formal-
ism @69#.
2. The system
We turn now to the particular system ~the LJ model! we3. The implications
If the calculations reported here are to represent more
than an exercise in equilibrium statistical mechanics then
their results need to be related to studies of the crystallization
process itself. It is not easy to do so in any satisfying way.
If one tries to make contact with observations on RGS one
is faced with a whole range of effects that are potentially
important. We have already noted the need to extend the
framework to include quantum effects—clearly important for
the lighter RGS and essential in understanding the behavior
of Helium. But the link with real systems is also complicated
by two other manifestations of the fcc-hcp knife edge on
which the LJ model sits. That knife edge is presumably re-
flected in the sensitivity of the observed structure to impuri-
ties @71#; it also means that effects that are missing from the
two-body central force LJ phenomenology ~and that are not
crucial to other aspects of the RGS behavior! need to be
addressed @14#. As a more limited objective one may try to
make contact with simulations of the crystallization process-10
Page 11
LATTICE-SWITCH MONTE CARLO METHOD: . . . PHYSICAL REVIEW E 65 036710identifies the orthodox lattice convoluted with the orthodox
basis. For the purposes of LSMC
$
RW
%
a
merely identifies one
particular configuration of structure a . The particular simplic-
ity ~perfect periodicity! of this configuration is not essential to
the LSMC procedure; the ‘‘lattice switch’’ is but a simple case
of a more general strategy that switches between configura-
tions of different phases @7#.
@20# We shall frequently suppress the a label when no ambiguity
results.
@21# Whenever we refer to the ‘‘difference’’ between the ~free! en-
ergies of the two structures we shall mean, specifically, the
value for fcc minus the value for hcp.
@22# The relevant lattice sums are defined in, e.g., @23#. The results
generated and used here are in agreement with those tabulated
in @24#.
@23# N. W. Ashcroft and N. D. Mermin, Solid State Physics ~Saun-
ders, Philadelphia, 1976!.
@24# T. H. K. Barron and C. Domb, Proc. R. Soc. London, Ser. A
227, 447 ~1955!.
@25# The ground-state energies were computed as parameterized
functions of density, in the limits of large nx and N; these
functions were then used in subsequent simulations.
@26# While this constraint deserves some thought @2# it is imple-
mented naturally and automatically within a MC framework.
@27# In the present work ~in contrast to our studies of hard spheres
@1,2,7#! we shall display factors of b5
@
kT
#
21 explicitly.
@28# See, for example, Ref. @24#; D. J. Lacks and G. C. Rutledge, J.
Chem. Phys. 101, 9961 ~1994!; R. G. Della Valle and E.
Venuti, Phys. Rev. B 58, 206 ~1998!.
@29# For the purposes of the harmonic theory the ‘‘lattice sites’’ are
assumed to be positions of classical equilibrium for the pre-
vailing mechanical constraints.
@30# See e.g., A. A. Maradudin, in Dynamical Properties of Solids,
edited by G. K. Horton and A. A. Maradudin ~North-Holland,
Amsterdam, 1974!, p. 18.
@31# More specifically f h is the harmonic contribution to the free-
energy density to within terms that do not depend upon the
phase label a and are thus irrelevant to the free-energy differ-
ence of interest.
@32# The origin of the ordinate scale in Fig. 2 is fixed by imple-
menting Eq. ~21! with the eigenvalues l j
a expressed in units of
e/s2.
@33# If one forms the abscissa for Fig. 2 from rx rather than nx each
f h function shows discontinuities at those rx values where nx
changes; the difference between the two functions shows dis-
continuities where either nx value change ~such as those seen
in the plot of the ground-state energy difference in Fig. 1!.
Parameterizing the cutoff by nx , and interpolating between the
values associated with the observable discrete values of this
variable, reveals essentially smooth underlying behavior.
@34# Z. W. Salsburg and D. A. Huckaby, J. Comput. Phys. 7, 489
~1971! evaluate the harmonic contribution to the free energies
of each phase, at the densities that minimize their respective
ground-state energies, using a potential including first- and
second-neighbor interactions.
@35# In simulations conducted in the NVT ensemble it is convenient
to factor out the difference between ground-state energies
DE0[NDe0 ~which is then independent of configuration! and
define M in terms of the excitation energies alone. In the NPT036710ensemble the ground-state energy difference is configuration
dependent ~varies with density! and there is less advantage to
making this separation.
@36# The sign convention adopted here has no deep significance.
@37# They lie near ~but not at! M50, see Sec. VII A.
@38# Here . . . identifies the appropriate ensemble variable V or P.
@39# One may see the use of multicanonical weights ~sampling with
them; then folding them out! as a way of ‘‘simply improving’’
~albeit infinitely! the statistical precision of an unbiased esti-
mator of the desired ratio.
@40# See, for example, the study of the freezing transition in @7#.
@41# There are other possibilities to which we return in Sec. VIII.
@42# In our studies of hard spheres @2# we treated the ratio of the
unit cell parameters c and a for the hcp structure as a simula-
tion variable. We were unable to resolve any departure from
the ideal close-packed limit c/a5A 83 . We did not reopen this
issue in the present studies: the c/a ratio is fixed throughout at
its ideal value.
@43# If the switch is required to accommodate a significant change
in some quantity other than the lattice type ~e.g., the volume!
the order parameter needs to be defined so as to carry an ex-
plicit phase label ~and it becomes less appropriate to refer to it
as an ‘‘order parameter’’!. The weight function h
@
M
#
then
splits into two distinct branches h
a
@
M
#
and the switch in-
volves a discontinuous change in the weight; this discontinuity
can be tuned to cancel the typical energy cost of the change in
the ‘‘other’’ quantity.
@44# By E(
$
uW
%
,
$
RW
%
a
) we mean the total configurational energy
E(
$
rW
%
) reparameterized ~for a configuration of phase a) by
$
uW
%
,a .
@45# G. R. Smith and A. D. Bruce, J. Phys. A 28, 6623 ~1995!; Phys.
Rev. E 53, 6530 ~1996!.
@46# The multicanonical weights are defined only to within an arbi-
trary additive constant; we fix that constant by the convention
that the macrostate of lowest probability has zero weight. The
weights may then be regarded as logarithmic measures of the
macrostate probabilities relative to that of the lowest-
probability macrostate. The height of the ‘‘peak’’ in the weight
function for a given phase then measures the probability bar-
rier to be surmounted, starting from the equilibrium states as-
sociated with that phase.
@47# Indeed for the LJ system @of the size, and under the conditions
defined in Fig. 5~a!# switches are occasionally possible even
without the aid of multicanonical weighting; in the correspond-
ing HS system such occurrences are some 1011 times less fre-
quent.
@48# Switches could also be launched successfully from states for
which DE˜ is negative but the scheme implemented here does
not need to locate such states and does not do so.
@49# See for example R. J. Hardy, J. Geophys. Res., @Space Phys.#
85, 7011 ~1980!.
@50# S. Somasi, B. Khomani, and R. Lovett, J. Chem. Phys. 113,
4320 ~2000!.
@51# The ~fractional! density difference at the phase boundary is at
most 231024.
@52# A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 61,
2635 ~1988!; 63, 1195 ~1989!.
@53# D. A. Kofke, J. Chem. Phys. 98, 4149 ~1993!.-11
basis. For the purposes of LSMC
$
RW
%
a
merely identifies one
particular configuration of structure a . The particular simplic-
ity ~perfect periodicity! of this configuration is not essential to
the LSMC procedure; the ‘‘lattice switch’’ is but a simple case
of a more general strategy that switches between configura-
tions of different phases @7#.
@20# We shall frequently suppress the a label when no ambiguity
results.
@21# Whenever we refer to the ‘‘difference’’ between the ~free! en-
ergies of the two structures we shall mean, specifically, the
value for fcc minus the value for hcp.
@22# The relevant lattice sums are defined in, e.g., @23#. The results
generated and used here are in agreement with those tabulated
in @24#.
@23# N. W. Ashcroft and N. D. Mermin, Solid State Physics ~Saun-
ders, Philadelphia, 1976!.
@24# T. H. K. Barron and C. Domb, Proc. R. Soc. London, Ser. A
227, 447 ~1955!.
@25# The ground-state energies were computed as parameterized
functions of density, in the limits of large nx and N; these
functions were then used in subsequent simulations.
@26# While this constraint deserves some thought @2# it is imple-
mented naturally and automatically within a MC framework.
@27# In the present work ~in contrast to our studies of hard spheres
@1,2,7#! we shall display factors of b5
@
kT
#
21 explicitly.
@28# See, for example, Ref. @24#; D. J. Lacks and G. C. Rutledge, J.
Chem. Phys. 101, 9961 ~1994!; R. G. Della Valle and E.
Venuti, Phys. Rev. B 58, 206 ~1998!.
@29# For the purposes of the harmonic theory the ‘‘lattice sites’’ are
assumed to be positions of classical equilibrium for the pre-
vailing mechanical constraints.
@30# See e.g., A. A. Maradudin, in Dynamical Properties of Solids,
edited by G. K. Horton and A. A. Maradudin ~North-Holland,
Amsterdam, 1974!, p. 18.
@31# More specifically f h is the harmonic contribution to the free-
energy density to within terms that do not depend upon the
phase label a and are thus irrelevant to the free-energy differ-
ence of interest.
@32# The origin of the ordinate scale in Fig. 2 is fixed by imple-
menting Eq. ~21! with the eigenvalues l j
a expressed in units of
e/s2.
@33# If one forms the abscissa for Fig. 2 from rx rather than nx each
f h function shows discontinuities at those rx values where nx
changes; the difference between the two functions shows dis-
continuities where either nx value change ~such as those seen
in the plot of the ground-state energy difference in Fig. 1!.
Parameterizing the cutoff by nx , and interpolating between the
values associated with the observable discrete values of this
variable, reveals essentially smooth underlying behavior.
@34# Z. W. Salsburg and D. A. Huckaby, J. Comput. Phys. 7, 489
~1971! evaluate the harmonic contribution to the free energies
of each phase, at the densities that minimize their respective
ground-state energies, using a potential including first- and
second-neighbor interactions.
@35# In simulations conducted in the NVT ensemble it is convenient
to factor out the difference between ground-state energies
DE0[NDe0 ~which is then independent of configuration! and
define M in terms of the excitation energies alone. In the NPT036710ensemble the ground-state energy difference is configuration
dependent ~varies with density! and there is less advantage to
making this separation.
@36# The sign convention adopted here has no deep significance.
@37# They lie near ~but not at! M50, see Sec. VII A.
@38# Here . . . identifies the appropriate ensemble variable V or P.
@39# One may see the use of multicanonical weights ~sampling with
them; then folding them out! as a way of ‘‘simply improving’’
~albeit infinitely! the statistical precision of an unbiased esti-
mator of the desired ratio.
@40# See, for example, the study of the freezing transition in @7#.
@41# There are other possibilities to which we return in Sec. VIII.
@42# In our studies of hard spheres @2# we treated the ratio of the
unit cell parameters c and a for the hcp structure as a simula-
tion variable. We were unable to resolve any departure from
the ideal close-packed limit c/a5A 83 . We did not reopen this
issue in the present studies: the c/a ratio is fixed throughout at
its ideal value.
@43# If the switch is required to accommodate a significant change
in some quantity other than the lattice type ~e.g., the volume!
the order parameter needs to be defined so as to carry an ex-
plicit phase label ~and it becomes less appropriate to refer to it
as an ‘‘order parameter’’!. The weight function h
@
M
#
then
splits into two distinct branches h
a
@
M
#
and the switch in-
volves a discontinuous change in the weight; this discontinuity
can be tuned to cancel the typical energy cost of the change in
the ‘‘other’’ quantity.
@44# By E(
$
uW
%
,
$
RW
%
a
) we mean the total configurational energy
E(
$
rW
%
) reparameterized ~for a configuration of phase a) by
$
uW
%
,a .
@45# G. R. Smith and A. D. Bruce, J. Phys. A 28, 6623 ~1995!; Phys.
Rev. E 53, 6530 ~1996!.
@46# The multicanonical weights are defined only to within an arbi-
trary additive constant; we fix that constant by the convention
that the macrostate of lowest probability has zero weight. The
weights may then be regarded as logarithmic measures of the
macrostate probabilities relative to that of the lowest-
probability macrostate. The height of the ‘‘peak’’ in the weight
function for a given phase then measures the probability bar-
rier to be surmounted, starting from the equilibrium states as-
sociated with that phase.
@47# Indeed for the LJ system @of the size, and under the conditions
defined in Fig. 5~a!# switches are occasionally possible even
without the aid of multicanonical weighting; in the correspond-
ing HS system such occurrences are some 1011 times less fre-
quent.
@48# Switches could also be launched successfully from states for
which DE˜ is negative but the scheme implemented here does
not need to locate such states and does not do so.
@49# See for example R. J. Hardy, J. Geophys. Res., @Space Phys.#
85, 7011 ~1980!.
@50# S. Somasi, B. Khomani, and R. Lovett, J. Chem. Phys. 113,
4320 ~2000!.
@51# The ~fractional! density difference at the phase boundary is at
most 231024.
@52# A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 61,
2635 ~1988!; 63, 1195 ~1989!.
@53# D. A. Kofke, J. Chem. Phys. 98, 4149 ~1993!.-11
Page 12
@54# An adaptation of an argument of Ref. @49# suggests that the
phase boundary temperature at density r will differ from the
harmonic prediction by an amount DT(r) that varies as
@
f
-
#
2/
@
f
9
#
3 where the derivatives ~of the LJ potential f) are
evaluated at the nearest-neighbor separation for the chosen
density. This is the source of the dashed-dotted line in Fig. 11.
@55# G. Jackson and F. van Swol, Mol. Phys. 65, 161 ~1988!.
@56# K. Singer and W. Smith, Chem. Phys. Lett. 140, 406 ~1987!.
@57# A. R. Philpott and J. M. Rickman, J. Chem. Phys. 94, 1454
~1991!.
@58# J. M. Rickman and A. R. Philpott, J. Chem. Phys. 95, 7562
~1991!.
@59# S. Y. Sheu, C. Y. Mou, and R. Lovett, Phys. Rev. E 51, R3795
~1995!.
@60# As an ~admittedly extreme! example, the difference between
the LJ fcc and hcp free energies at ~say! kT/e50.5 and p50 is
some four orders of magnitude smaller than the absolute ~ex-
cess! free energies of the two phases.
@61# S. Pronk and D. Frenkel, J. Chem. Phys. 110, 4589 ~1999! used
the integration methods of Ref. @5# to determine the ‘‘absolute’’
entropies of hcp and fcc phases of hard spheres and thence the
difference between these entropies, which is smaller by some
four orders of magnitude. They also performed studies using
the LS method: they report comparable accuracy for compa-
rable computational resource.
@62# A. Rahman and G. Jacucci, Nuovo Cimento D 4, 357 ~1984!.
@63# M. C. Moody, J. R. Ray, and A. Rahman, J. Chem. Phys. 84,
by a gradual shear ~mutual translation! of lattice planes.
@67# The similarities between LS and the work of Ref. @63#—not so
apparent in the LS formulation appropriate for hard spheres—
become more obvious when LS is reformulated to deal with
‘‘soft potentials.’’
@68# C. Jarzynski ~unpublished!, has provided a general formulation
of this idea.
@69# D. M. Ceperley and E. L. Pollock, in Monte Carlo Methods in
Theoretical Physics, edited by S. Caracciolo and A. Fabrocini
~ETS, Pisa, 1991!.
@70# By switching between ensembles with truncated and untrun-
cated potentials we investigated the effects of the truncation of
the potential @Eq. ~7!# on anharmonic effects. As one might
expect these effects are similar in scale to that apparent in the
harmonic results summarized in Ref. @2#. We found that the
consequences for the phase diagram were small compared to
the statistical uncertainties shown in Fig. 11.
@71# For example, Argon containing small concentrations of oxygen
crystallizes as hcp: C. S. Barrett and L. Meyer, J. Chem. Phys.
42, 107 ~1965!; see also T. Bricheno and J. A. Venables, J.
Phys. C 9, 4095 ~1975!. Note that ~in contrast to the implica-
tions of Fig. 11! these results suggest that the stability of hcp
relative to fcc is enhanced by raising the temperature.
@72# S. Alexander and J. P. McTague, Phys. Rev. Lett. 41, 702
~1978!.
@73# P. R. tenWolde, M. J. Ruiz-Montero, and D. Frenkel, Phys.
Rev. Lett. 75, 2714 ~1995!.
A. N. JACKSON, A. D. BRUCE, AND G. J. ACKLAND PHYSICAL REVIEW E 65 0367101795 ~1985!.
@64# B. Kuchta and R. D. Etters, Phys. Rev. B 47, 14 691 ~1993!.
@65# S.-C. Mau and D. A. Huse, Phys. Rev. E 59, 4396 ~1999!.
@66# We refer to the ‘‘shear implementation’’ of Ref. @65#. In this
case the ‘‘physical’’ path comprises configurations generated036710@74# W. C. Swope and H. C. Andersen, Phys. Rev. B 41, 7042
~1990!.
@75# Following quenches to kT/e.0.45, which is close to the fcc-
hcp phase boundary suggested in Fig. 11.
@76# S. Auer and D. Frenkel, Nature ~London! 409, 1020 ~2001!.-12
phase boundary temperature at density r will differ from the
harmonic prediction by an amount DT(r) that varies as
@
f
-
#
2/
@
f
9
#
3 where the derivatives ~of the LJ potential f) are
evaluated at the nearest-neighbor separation for the chosen
density. This is the source of the dashed-dotted line in Fig. 11.
@55# G. Jackson and F. van Swol, Mol. Phys. 65, 161 ~1988!.
@56# K. Singer and W. Smith, Chem. Phys. Lett. 140, 406 ~1987!.
@57# A. R. Philpott and J. M. Rickman, J. Chem. Phys. 94, 1454
~1991!.
@58# J. M. Rickman and A. R. Philpott, J. Chem. Phys. 95, 7562
~1991!.
@59# S. Y. Sheu, C. Y. Mou, and R. Lovett, Phys. Rev. E 51, R3795
~1995!.
@60# As an ~admittedly extreme! example, the difference between
the LJ fcc and hcp free energies at ~say! kT/e50.5 and p50 is
some four orders of magnitude smaller than the absolute ~ex-
cess! free energies of the two phases.
@61# S. Pronk and D. Frenkel, J. Chem. Phys. 110, 4589 ~1999! used
the integration methods of Ref. @5# to determine the ‘‘absolute’’
entropies of hcp and fcc phases of hard spheres and thence the
difference between these entropies, which is smaller by some
four orders of magnitude. They also performed studies using
the LS method: they report comparable accuracy for compa-
rable computational resource.
@62# A. Rahman and G. Jacucci, Nuovo Cimento D 4, 357 ~1984!.
@63# M. C. Moody, J. R. Ray, and A. Rahman, J. Chem. Phys. 84,
by a gradual shear ~mutual translation! of lattice planes.
@67# The similarities between LS and the work of Ref. @63#—not so
apparent in the LS formulation appropriate for hard spheres—
become more obvious when LS is reformulated to deal with
‘‘soft potentials.’’
@68# C. Jarzynski ~unpublished!, has provided a general formulation
of this idea.
@69# D. M. Ceperley and E. L. Pollock, in Monte Carlo Methods in
Theoretical Physics, edited by S. Caracciolo and A. Fabrocini
~ETS, Pisa, 1991!.
@70# By switching between ensembles with truncated and untrun-
cated potentials we investigated the effects of the truncation of
the potential @Eq. ~7!# on anharmonic effects. As one might
expect these effects are similar in scale to that apparent in the
harmonic results summarized in Ref. @2#. We found that the
consequences for the phase diagram were small compared to
the statistical uncertainties shown in Fig. 11.
@71# For example, Argon containing small concentrations of oxygen
crystallizes as hcp: C. S. Barrett and L. Meyer, J. Chem. Phys.
42, 107 ~1965!; see also T. Bricheno and J. A. Venables, J.
Phys. C 9, 4095 ~1975!. Note that ~in contrast to the implica-
tions of Fig. 11! these results suggest that the stability of hcp
relative to fcc is enhanced by raising the temperature.
@72# S. Alexander and J. P. McTague, Phys. Rev. Lett. 41, 702
~1978!.
@73# P. R. tenWolde, M. J. Ruiz-Montero, and D. Frenkel, Phys.
Rev. Lett. 75, 2714 ~1995!.
A. N. JACKSON, A. D. BRUCE, AND G. J. ACKLAND PHYSICAL REVIEW E 65 0367101795 ~1985!.
@64# B. Kuchta and R. D. Etters, Phys. Rev. B 47, 14 691 ~1993!.
@65# S.-C. Mau and D. A. Huse, Phys. Rev. E 59, 4396 ~1999!.
@66# We refer to the ‘‘shear implementation’’ of Ref. @65#. In this
case the ‘‘physical’’ path comprises configurations generated036710@74# W. C. Swope and H. C. Andersen, Phys. Rev. B 41, 7042
~1990!.
@75# Following quenches to kT/e.0.45, which is close to the fcc-
hcp phase boundary suggested in Fig. 11.
@76# S. Auer and D. Frenkel, Nature ~London! 409, 1020 ~2001!.-12
Sign up today - FREE
Mendeley saves you time finding and organizing research. Learn more
- All your research in one place
- Add and import papers easily
- Access it anywhere, anytime
Start using Mendeley in seconds!
Readership Statistics
6 Readers on Mendeley
by Discipline
17% Chemistry
by Academic Status
33% Post Doc
17% Lecturer
17% Researcher (at an Academic Institution)
by Country
33% United Kingdom
17% Sweden
17% China
Groups
 |