n. J
Ed
A
al
-sw
y
at
s
ep
re
et
eth
ow
bi
e
PHYSICAL REVIEW E JANUARY 2000VOLUME 61, NUMBER 1this problem @1#. But it is clear that if one desires a technique
that is both generally applicable and reliable ~that is, has
quantifiable uncertainties! one must look to the Monte Carlo
~MC! method @2#, the standard computational tool for deal-
ing with many-body systems @3#.
The application of MC methods to the study of phase-
behavior presents a generic problem @4,5#: the free energy of
a phase cannot be expressed ~in a practically useful form! as
a canonical average over the associated configurations; free-
energy-estimation inevitably entails simulations that visit a
substantially wider spectrum of configurations, which to-
gether form a path through configuration space @6#. The stra-
tegic choices to be made concern the path itself—ultimately,
the physical character of the additional configurations
sampled—and the sampling procedure.
An acceptable path will fall into one or other of two
categories—we will call them reference-state and interphase
paths. A reference-state path links ~comprises sets of con-
figurations that interpolate between! the configuration space
associated with each phase to the configuration space asso-
ciated with some reference system @7# whose free energy is
known. An interphase path links the configuration space of
each point @the integration method ~IM!# or the difference
between the free energies of adjacent points ~the overlap
method!. The single-stage approach involves, in essence one
simulation exploring the entire path.
There are very many ways in which one can respond to
these strategic choices. Many of them are represented in the
large literature devoted to this problem @8–14#. But all of
them, in our view, lack one or more of the characteristics
~generality, transparency, precision! of a definitively satisfac-
tory solution to such a fundamental and simply posed prob-
lem.
In seeking that solution it seems to us there are good a
priori grounds for favoring an inter phase path, explored by
single-stage sampling. The prejudice on the choice of path
reflects the fact that, in using a reference state path, one has
to determine, separately, the absolute free energies of each
phase. These absolute free energies are typically very
large—arbitrarily so in the vicinity of a phase boundary—
compared to the quantity ~their difference! which is actually
of interest. In contrast, using an interphase path allows one to
focus directly on this quantity. The a priori preference for a
single-stage sampling rests on the transparency with whichLattice-switch Mo
A. D. Bruce, A. N. Jackson, G
Department of Physics and Astronomy, The University of
~Received 20
We present a Monte Carlo method for the direct ev
two crystal structures. The method is built on a lattice
structure onto a candidate configuration of the other b
while keeping the displacements with respect to the l
configurations is biased, multicanonically, to favor path
Carlo switch to the candidate configuration will be acc
efficiently sampled in a single process, and the diffe
measured probabilities. We explore and exploit the m
hard spheres. We show that the efficiency of the m
conserves correlated microstructure. We also show h
finds gateway arrangements which fulfill the sampling
differences between the free energies of the two clos
density and the constant pressure ensembles.
PACS number~s!: 05.10.Ln, 65.50.1m, 64.70.Kb
I. INTRODUCTION
Let us pose the problem. We are presented with a material
whose chemical composition is known; we are provided with
a model of the interatomic interactions; and we have identi-
fied two candidate crystalline structures. How should we pro-
ceed to determine which structure will be favored under
given conditions? Of course, equilibrium statistical mechan-
ics tells us what we must do, in principle: the favored struc-
ture will be that which has the greater a priori probability or
configurational weight or, in equivalent thermodynamic par-
lance, lower free energy. Thus the task is to compare the
configurational weights of ~determine the difference between
the free energies of! the candidate structures.
A variety of approximate strategies exist for addressingPRE 611063-651X/2000/61~1!/906~14!/$15.00te Carlo method
. Ackland, and N. B. Wilding
inburgh, Edinburgh EH9 3JZ, Scotland, United Kingdom
ugust 1999!
uation of the difference between the free energies of
itch transformation that maps a configuration of one
‘‘switching’’ one set of lattice vectors for the other,
tice sites constant. The sampling of the displacement
leading to gateway arrangements for which the Monte
ted. The configurations of both structures can then be
nce between their free energies evaluated from their
hod in the context of extensive studies of systems of
od is controlled by the extent to which the switch
, microscopically, the procedure works: the system
as intelligently. We establish, with high precision, the
packed structures ~fcc and hcp! in both the constant
one phase to that of the other. Both categories of path em-
brace many further subcategories. Thus, a reference-state-
path may run through a space of thermodynamic coordinates
or through a space of model parameters. An interphase path
may be ‘‘physically motivated’’—modeling authentically the
configurations through which a system actually passes in the
course of a phase transformation or it may be ‘‘computation-
ally motivated’’ ~‘‘nonphysical’’!—designed, pragmatically,
to deliver a result.
The sampling procedures used to explore the chosen path
also fall, broadly, into one or other of two categories – we
will call them multistage and single stage. The multistage
approach entails a series of independent simulations each of
which explores a different point on the path; the simulations
may determine simply the derivative of the free energy at906 ©2000 The American Physical Society

Nhanced to provide a range of extended sampling
techniques—multicanonical @15#, expanded ensemble @16#,
and simulated tempering @17#. These methods ~whose origins
can be traced back to much earlier pioneering work @18#!
allow one to construct a MC procedure that will traverse
virtually any desired path through configuration space. Here
we adopt the multicanonical framework. In this framework,
the desired path is represented as a discrete series of mac-
rostates, defined by some chosen macroscopic property @6#;
in the multicanonical sampling procedure each macrostate is
visited with a probability that is enhanced, or diminished,
with respect to its canonical value, by an amount that is
controlled by a multicanonical weight; the set of weights is
constructed so that, while the canonical probabilities vary
vastly over the path, the multicanonical probabilities are es-
sentially constant, allowing the whole path to be negotiated
in one simulation.
The core issue is, then, the design of the interphase
path—at heart, the choice of an appropriate order parameter
@19#. The choice is important: it determines, implicitly @20#,
the nature of the configurations sampled during the inter-
phase traverse, the MC time required for that traverse, and
thence the statistical quality of the final results.
Outside the context of structural-phase behavior—in the
case of liquid-gas phase behavior, for example—the choice
is clear and a multicanonical strategy is securely in place.
ready represented in the microscopic dynamics of a single
phase. This approach is illustrated schematically in Fig. 1~a!.
It has been successfully used in studies of phase behavior in
ferromagnets @22#, fluids @23#, and lattice gauge theories
@24#.
In the context of structural phase behavior it is clear that
this kind of strategy will seldom be fruitful @25,26#. In such
systems a traverse through an inhomogeneous two-phase
~necessarily noncrystalline! region will involve substantial,
physically slow, restructuring—vulnerable to further ergodic
traps, and compounding the intrinsic slowness of the multi-
canonical sampling process. To the two general a priori pref-
erences expressed above we thus add a third, specific to the
structural context: the interphase path should comprise mac-
rostates that are single phase and crystalline. This paper
shows how to identify, build, and exploit a path of this kind.
The key ideas are simple. In any crystalline configuration
each atomic position coordinate may be expressed as the sum
of a lattice vector and a displacement vector. The configura-
tion space associated with each structure, individually, may
be explored by standard MC procedures which stochastically
update the displacement vectors while keeping the lattice
vectors constant. In principle the passage from one phase to
the other may be accomplished by a lattice switch ~LS! in
which one entire set of lattice vectors is replaced with the
other, while the displacement vectors are held fixed. For-the associated uncertainties ~error bounds! are prescribed.
We shall return to these points in Sec. V. With these strategic
choices made, one is left with two tasks—one conceptual
~designing the interphase path! and the other practical ~for-
mulating the sampling algorithm!.
The practical issue is relatively easily addressed. In recent
years, the Monte Carlo toolkit has been significantly en-
PRE 61 LATTICE-SWITCH MOThe order parameter is identified with that—the density—
associated with the accompanying critical point. The result-ing interphase configurations are then generically inhomoge-
neous, comprising two coexisting regions ~one of each
phase!, separated by an interface. On this path, it is the free
energy cost of this interface that provides the ergodic barrier
which has to be surmounted by multicanonical weighting
@21#. The passage along the path ~the motion of the interface!
involves processes which differ only in scale from those al-
FIG. 1. Schematic representations of the dif-
ferent ways in which multicanonical sampling
methods can be used to achieve interphase cross-
ing. In the conventional approach ~a! the sam-
pling algorithm is biased so as to enhance the
probability of the mixed phase states lying along
a path ~the heavy dark line! linking the two re-
gions of configuration space. In the lattice-switch
method ~b! the bias is constructed so as to en-
hance the probability of the subsets of states ~the
white islands!, within the single-phase regions,
from which the switch operation ~the large
dashed arrow! will be accepted.
907TE CARLO METHODmally this LS can be incorporated into the MC procedure
simply by treating the lattice type as an additional stochastic

KLvariable. In practice this interphase ‘‘path’’ ~blind leap! will
not work. Implemented this way the LS will map a ‘‘typi-
cal’’ configuration of one structure onto an ‘‘untypical’’
~high-energy! conjugate configuration @27#; the associated
MC step will generally be rejected. To make it work the LS
needs to be extended to include two segments of ‘‘path’’
~each lying entirely within one phase! which connect the sets
of equilibrium configurations with the special configurations
~we will call them gateway configurations! from which a
successful LS can be initiated @28#. These path-segments
may be labeled by an ‘‘order parameter’’ which measures the
mismatch between the energies of the configurations linked
by LS. This order parameter has a high value for the equi-
librium configurations, lying at one end of a path segment:
these configurations are not energy matched @29# to their
conjugates. It has a low value for the gateway configurations
at the other end: gateway configurations ~whatever other at-
tributes they may have! are necessarily energy matched to
their conjugates. Multicanonical weights are attached to the
macrostates of this order parameter, so that the multicanoni-
cal sampling procedure explores both path segments evenly,
surmounting the probabilistic barrier which, in this case, re-
flects the smallness of the statistical weight of the gateway
configurations. Together, the multicanonical sampling and
the lattice switch provide a configuration space ‘‘look and
leap’’ @Fig. 1~b!# which visits both phases while remaining at
all times crystalline.
The LS method was introduced by us and described in
outline form in an earlier brief communication @30#. Since
that time it has been applied by two other groups @14,31#.
The present paper has three principal objectives.
The first objective ~with which we have already engaged
in the preceding discussion! is to explain the core idea more
fully: the ‘‘idea’’ ~biased sampling to facilitate a global co-
ordinate change! represents, we believe, a significantly new
form of extended sampling, which merits further exposure.
Our second objective is to achieve a deeper understanding
of how the process works—in particular the implications of
the form chosen for the LS operation adopted ~it is not
unique! and the microscopic character of the gateway con-
figurations which the system locates in response to the mul-
ticanonical weighting, tailored to support that operation. We
show that the efficiency of the LS operation depends signifi-
cantly on the extent to which it conserves correlated micro-
structure. And we find that the gateway configurations have
features which reflect the specific nature of the lattice-switch
transformation we adopt, in a microscopically intelligible
~even intelligent! way.
Our third objective is to extend our study of the phase
behavior of hard spheres. This problem is of enduring inter-
est, displaying a richness that belies the simplicity of the
model itself @32#. The relative stability of the two closed-
packed ~fcc and hcp! structures is particularly finely bal-
anced: the entropy difference @33# is so small ~smaller than
the entropy change at freezing by of order 1023) that it can
easily be lost in statistical uncertainties. Discrepancies ~large,
in relative terms! between a recent IM study @10# and its
predecessors @9# provided the motivation for our develop-
ment of the LS method. In this paper we present results in
908 BRUCE, JACKSON, ACthe constant-density ensemble, both near the melting density
and at the close-packed limit. In so doing we resolve thediscrepancy between the results, near melting, reported in
our initial study @30# and those—also using LS—reported
recently by Pronk and Frenkel @31#: the fault was ours, stem-
ming from a failure to recognize the consequences of center-
of-mass drift. We also show that the method can be extended
straightforwardly—in this case at least—to the constant-
pressure ensemble.
The paper is structured as follows. Section II sets out the
theoretical framework. We define the model, the competing
structures, and the associated configurational weights: in the
case of hard sphere systems the latter are purely entropic. We
identify an appropriate form of lattice-switch transformation:
here, it is designed to capitalize on the close-packed layers
common to both structures. To bias the displacement sam-
pling we need to define an appropriate measure of the ‘‘en-
ergy cost’’ of the lattice switch; we will see that the number
of pairs of overlapping spheres created by the transformation
fulfills this role simply and effectively. The efficiency of the
method also potentially depends on the choice of representa-
tion of both the lattice-to-lattice mapping and the particle
displacements: we discuss the principles involved in the
choice of representation. Section III provides computational
implementation details, including the procedures used to
evolve an appropriate multicanonical sampling distribution.
Section IV contains our results. Finally, in Sec. V, we offer
our conclusions in relation to both the hard sphere system
and the lattice-switch method.
II. FORMULATION
A. The model system
We consider a system of N particles, of spatial coordi-
nates
$
rW
%
, confined within a volume V, and subject to peri-
odic boundary conditions. The interactions are those of hard
spheres of diameter D; the configurational energy is of the
form
E
~
$
rW
%
!5
H
0 if ri j>D ;i , j ,
‘ otherwise, ~1!
where ri j5urW i2rW ju. The total configurational weight associ-
ated with this system is
V
~
N ,V !5
)
i
F
E
V
drW i
G
)
^
i j
&
Q
~
ri j2D !, ~2!
where Q(x)[1(0) for x>0(,0), and the product on
^
i j
&
extends over all particle pairs. The associated entropy den-
sity is
s
~
N ,V ![
1
N ln V~N ,V !. ~3!
We are concerned with the entropy of specific phases ~the
two familiar crystalline close-packed structures! of this sys-
tem. In general, the entropy of a phase measures the weight
of the configurations satisfying some constraint that is char-
acteristic of that phase. It is therefore necessary in principle
~although in practice the issue is typically skirted! to formu-
PRE 61AND, AND WILDINGlate a constraint that identifies a configuration as ‘‘belonging
to’’ a given crystalline phase. One can do so—very naturally,

Nand in the traditions of lattice dynamics @34#—by decompos-
ing the particle position coordinates into a sum of ‘‘lattice’’
and ‘‘displacement’’ vectors:
rW i5RW i
a
1uW i . ~4!
Here
$
RW
%
a
[RW i
a
,i51, . . . ,N is a set of fixed ~configuration-
independent! vectors associated with the crystalline structure
labeled a . We will refer to them as ‘‘lattice vectors.’’ But we
use this term a little loosely: more precisely, we mean the set
of vectors identified by the orthodox crystallographic lattice,
convolved with the orthodox basis @35,36#. The other vectors
$
uW
%
represent displacements with respect to the ‘‘lattice’’
sites; the symmetry of the structures of interest here ensures
that these displacements have zero ensemble average. This
framework provides us with a number of ways of identifying
the configurations to be associated with structure a . First,
one might adopt the criterion that all particle displacements,
with respect to the associated lattice sites, lie within some
nominated spatial cutoff, chosen to be sufficiently large that
the results are independent of its specific value. This crite-
rion has the merit that it does not stray beyond the concepts
of equilibrium statistical mechanics. Alternatively one might
identify the relevant configurations as the set that are acces-
sible from some member of the set ~the perfect crystalline
state, for example! within some nominated temporal cutoff.
The merit of this choice is that it is a quasiformal expression
of what, in practice, computer simulation attacks on this
problem actually do, albeit implicitly: the free energy as-
signed to a phase ~in, for example, IM-based studies! repre-
sents the weight of the configuration space sampled on the
time scale of the simulation. The result should be indepen-
dent of that time scale provided it ~the scale! is long enough
that the configuration space of each structure is effectively
sampled, but still short compared to interphase crossing
times. Whichever view one takes ~in practice we adopt the
latter: see Sec. III A! one may write, for the configurational
weight associated with structure a
V
~
N ,V ,a!5
)
i
F
E
a
duW i
G
)
^
i j
&
Q
~
ri j2D !, ~5!
where *
a
signifies integration subject to the chosen configu-
rational constraint.
In the thermodynamic (N→‘) limit, the associated en-
tropy density
s
~
N ,V ,a![
1
N ln V~N ,V ,a! ~6!
depends only on the particle number density, which we write
in the dimensionless form
r
˜
[
r
rCP
[
N/V
A2/D3
, ~7!
where rCP , the number density at close packing, provides
the natural scale. The range of interest to us here extends
PRE 61 LATTICE-SWITCH MOfrom the melting density r˜.0.736 @37# through to the close-
packed limit r˜51.In the close-packed limit the configurational integral @Eq.
~5!# may be rewritten @38# as the product of two terms:
V
~
N ,V ,a!5V0~N ,V !Va ~8!
The first term here is defined by
V0~N ,V !5F
De
12eG
3N
~9a!
with
e[12r˜ 1/3. ~9b!
The associated contribution to the entropy is logarithmically
divergent in the close-packed limit @39#, but independent of
the phase. The second contribution to the configurational in-
tegral is defined by
V
a
5
)
i
F
E
a
duW i
G
)
^
i j
&
nn
Q
~
ui j
uu
11 !@11O~e!# , ~10a!
where @40#
uW i j[uW i2uW j[ui j
uu nˆ i j
a
1uW i j
’
~10b!
while nˆ i j
a is a unit vector from lattice site j to nearest neigh-
bor lattice site i. The associated contribution to the entropy is
finite, but depends on the phase through the geometry of the
nearest neighbor vectors. It may be visualized as that of a set
of hard dodecahedra @41#.
Now let us recall that the quantity of immediate interest is
the difference between the entropy densities of the two
phases. It may be written as
Ds
ab
[s
~
N ,V ,a!2s~N ,V ,b!5
1
N ln Rab~N ,V !, ~11!
where
R
ab
~
N ,V ![
V
~
N ,V ,a!
V
~
N ,V ,b! 5
P
~
auN ,V !
P
~
buN ,V ! . ~12!
Here P(auN ,V) is the probability that a system, free to ex-
plore the joint configuration space of the two structures ~and
visiting configurations with the appropriate probabilities—all
equal in this case! will be found in some configuration of
structure a .
In the constant density ensemble, then, the computational
task is to determine the ratio defined by Eq. ~12!. In the
constant pressure ensemble we require the ratio R
ab
(N ,P!)
of the partition functions
Z
~
N ,P!,a!5
E
dVV
~
N ,V ,a!e2P
!V
, ~13!
where P! is a measure of the pressure @42#. The associated
thermodynamic potential is the Gibbs free energy density
defined by
909TE CARLO METHODg
~
N ,P*,a![2
1
N ln Z~N ,P
!
,a! ~14!

KLso that, in analogy with Eq. ~11!,
Dg
ab
[g
~
N ,P!,a!2g~N ,P!,b!52
1
N ln Rab~N ,P
!
!.
~15!
B. The lattice-switch method
The two close-packed structures of interest here are
shown schematically in Fig. 2. In principle there are many
transformations which will map one set of lattice vectors into
the other; we shall consider the criteria guiding the choice in
Sec. II C. The mapping used in most of the work reported
here is shown schematically in Fig. 3. This scheme exploits
the fact that the two structures differ only in respect of the
stacking pattern of the close-packed planes. A suitable trans-
formation can then be constructed that entails, simply, trans-
lating appropriate close-packed planes. By ‘‘translate’’ we
mean, more precisely, ‘‘relocate at a position defined by an
FIG. 2. Schematic representations of the two close-packed
structures. The structures differ only in regard to the stacking
pattern of the close-packed (x-y) planes which are of the form
ABCABC••• for fcc ~upper! and ABAB••• for hcp ~lower!. The
vector labeled tW is instrumental in defining the LS operation, shown
in Fig. 3.
FIG. 3. The LS transformation applied to the perfect-crystal
configuration. The diagram shows 6 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 hard
spheres located at the sites of the two close-packed structures. In
this realization of the fcc→hcp lattice switch, the top pair of planes
are left unaltered, while the other pairs of planes are relocated by
910 BRUCE, JACKSON, ACtranslations, specified by the vectors 2 tW ~white arrows! and tW
~black arrows!. The vector tW is identified in Fig. 2.appropriate translation vector’’: one should not think of the
planes as ‘‘sliding through’’ the intermediate states.
Figure 3 shows the application only to the perfect-crystal
configurations where energy-matching is guaranteed @43#. In
general ~that is, for ‘‘typical’’ configurations: see Fig. 4 for
an example! the two configurations related by the LS opera-
tion will not be energy matched: since adjacent planes are
translated differently, the translations may—indeed, with
overwhelming probability will—map a realizable configura-
tion ~of one structure!, in which there are no overlapping
spheres, onto an unrealizable configuration ~of the other! in
which there are overlaps. A MC lattice switch ‘‘move’’ will
be rejected for most configurations. But not quite all: gate-
way configurations ~configurations that are energy matched
@29# to their conjugates! must exist, in significant measure. In
particular, it is clear on grounds of continuity that configu-
rations ‘‘close enough’’ to perfect-crystal form must fall into
this category. One might therefore choose these ‘‘small-
displacement’’ configurations to act as the gateway states,
and define a multicanonical weighting procedure accord-
ingly. However, one can avoid having to make this explicit
choice, and, instead, let the system find gateway configura-
tions itself. To do so we must define a measure of the mis-
match between the energies of the configurations linked by
the transformation.
In the present context that mismatch is quantified by the
number of pairs of overlapping spheres created by the trans-
formation. To that end let M (
$
uW
%
,a) denote the number of
overlapping pairs associated with the displacements
$
uW
%
within the structure a . And define @44#
M
~
$
uW
%
![M ~
$
uW
%
,hcp!2M ~
$
uW
%
,fcc!. ~16!
Since M (
$
uW
%
,a) will necessarily be zero for any realizable
set of displacements of structure a , the overlap order param-
eter M is necessarily >0(<0) for realizable configurations
of the fcc ~hcp! structure. Figure 4 provides a concrete ex-
ample. The gateway configurations may then be identified
FIG. 4. The LS transformation applied to a ‘‘typical’’ configu-
ration. The crosses identify the ‘‘lattice sites’’; the small circles
locate the sphere centers in this configuration of displacements
$
uW
%
.
This configuration is realizable ~gives no overlaps! in the fcc struc-
ture; under the LS transformation it is mapped onto an ~unrealiz-
able! hcp configuration with three overlapping pairs of hard spheres
~shown with dashed boundaries!. Thus, for this configuration, the
overlap order parameter M(
$
uW
%
)53 @Eq. ~16!#.
PRE 61AND, AND WILDING~without prejudging their microscopic character! as the set of
configurations for which M50: a displacement pattern
$
uW
%

Nfor which M50 is realizable in both structures ~energy
matched!. A LS MC step initiated from an M50 configu-
ration will be accepted; if initiated from outside this set of
configurations it will be rejected.
The sampling algorithm must thus be multicanonically
customized so as to enhance the probability along a notional
line in M space, extending from the ‘‘equilibrium’’ M val-
ues ~reflecting the number of overlaps created by a LS acting
on a typical configuration! through to the M50 gateway
configurations. This aim is realized by augmenting the sys-
tem energy function Eq. ~1!:
E
~
$
rW
%
!→E~
$
rW
%
!1h„M~
$
uW
%
!…[E˜ ~
$
rW
%
! ~17!
where h(M),M50,61,62••• constitute a set of multica-
nonical weights @15#. These weights need to be chosen so as
to allow the system to access the M50 gateway configura-
tions, and thence ~through the LS! the full joint configuration
space of the two structures. The desired ratio of configura-
tional weights, which reflects the canonical distribution
P(MuN ,V) @Eq. ~12!# may then be estimated from the mea-
sured multicanonical distribution, P(MuN ,V ,
$
h(M)
%
) with
the identification
Rfcc,hcp~N ,V !5
(M.0
P
~
MuN ,V !
(M,0
P
~
MuN ,V !
5
(M.0
P
~
MuN ,V ,
$
h
~
M!
%
!eh(M)
(M,0
P
~
MuN ,V ,
$
h
~
M!
%
!eh(M)
.
~18!
Here the exponential reweighting of the multicanonical dis-
tribution folds out the bias associated with the weights,
whose residual effects are then simply as desired—the re-
moval of the ergodic barrier between the two branches of the
distribution.
C. Representations: tuning the lattice switch
We have presented the LS method in its simplest
realization—the one we have used for most of the studies
reported here. We now outline two important respects @45# in
which some degree of generalization is possible, and may be
desirable, in subsequent applications. Both involve the
choice of representation of the LS transformation.
We have already alluded to the first point: there are many
forms of lattice-to-lattice mapping. It is clear that the effi-
ciency of the method will depend significantly upon the map-
ping chosen. Evidently the choice should be made so as to
match up, as closely as possible, the energy of the two con-
figurations it links. In the context of hard spheres this aim is
realized by choosing the mapping which gives the smallest
equilibrium overlap count ~mean uMu value!, which gives a
measure of the entropic barrier that has to be negotiated by
PRE 61 LATTICE-SWITCH MOthe multicanonical procedure. The smaller this barrier, the
shorter is the path to the gateway configurations from whicha successful LS may be launched. Since the multicanonical
simulations traverse this path only slowly ~essentially diffu-
sively, at best! the gains here are potentially substantial. It is
intuitively clear that the scheme described above will fulfill
this criterion well: in this representation, the LS translates
close-packed planes bodily, so it can create overlaps only
between spheres associated with different planes. But it is
useful to explore other schemes—partly to check that there is
no significantly better alternative, but principally to under-
stand the different factors that control the efficiency. We
have done so; the results are to be found in Sec. IV A.
There is a second—less obvious—generalization of the
framework. In the simple realization, the particle positions
are written in the ‘‘lattice plus displacement’’ representation
provided by Eq. ~4!. The LS operation then maps a configu-
ration of one structure onto a configuration of the other with
the same set of displacements. This is unnecessarily restric-
tive. More generally we are at liberty to write, in place of Eq.
~4!,
rW5RW a1Ta•uW , ~19!
where rW ,RW a, and uW are now column vectors with 3N ele-
ments and Ta is a 3N33N non-singular matrix, whose form
~possibly
$
uW
%
-dependent! is at our disposal. Equation ~5! is
then replaced by
V
~
N ,V ,a!5
)
i
F
E
a
duW i
G
•det Ta
)
^
i j
&
Q
~
ri j2D !. ~20!
From the standpoint of the ~standard! single-phase part of the
MC procedure, this change in representation is equivalent to
changing the form of the configurational energy:
E
~
$
rW
%
!→E~
$
rW
%
!2ln@det Ta# , ~21!
This change introduces some computational overheads,
which could be substantial if the T transformation is not
local. The potential pay off lies in the LS part of the MC
procedure. One might hope to be able to tune the form of the
T matrix so that ‘‘typical’’ configurations of the one struc-
ture are mapped ~by LS! into ‘‘typical’’ configurations of the
other. In the case of the hard sphere problem, however, our
results ~Sec. IV A! suggest that there is little to be gained
here by this kind of tuning.
III. IMPLEMENTATION
A. Monte Carlo procedures
First we consider the procedure for MC sampling of the
particle displacements, for a given structure ~set of lattice
vectors!. As discussed in Sec. II A this sampling should, in
principle, satisfy some appropriate configurational constraint
@46#. In our original studies @30# we chose to implement this
constraint explicitly, through our sampling distribution: can-
didate displacements were drawn from a flat ~‘‘top hat’’!
distribution. This procedure can be made to work. But the
constraint explicitly breaks the translational invariance; and
one must deal with the consequences. In particular the con-
911TE CARLO METHODfigurational integral effectively being evaluated then depends
upon the location of the center of mass and thence upon the

KLtop-hat cutoff; this dependence sets in when the displace-
ment acquired by the center of mass, in the course of its slow
diffusive motion, becomes comparable with the top-hat cut-
off. One can avoid this problem simply by fixing the center
of mass. Our failure to do so in Ref. @30# led to results which
differ significantly from those we present here. In the studies
reported here we have chosen the ‘‘implicit’’ realization of
the configurational constraint ~practically, but not conceptu-
ally equivalent to ignoring it! which rests ~Sec. II A! on time
scales. Spheres were chosen at random, and trial changes to
the current displacement drawn from a uniform distribution.
The displacement update is accepted according to the Me-
tropolis prescription @3#
pa~$uW %→$uW 8%!5min$1,exp@2DE˜ ~$rW%!#% ~22!
where E˜ (
$
rW
%
) is defined in Eq. ~17! @47#. In addition to the
constraint that the update should yield a realizable configu-
ration of the current phase, this acceptance probability re-
flects the chosen weights which are defined ~Sec. III B ex-
plains how! on the space of the overlap order parameter M
@Eq. ~16!#. To minimize the computational time spent deter-
mining how a proposed move affects the value of M we
used a local overlap array, holding information on which
neighbors of a given sphere currently overlap with that
sphere in the conjugate configuration generated by a LS.
The representation of the close-packed limit provided by
Eq. ~10a! can be handled with only minor amendments: the
constraint ri j.D identifying realizable configurations is re-
placed by a constraint on the scaled displacement-difference
coordinates ui j
uu
.21. The overlap order parameter ~measur-
ing the number of times the hard sphere constraint is violated
in the conjugate configuration! is redefined accordingly. In
this limit particle ‘‘interactions’’ ~encounters! may occur
only between immediate neighbors. At other densities we
allowed for the possibility of encounters between nominal
second neighbors. We found, however, that although the
number of such encounters grows rapidly with the approach
to the melting density, the consequences for the relative en-
tropy of the two structures is insignificant under the condi-
tions studied here @48#.
In addition to particle moves the constant-pressure simu-
lations require updates of the simulation-cell parameters. In
such an update ~implemented on average once per sweep! a
trial set of cell parameters are selected, and accepted with
probability @49#
pa~V→V8!5min$1,exp@2DE˜ ~$rW%!2P!DV
1N ln
~
V
8
/V !#
%
, ~23!
where V
8
is the volume associated with the trial parameters.
Note that this kind of update—a dilation—changes E˜ (
$
rW
%
)
both trivially ~so as to forbid moves causing ‘‘real’’ over-
laps! and more subtly through changes in the count of the
overlaps in the conjugate configuration. A volume update
thus requires recalculation of the entire local overlap array.
Now consider the lattice switch. The switch may be
viewed as an updating of the ‘‘lattice’’ type a , regarded as a
912 BRUCE, JACKSON, ACstochastic variable. The prescription for such an update is
quite simple. After every particle update the value of M ischecked ~it is already known!. The LS is performed if ~and
only if! the gateway condition M50 is satisfied.
B. Calculating the weights
The determination of an acceptable set of multicanonical
weights @15# can be accomplished in a number of ways—
none, seemingly, entirely systematic. We describe briefly the
techniques we have used in the present study. Figure 5 pro-
vides some illustration. For further details and references to
other work the reader is referred to Refs. @15,26,50,51#.
The simplest method is the visited-state ~VS! technique
@50#. In this approach a suitable set of weights is evolved
through an iterative process ~Fig. 5!, the next set of weights
depending upon the distribution of the ~overlap! order pa-
rameter over the macrostates visited using the current set of
weights. This process is repeated until the weights yield a
distribution P(MuN ,V ,
$
h(M)
%
) that is effectively flat. This
method proved quite adequate for our smallest system.
For larger systems, however, we found it more efficient to
appeal to the transition probability ~TP! method @50#. In the
simplest realization of this method the simulation is initiated
from a ‘‘cold’’ ~zero displacement! configuration ~a member
of the M50 macrostate! for one structure. In the course of
its subsequent evolution towards equilibrium for that struc-
ture the numbers of transitions between different M mac-
rostates are recorded, and subsequently used to construct an
estimator of the macrostate-transition-probability matrix.
This TP matrix can be used to estimate the macrostate prob-
ability distribution and thence to provide an estimate for a set
of weights ~Fig. 5!, which can in turn be refined via VS. For
our intermediate size system this method worked well.
FIG. 5. Illustration of the weight-generation process, for a sys-
tem of N5216 hard spheres. The points marked VS are the results
of the first 3 iterations of the visited-states algorithm, initiated from
an fcc equilibrium state. The points marked TP emerge from one
application of the transition probability method. The solid line
shows a refined ~usable! set of weights.
PRE 61AND, AND WILDINGIn the case of our largest system we found it necessary to
modify the method somewhat, so as to limit the rate at which

Nthe simulation passes through M space. One way of doing
this is to constrain the system to macrostates with overlap
order parameters below some ‘‘barrier’’ value, which is
gradually incremented ~moved ‘‘out’’ in M space!, at inter-
vals of the order of the equilibration time.
By fiat the two structures have the same weights for M
50. In principle, the weights associated with the two struc-
tures for nonzero uMu are different @i.e., h(M)Þh
(2M)], and have to be evolved separately. In practice the
weights of the two structures are very similar—a reflection
of the similarity of the entropies of the two phases. Conse-
quently, one set of weights provides an excellent first ap-
proximation to the other, for refinement by VS.
C. Simulation details
The specific form of the LS operation we have chosen
~Fig. 3! imposes restrictions on the geometry of the system
simulated: with normal periodic boundary conditions the sys-
tem must comprise integral multiples of 6 close-packed
planes. It is possible to avoid this restriction by using more
elaborate boundary conditions @14#, but we chose to avoid
this complication and simulate systems comprising 63
5216, 12351728, and 18355832 spheres. Simulations were
performed at two densities, namely, @see Eq. ~7!# r˜
50.7778 @47#, and r˜51, the close-packed limit.
The maximum step size for displacement updating was
chosen so as to minimize the autocorrelation time of the
overlap order parameter @Eq. 16#. We found a maximum step
size of 0.13D produced the best results at r˜50.7778, while a
value close to unity was found to be appropriate in the close-
packed limit, in the representation ~and scaled units! given in
Eq. ~10a!.
A significant proportion of our simulation time was de-
voted to the process of weight-determination. For our largest
system we used 106 Monte Carlo sweeps ~MCS! to generate
a first ~TP! estimate of the weights, with a further 53106
MCS devoted to weight-refinement using VS.
The free-energy differences of interest were then deter-
mined by further simulations in the multicanonically
weighted ensemble. For each system ~density, and size! we
performed a series of runs each long on the scale of the
autocorrelation time of the overlap order parameter. Each of
these runs then provides an independent estimate of the
~logarithm of the! probability ratio required @Eq. ~18!#. The
standard deviation of these estimates provides a basis for
assigning an associated statistical uncertainty. Implementing
this stage required simulation times ranging from ;2.5
3108 MCS for N5216 to ;43107 MCS for N55832.
IV. RESULTS
A. The effects of the representation
As discussed in Sec. II C the LS operation can be imple-
mented with different choices of representation of the lattice
mapping or the particle displacements @45#. The efficiency of
a lattice mapping is measured ~inversely! by the equilibrium
overlap count. Table I shows results for a variety of map-
PRE 61 LATTICE-SWITCH MOpings, chosen to expose the different factors that control the
mapping efficiency. Mapping number 1 is the one describedin Fig. 3, and used throughout this work: the notation (0,
2 tW ,1 tW) signifies that the three pairs of planes counting from
the top of Fig. 3 are translated respectively by 0,2 tW , and 1 tW .
A similar convention is used to label mappings 2 and 3. In
mapping 4 ~‘‘random-plane’’! an hcp configuration is gener-
ated by taking an fcc configuration and restacking its close-
packed planes in a random order, in an hcp pattern. In map-
ping 5 ~‘‘random-site’’! an hcp configuration is generated by
mapping the particle displacements in an fcc configuration
randomly on to the sites of an hcp lattice.
The random-site mapping ~number 5! shows the largest
overlap count. One can account for its value, rather well, by
regarding the particle displacements as isotropic, Gaussian,
and independent of structure @52#, and estimating the prob-
ability that two particles associated with nearest-neighbor
sites, and with displacements drawn randomly from this dis-
tribution, will overlap.
Using the random-plane mapping ~number 4! cuts the
overlap count by a factor of ~a little more than! 2 with re-
spect to random site. This efficiency gain simply reflects the
fact that of the 6N potential overlaps between near-
neighbors, only the 3N associated with neighbors in different
~but adjacent! planes can now contribute.
Mapping 3 simply generates one fcc configuration from
another ~it is useful only because it is informative!: its over-
lap count is cut by a further factor of 2. This reflects the fact
that this mapping ~similar to mappings 1 and 2! moves close-
packed planes in pairs, thus guaranteeing no overlaps be-
tween the two members of each pair.
Mappings 2 and 1 show further—smaller but still practi-
cally useful—cuts in the overlap count. The origin of these
gains is more interesting. It is clear that they must reflect the
size of the translation vector used: mappings 1 through 3
differ only in this respect. This vector controls the extent of
the shear which the mapping introduces between successive
pairs of planes. The following interpretation seems reason-
able. The displacement patterns in adjacent planes will be
correlated to some extent, with undulations in one surface
~the z components of the displacements! matched to undula-
tions in its neighbor. The smaller the shear, the more closely
these undulations will remain matched to one another ~in the
conjugate configuration!, and the smaller the overlap count.
With increasing shear, this advantage is lost and the behavior
should ~and indeed does! approach the limit ~one quarter of
the overlap count for mapping 5! one would expect in the
absence of such correlations. The fact that this ‘‘approach’’
TABLE I. The efficiency of different lattice mappings ~for N
51728 and r˜50.7778), as measured by the number of overlaps
~per sphere! that they generate. Refer to the text for details.
mapping description effect m5M/N
1 (0,2 tW ,1 tW) fcc→hcp 0.150~1!
2 (0,2 tW ,22 tW) fcc→hcp 0.183~1!
3 (0,3 tW ,23 tW) fcc→fcc 0.194~1!
4 random-plane fcc→hcp 0.373~2!
5 random-site fcc→hcp 0.820~3!
913TE CARLO METHODis already apparent in the performance of mapping 2 is con-
sistent with the fact that the measured correlation length of

˜nt
tothe surface undulations at the density concerned is found to
be close to the magnitude of the translation vector tW .
These results help to clarify the factors which control the
overlap count of the mapping ~number 1! we have actually
used. It is tempting to attribute the overlaps to the fact that
the LS (fcc→hcp, say! maps each particle from an environ-
ment in which adjacent close-packed planes have different
stacking labels ~A and C, say! to one in which they have the
same label ~C, say!. The results for mappings 1–3 show that
it would be misleading to think this way. The overlaps sim-
ply reflect the numbers of particles that ‘‘see’’ a new adja-
cent close packed plane ~irrespective of its label!, and the
extent to which it is ‘‘new.’’ This is the reason for the simi-
packed planes in a system of 216 spheres at r50.7778, in the
equilibrium hcp and fcc macrostates, and in the gateway (M50)
macrostate. The separation is measured with respect to the equilib-
rium separation d0 and is expressed in units of the sphere separation
d @57#.
FIG. 7. ~a! The mean value of the separation d between adjace
macrostates of different M. The separation is measured with respect
~see Fig. 3! are translated together by LS; category ~ii! planes are trans
the c/a ratio @58# in a constant-pressure ensemble @with the same paramlibrium fcc, equilibrium hcp, and gateway (M50) regions.
The macrostates corresponding to the equilibrium crystal
structures have similar, near-Gaussian, d distributions. In
contrast, for the gateway macrostate the distribution is bimo-
dal: in this macrostate, some planes are systematically
moved closer to one another, while ~in equal measure! others
are shifted apart. On closer examination one finds that it is
the planes which are translated together by the LS @e.g., the
pair of planes marked ~i! in Fig. 3# that fall into the first
category, while the planes that are translated differently by
the LS @e.g., the pair of planes marked ~ii! in Fig. 3# fall into
the second. The evolution, with M , of the mean plane sepa-
ration ~for both categories! is shown in Fig. 7~a!. The behav-
ior thus unearthed is entirely reasonable. The LS operation
can only create overlaps between neighboring planes which
are translated by different amounts ~sheared with respect to
one another!. The algorithm resolves the task set by the bias
towards M50 by moving these pairs of planes ~the ones
vulnerable to overlaps! further apart, at the expense of a
close-packed planes in a system of 216 spheres at r˜50.7778, for
the equilibrium separation d in units of d @57#. Category ~i! planeslarity between the overlap counts for the two structures ~Sec.
IV C!. It shows, moreover, that any simple @53# tuning of the
displacement representation ~the choice of T matrix! is likely
to be of no advantage here @54#.
B. How it works: the gateway configurations
A LS operation will work ~be accepted as a MC move!
only when launched from a small subset of the configura-
tions actually visited: these, by definition, are the ‘‘gateway
configurations.’’ As noted earlier, one could identify a priori
configurations ~those characterized by ‘‘small enough’’ dis-
placements! which fall into this set. But we have elected,
rather, to let the system ~the algorithm! identify them on the
basis of their defining characteristic—that they have zero
overlap order parameter M @28#. It is then interesting to
investigate the microscopic characteristics of the configura-
tions picked out by this constraint. Figure 6 shows the dis-
tribution of the separation d between adjacent close-packed
(x-y) planes @55#, for M macrostates corresponding to equi-FIG. 6. Distribution of the separation d between adjacent close-
914 PRE 61BRUCE, JACKSON, ACKLAND, AND WILDING0
lated through different amounts by LS. ~b! The evolution with M of
eters as ~a!#. The horizontal line marks the ideal-close-packed value.

Ncompression of the others @56#. In simulations conducted at
constant pressure this effect ~still present! is supported by a
second. Figure 7~b! shows that the algorithm now exploits
the additional degrees of freedom ~the shape of the simula-
tion cell! to locate gateway states with values of the c/a ratio
enhanced above the ideal close-packed value @58#. Again, the
advantages with respect to the switch are clear.
It is tempting to say that the sampling is intelligent. In any
event it is clear that the algorithm locates and utilizes con-
figurations which it would be difficult to exploit explicitly in
the design of the switch operation.
C. Entropies of crystalline structures
The essential output of a LS simulation is in the form of
the normalized probability distribution of the overlap-order
parameter, reweighted to remove the bias in the multicanoni-
cally weighted distribution actually measured. Figure 8
shows the results for this distribution ~at r˜50.7778) for
three different N values. As one would expect the distribu-
tions each comprise two peaks ~one associated with each
phase! each of which is nearly Gaussian @59# and sharpens
with increasing N @60#. Note the close correspondence be-
tween the equilibrium overlap counts for the two structures.
This result is not required by definition, or any obvious sym-
metry. Rather it should be seen as a further manifestation
~the smallness of the entropy difference between the phases
is the prime one! of the similarity of the local particle envi-
ronments in the two structures.
The relative weights of the two peaks is a direct measure
of the difference between the entropies of the two structures
@Eqs. ~11!, ~12!, ~18!#. Since the entropies are extensive the
ratio of the peak weights grows exponentially with N @61#;
FIG. 8. The probability distribution of the overlap order param-
eter per particle, m[M/N , for systems of three different N values
at r˜50.7778. The lines provide Gaussian guides to the eye; the
statistical uncertainties on the data points are smaller than the sym-
bol size. The entropy difference is identified from the logarithm of
the ratio of the integrated weights of the two peaks. The hcp peak
for the largest system is not visible on this scale.
PRE 61 LATTICE-SWITCH MOthe fact that ~in this case, at least for our smaller systems! the
two peaks can even be displayed on the same scale is areflection of the exceptionally delicate balance between the
two entropy densities.
Figure 8 allows one to see that fcc is the thermodynami-
cally preferred structure. This conclusion is expressed quan-
titatively in the results gathered in Table II. Our results at
r
˜
50.7778 correct those of our earlier work @30#, as ex-
plained in Sec. III A. They are in full accord with the results
~both LS- and IM-based! reported by Pronk and Frenkel
@31#. The close correspondence between the results for N
51728 and N55132 confirms that the former system is al-
ready representative of the thermodynamic limit. Table II
also shows the results of our studies at the close-packed
limit, using the hard-dodecahedron representation @Eq.
~10a!#. Our results seem at variance with the IM-based result
of Woodcock @64#, even allowing for the large uncertainty
attached to that result. They are close to those ~based on LS!
reported by Mau and Huse @14#. But the differences ~for the
smaller systems, particularly! appear to be statistically sig-
nificant @65#. Figure 9 gives an alternative view of these
results. It utilizes the parametrization of the measured pres-
sure difference between the two phases provided by Speedy
@66# to determine the entropy difference, as a function of
density, given the entropy difference at a chosen reference
density; we have used the results of the present work at r˜
50.7778.
Table III shows the results of our studies in the constant
pressure ensemble. The quantity of interest here is the differ-
ence between the Gibbs free energy densities at the chosen
pressure, which follows from the relevant distribution with
the aid of Eq. ~15!. In fact the Gibbs free energy density
TABLE II. The difference in the entropy densities of the fcc and
hcp structures, Ds[Ds fcc,hcp @Eq. ~11!#; the associated uncertainties
are in parenthesis. The results of the present work ~PW! supercede
those of Ref. @30#. The results of Ref. @64# supercede those of Ref.
@10#. IM stands for integration method; SM is the lattice shear
method of Refs. @14,63#.
r/rCP N Ds(1025kB) Method Ref.
0.731 512 85 ~10! SM @14#
0.736 12000 230 ~100! IM @64#
0.736 12096 87 ~20! IM @62#
0.739 512 90 ~4! LS @14#
0.7778 216 132 ~4! LS @31#
0.7778 1728 112 ~4! LS @31#
0.7778 1728 113 ~4! IM @31#
0.7778 216 133 ~3! LS PW
0.7778 1728 113 ~3! LS PW
0.7778 5832 110 ~3! LS PW
1.0 12000 260 ~100! IM @64#
1.0 512 110 ~20! SM @14#
1.0 64 91 ~5! LS @14#
1.0 216 107 ~4! LS @14#
1.0 512 119 ~3! LS @14#
1.0 1000 113 ~4! LS @14#
1.0 216 131 ~3! LS PW
1.0 1728 125 ~3! LS PW
915TE CARLO METHODdifference Dg for a given pressure, and the entropy density
difference Ds at a physical density that is the thermodynamic

KLconjugate of that pressure for one of the phases, differ ~in
magnitude @67#! by terms that are second order in the pres-
sure difference between the two phases. That pressure differ-
ence is extremely small @66#, as is the difference between the
measured densities of the two structures ~Table III!. In these
circumstances one would expect the magnitude of Dg to fall
on the Ds plot in Fig. 9; and indeed, within the residual
uncertainties, it does.
V. DISCUSSION: REVIEW AND PROSPECTS
In the work described here we have been concerned both
with a system of long-standing interest—the hard sphere
crystal—and a method —lattice-switch Monte Carlo— with
potentially wide applicability. We divide our concluding dis-
cussion accordingly.
The full agreement between the present work and that of
Ref. @31# leaves little doubt that the equilibrium entropy dif-
ference between the two close-packed structures has finally
been established securely and with high precision—at least
at one density. Although a small discrepancy with respect to
the results of Ref. @14# remains, the accord of our close-
packed limit results with those established using pressure
difference measurements @66# suggests that the curve in Fig.
9 provides a relatively complete and trustworthy picture of
the density dependence.
Notwithstanding the simplicity of the model, these results
do have implications for experimentally realizable systems.
The immediate relevance to atomic systems is tenuous @68#,
but the model has been widely used to account for the be-
havior of assemblies of ‘‘hard,’’ ‘‘spherical’’ colloidal par-
ticles @32#. Since the predicted entropy-density difference is
so small there are potentially many ways ~residual interac-
tions between the spheres; polydispersity! in which the ap-
plicability of the theory may be compromised. But, of these,
FIG. 9. The difference in the entropy densities of the fcc and
hcp structures, Ds[Ds fcc,hcp @Eq. ~11!#, as a function of reduced
density r˜ . The data points are as given in Table II. The solid line is
the result of an integration of the pressures of the phases @66#. Note
that this line passes through our result at r˜50.7778 by construc-
tion.
916 BRUCE, JACKSON, ACit seems that the most significant issues to be addressed are
to do with scales—length and time.First, the lengths. In the experiments reported in Ref. @69#
the colloidal particles have diameters of order 1027 m and
the samples comprise crystallites with linear dimensions of
order 1025 m. The number of particles in such a crystallite
(N;106) is large compared to those in our simulation,
which is ~as we have seen! sufficient to allow us to deduce
properties of the thermodynamic limit. But it is not large
enough to guarantee that the behavior displayed will actually
be that of the thermodynamic limit. To see this—and its
principal implications—one needs to consider the stability of
the perfect fcc crystal with respect to hcp-type stacking
faults. Following Ref. @69# we may introduce a parameter a
@70# measuring the probability that a chosen close-packed
plane sits within an fcc environment as distinct from the hcp
environment. A simple argument ~Appendix A! using the
pseudospin parametrization of stacking patterns provided in
Ref. @14# then yields the result
a5
1
2 S 11tanhF
N
’
Ds
2 G D , ~24!
where N
’
is the number of particles in a close packed layer
and Ds ~a function of r˜ ) is the fcc-hcp entropy difference
per particle, as given in ~and in the units of! Table II. The
thermodynamic ideal (a51) is thus realized only to the ex-
tent that N
’
Ds is large compared to unity. For the length
scales given above, N
’
Ds.1. The obvious implications are
qualitatively consistent with the observations reported in
Ref. @69# which show a values ~deduced from Bragg scat-
tering intensities! ranging from 0.5 @signaling essentially
random-hexagonal-close packing, ~rhcp!# through to a50.8.
The observed spread in a values reflects, presumably, the
issue of time scales. The smallness of the entropy difference
~which supplies the kinetic driving force towards the equilib-
rium state! suggests that the equilibrium behavior will be
observed only in samples which are grown sufficiently
slowly and ~or! given sufficient time for subsequent anneal-
ing @71#. The results of Ref. @69# do indeed suggest a corre-
lation between observed a value and the slowness of the
growth process. Experiments done in microgravity @72#,
where growth processes are greatly accelerated, yield essen-
tially randomly close-packed crystals.
Now let us turn to the lattice-switch method. There are
two questions here. One, does the method represent a signifi-
cant advance with respect to existing methods? Two, is it
generally applicable?
The main alternative method ~the benchmark against
which others need to be assessed! is probably integration
TABLE III. The difference in the gibbs free energy densities of
the fcc and hcp structures Dg[Dg fcc,hcp @Eq. ~15!#; the associated
uncertainties are in parenthesis. P! gives the pressure @42# in units
of kBT/D3.
P!(D23)
r
˜ hcp r˜ fcc N Dg(1025kBT)
14.58 0.7776~1! 0.7775~1! 216 2113 ~4!
14.58 0.7770~3! 0.7774~2! 1728 2112 ~3!
PRE 61AND, AND WILDINGalong a reference path, of which the work reported in Ref.
@31# represents, to our knowledge, the most refined example.

If one compares the two techniques ~LS and IM! on the basis the Royal Society of Edinburgh. A.D.B. acknowledges help-
PRE 61 917LATTICE-SWITCH MONTE CARLO METHODof precision-for-computational-buck there seems to be no
clear winner in the hard-sphere studies to date: Ref. @31#
reports calculations using both methods that achieve compa-
rable levels of precision on the basis of comparable compu-
tational time. But one should note that the entropy difference
ultimately determined is some four orders of magnitude
smaller @73# than the separate entropies of the two phases,
determined by IM. One can see this as a testimony to the
care with which the recent IM studies have been carried out;
or ~as suggested in Sec. I! as a strong indicator that another
approach using an interphase path is called for. There are
also two other counts—both somewhat subjective—on
which we suggest that the LS approach is superior. First, it
seems to us relatively illuminating ~by comparison with IM!
to read-off the result for a free energy difference directly
from a figure similar to Fig. 8 which shows what it means.
Secondly it also seems to us that LS wins in regard to the
transparency of the uncertainties to be attached to its results.
The LS error bounds represent purely statistical uncertainties
associated with the measurement of the relative weights of
two distribution peaks. The IM error bounds have to aggre-
gate the uncertainties associated with different stages of the
integration process.
As regards the second question, we expect that the
method will, with appropriate extensions, be widely appli-
cable. The first extension must clearly be to accommodate
soft potentials. The LS operation will then need gateway
configurations in which the energies of the two structures
~measured with respect to their ground-state energies—or in-
deed any fixed reference energy! are closely matched @43#.
The ‘‘overlap order parameter’’ will need to be redefined
accordingly. With no more than this degree of elaboration
the method should be applicable immediately to investigate
the widespread ‘‘competition’’ between fcc and hcp ordering
in the phase behavior of the elements @74#.
More generally, moving beyond the space of fcc-hcp
structures, the choice of lattice-to-lattice mapping will re-
quire some thought. Mappings which preserve the relative
positions of significant subsets of the particles ~the analog of
the close-packed planes! are likely to be optimal. The license
to choose ones representation of the displacements ~Sec.
II C! may also prove useful. Simple transformations @53# will
help if the mapping takes particles between environments in
which the spectrum of single-particle displacements is sig-
nificantly different. In such cases one might envisage using a
MC-annealing procedure to refine the choice of representa-
tion. The use of normal coordinates has some advantages
here—but possibly not enough to offset the fact that the in-
teraction potential is nonlocal when expressed in Fourier
space.
ACKNOWLEDGMENTS
A.N.J. acknowledges the support of the EPSRC. N.B.W.
acknowledges the financial support of the Royal Society
~Grant No. 19076!, the EPSRC ~Grant No. GR/L91412! andful discussions with Dr. Mike Evans.
APPENDIX A: DISPLACEMENT ENTROPY VERSUS
STACKING ENTROPY
Consider a system of N hard spheres arranged in N
i
close-
packed layers of N
’
particles. Following Ref. @14# one may
conveniently index each of the close-packed layers with a
pseudo-spin ~Ising-like! variable s , where s i51 signifies
that layer i has an fcc environment ~the two immediately
adjacent layers are not aligned with one another! while s i
521 implies an hcp environment ~adjacent planes are
aligned with one another!. The probability of a particular
stacking sequence
$
s
%
~if these variables may properly be
regarded as annealed! then satisfies
ln P
~
$
s
%
uN ,V !5S~N ,V ,
$
s
%
!1const, ~A1!
where S(N ,V ,
$
s
%
) measures the entropy associated with the
configurations ~displacements! consistent with the particular
structure
$
s
%
. Following Ref. @14# this entropy ~we will refer
to it here as ‘‘displacement entropy’’! can usefully be written
in the form of an expansion:
S
~
N ,V ,
$
s
%
!5Ns01N’h(
i
s i1N’J(
^
i j
&
s is j1••• .
~A2!
The expansion is effectively ordered in the range of the en-
tropic interlayer ‘‘interactions’’: the ellipsis represents con-
tributions from interactions ~microscopically, displacement-
displacement correlation functions! extending over more
than 4 layers. The analysis of Ref. @14# indicates that the
series converges quickly, except close to melting. If we ne-
glect the interaction terms altogether we may make the iden-
tification
h5
1
N
’
N
i
@
S
~
N ,V ,
$
s511
%
!2S~N ,V ,
$
s521
%
!#
5
Ds fcc,hcp
2 ~A3!
and, from Eq. ~A1!,
^
s
&
5
1
N (
$
s
%
,i
P
~
$
s
%
uN ,V !s i5tanh@N’h#
5tanh
F
N
’
Ds fcc,hcp
2 G ~A4!
from which Eq. ~24! follows. The correspondence with a 1D
paramagnet is clear. The familiar competition ~between ori-
entation energy and entropy! is played out here as a compe-
tition between displacement entropy and stacking entropy,
with N
’
playing the role of an inverse temperature.

KL@1# G. Jacucci and N. Quirke, in Complete Simulation in the Phys-
ics and Chemistry of Solids edited by C.R.A. Catlow and W.
C. Mackrodt, Vol. 166 of Lecture Notes in Physics, ~Springer,
Berlin, 1982!, p. 38.
@2# Most, perhaps all, of the methods referred to here may also be
realized within the framework of molecular dynamics.
@3# K. Binder, Rep. Prog. Phys. 60, 487 ~1997!.
@4# D. Frenkel, in Molecular-Dynamics Simulation of Statistical-
Mechanical Systems, edited by G. Ciccotti and W.G. Hoover
~North-Holland, Amsterdam, 1986!, p. 151.
@5# 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.
@6# A path is defined by a sequence of values of either some mac-
roscopic observable or some parameter l , which controls a
thermodynamic field or a model parameter.
@7# In some instances one needs to make use of separate reference
systems for each phase.
@8# A full survey is not appropriate here; we note only some ex-
amples. References @9# and @10# both use multistage ~integra-
tion! methods involving reference paths—the former to an Ein-
stein solid, and the latter to a single-occupancy-cell model of
an ideal gas–to study the fcc-hcp phase behavior of hard
spheres. Reference @11# describes a study of the fcc-bcc phase
behavior of Lennard-Jones systems, using a multistage ~over-
lap! method applied to a nonphysical interphase path, with the
l parameter @6# indexing a configurational energy that inter-
polates between those of the two structures. This method has
been widely used ~see, e.g., Ref. @12#!. Reference @13# de-
scribes a physical interphase path linking NaCl and CsCl struc-
tures, but does not utilize the path for a free-energy-difference
calculation. The ‘‘lattice-shear’’ method described in Ref.
@14#, and applied to hard spheres, is a natural refinement of
Ref. @13# in which a single-stage sampling procedure is used to
explore a path constructed out of a series of mutually sheared
lattices that interpolate between the lattices of interest.
@9# D. Frenkel and A.J.C. Ladd, J. Chem. Phys. 81, 3188 ~1984!.
@10# L.V. Woodcock, Nature ~London! 384, 141 ~1997!.
@11# A. Rahman and G. Jacucci, Nuovo Cimento D 4, 357 ~1984!.
@12# M.C. Moody, J.R. Ray, and A. Rahman, J. Chem. Phys. 84,
1795 ~1985!; B. Kuchta and R.D. Etters, Phys. Rev. 47, 14 691
~1993!.
@13# M. Parrinello and A. Rahman, J. Phys. ~Paris! 42, 511 ~1981!.
@14# S.-C. Mau and D.A. Huse, Phys. Rev. E 59, 4396 ~1999!.
@15# B.A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 ~1992!; B.A.
Berg, J. Stat. Phys. 82, 323 ~1996!.
@16# A.P. Lyubartsev, A.A. Martsinovski, S.V. Shevkunov, and
P.N. Vorontsov-Velyaminov, J. Chem. Phys. 96, 1776 ~1992!.
@17# E. Marinari and G. Parisi, Europhys. Lett. 19, 451 ~1992!.
@18# G.M. Torrie and J.P. Valleau, Chem. Phys. Lett. 28, 578
~1974!.
@19# Since, by definition, the end points of an interphase path lie in
different phases, the associated macroscopic property may rea-
sonably be described as an order parameter.
@20# The nature of the configurations visited along a path is not
generally predictable a priori from the choice of order param-
eter. One may say only that the configurations sampled at a
given point on the path will be those which have measurable
canonical probabilities conditioned on that macrostate, and
which are, moreover, accessible on the relevant time scales.
@21# The interfacial free energy emerges as a by-product—and in
918 BRUCE, JACKSON, ACsome cases @22# may actually be the principal focus of interest.
@22# B.A. Berg, U.H.E. Hansmann, and T. Neuhaus, Z. Phys. B:Condens. Matter 90, 229 ~1993!.
@23# N.B. Wilding, Phys. Rev. E 52, 602 ~1995!.
@24# B. Grossmann, M. Laursen, T. Trappenberg, and U.J. Wiese,
Phys. Lett. B 293, 175 ~1992!.
@25# It might work acceptably well if the dynamics of the interface
between the two phases is favorable: systems with martensitic
phase transitions may fall into this category: Z. Nishiyama,
Martensitic Transformations ~Academic, New York, 1978!.
Note also that the special case in which the structural phase
transition involves no change of symmetry can be handled
within the standard multicanonical framework: see Ref. @26#.
@26# G.R. Smith and A.D. Bruce, Phys. Rev. E 53, 6530 ~1996!.
@27# The conjugate of a given configuration is the configuration
associated with the same set of displacements attached to the
other set of lattice vectors.
@28# The general defining characteristic of the gateway configura-
tions is that a LS operation, launched from within this set, will
be accepted with a probability sufficient to make the attempt
worthwhile.
@29# A configuration is energy matched ~to its conjugate! if the
difference between its energy and that of its conjugate is small,
on the scale of kBT .
@30# A.D. Bruce, N.B. Wilding, and G.J. Ackland, Phys. Rev. Lett.
79, 3002 ~1997!.
@31# S. Pronk and D. Frenkel, J. Chem. Phys. 110, 4589 ~1999!.
@32# P.N. Pusey, in Liquids Freezing and the Glass Transition, ed-
ited by J.P. Hansen, D. Leveesque, and J. Zinn-Justin
~Elsevier, Amsterdam, 1991! p. 763.
@33# In its most general context the LS is the basis for estimating
the differences between the free energy ~Helmholtz or Gibbs!
of two crystalline phases; in the case of hard sphere systems, at
constant density, this free energy difference is purely entropic.
@34# M. Born and K. Huang, Dynamical Theory of Crystal Lattices
~Clarendon, Oxford, 1968!.
@35# The hcp ‘‘lattice’’ is not a ‘‘Bravais lattice’’; see e.g., Ref.
@36#, p. 79.
@36# N.W. Ashcroft and N.D. Mermin, Solid State Physics ~Saun-
ders, Philadelphia, 1976!.
@37# W.G. Hoover and F.H. Ree, J. Chem. Phys. 49, 3609 ~1968!.
@38# W.G. Rudd, Z.W. Salsburg, A.P. Yu, and F.H. Stillinger, J.
Chem. Phys. 49, 4857 ~1968!.
@39# The divergence is an artifact of the classical character of the
model.
@40# As detailed in Ref. @38# the argument leading to Eq. ~10a!
actually entails a rescaling of the displacement coordinates by
an e-dependent factor.
@41# The ‘‘visualization’’ exercise has its limitations. The dodeca-
hedra ~and the lattice spacing! should be thought of as ‘‘infi-
nitely large’’ compared to the mean separation of adjacent
faces.
@42# We subsume a factor of kBT into the ‘‘pressure’’ P!.
@43# If the interparticle potential is not of the hard-sphere form @Eq.
~1!# the ‘‘perfect crystal’’ configurations ~classical ground
states! of the two structures will generally have different ener-
gies. But one may handle the effects of this energy mismatch
quite simply, by attaching different multicanonical weights to
the gateway configurations of the two structures.
@44# The sign convention here has no deep significance. Defining
PRE 61AND, AND WILDINGM so that it has different signs in the two phases simply
allows us to make a visually clear distinction ~Fig. 8! between

the contributions which each phase makes to the M distribu-
tion.
@45# In fact, within the framework of periodic boundary conditions,
the second generalization subsumes the first. Thus the various
lattice-to-lattice mappings discussed in Sec. IV A can all be
thought as a single mapping but with different hcp T matrices
@Eq. ~19!#, chosen to interchange appropriate displacements.
@46# In the absence of any phase-defining constraint on the configu-
52d0 /D˜ with d0 the mean separation between close-packed
planes and D˜ the mean separation between neighboring sphere
centers in close packed planes. In the hcp structure c/a5A 83
51.63299 if the packing is ideal @36#; in the fcc structure it
has this value by symmetry. In our constant-pressure simula-
tions ~for N51728,P!514.58) we found
@
c/a
# fcc
51.6333(3) and
@
c/a
#hcp51.6332(3). Evidently, whatever
difference there is between the true values of these quantities is
PRE 61 919LATTICE-SWITCH MONTE CARLO METHODrational integral specified in Eq. ~5!, the integral extends over
all configurations compatible with the boundary conditions.
@47# This density value was chosen to coincide with one of those
studied in Ref. @9#.
@48# This conclusion is in qualitative accord with that of Ref. @14#,
although we find the effects of next-neighbor encounters to be
somewhat smaller than is reported there.
@49# D. Frenkel, in Computer Modelling of Fluids, Polymers and
Solids, edited by C.R.A. Catlow, C.S. Parker, and M.P. Allen
~Kluwer Academic, Dordrecht, 1990!, p. 83.
@50# G.R. Smith and A.D. Bruce, J. Phys. A , 28, 6623 ~1995!.
@51# M. Fitzgerald, R.R. Picard, and R.N. Silver, Europhys. Lett.
46, 282 ~1999!.
@52# These approximations are surprisingly good, except close to
melting. See inter alia D.A. Young and B.J. Alder, J. Chem.
Phys. 60, 1254 ~1974!; J. Piasecki, and L. Peliti, J. Phys. A 26,
4819 ~1993!; R. Ohnesorge, H. Lo¨wen, and H. Wagner, Euro-
phys. Lett. 22, 245 ~1993!.
@53# By ‘‘simple’’ we mean a T matrix that preserves locality,
merely acting to rotate or reflect individual displacements.
Nonlocal transformations which mix the displacement vectors
might reduce the overlap count further.
@54# The fact that the two phases have such similar entropies tells
us that in principle there must be some representation of the
LS in which virtually all realizable configurations of one phase
transform into realizable configurations of the other: for this
representation the white islands of Fig. 1~b! would fill the
single phase regions; the effective length of the interphase path
would be short, delivering results of essentially arbitrary pre-
cision. We think that this representation is not ‘‘simple.’’
@55# More precisely, the difference between the z components of
the center of mass coordinates of these planes.
@56# The gateway configurations have further distinctive features
~with respect to their equilibrium counterparts!: the root-mean-
square particle displacement, in the z direction, is reduced; and
the in-plane correlation length of these displacements ~the un-
dulations of the close packed planes! is enhanced along the tW
direction.
@57# We define d5D(r˜21/321), giving the shortest distance be-
tween two sphere surfaces in the ideally close-packed perfect
crystal configuration.
@58# We define a c/a ratio, for both hcp and fcc structures by c/anot easily resolved by single-phase averages.
@59# A Gaussian parametrization of the peaks provides a first ap-
proximation to a set of multicanonical weights—but one which
significantly underestimates the entropic barrier against the LS.
@60# Note that in Fig. 8 the horizontal axis is the intensive variable
m[M/N .
@61# The fact that the peak-weight difference diverges exponen-
tially fast with N presents no computational problem: the mul-
ticanonical procedure is designed to cope with ~and quantify!
differences of this scale.
@62# P.G. Bolhuis, D. Frenkel S.-C. Mau, and D.A. Huse, Nature
~London! 388, 235 ~1997!.
@63# The ‘‘shear’’ implementation described by Mau and Huse @14#
is ~they report! substantially harder to implement than the
switch method which we discuss here, and which they use in
most of their work.
@64# L.V. Woodcock, Nature ~London! 388, 236 ~1997!.
@65# To provide one independent check of the close-packed limit
code we determined the amplitude of the leading nontrivial
term in the expansion of the pressure about the close-packed
limit @38# in the fcc structure, finding agreement to 4 signifi-
cant figures with results reported in Ref. @66#.
@66# R.J. Speedy, J. Phys.: Condens. Matter 10, 4387 ~1998!.
@67# But Ds and Dg are of opposite sign.
@68# Y. Choi, T. Ree, and F.H. Ree, J. Chem. Phys. 99, 9917
~1993!. show that predictions for the phase diagram of a Len-
nard Jones solid depend extremely sensitively on the fcc-hcp
hard-sphere entropy difference.
@69# P.N. Pusey, W. van Megen, P. Bartlett, B.J. Anderson, J.G.
Rarity, and S.M. Underwood, Phys. Rev. Lett. 63, 2753
~1989!.
@70# Note our double use of a as both a stacking-type probability
and a structure label.
@71# Reference @31# estimates the rhcp-to-fcc annealing time for a
hard-sphere-colloidal crystal of linear dimension 1023 m to
lie in the range ‘‘months to years.’’
@72# J. Zhu, M. Li, R. Rogers, W. Meyer, R.H. Ottewill, W.B.
Russel, and P.M. Chaikin, Nature ~London! 387, 883 ~1997!.
@73# See the full accounts of the IM method provided in the original
studies @9#.
@74# D.A. Young, Phase Diagrams of the Elements ~University of
California, Berkeley, 1991!.