Lattice-switch Monte Carlo simulation for binary hard-sphere crystals
A. N. Jackson and G. J. Ackland*
SUPA, School of Physics, University of Edinburgh, Edinburgh EH9 3JZ, Scotland, United Kingdom

Received 3 August 2007; published 7 December 2007
We show how to generalize the lattice-switch Monte Carlo method to calculate the phase diagram of a binary
system. A global coordinate transformation is combined with a modification of particle diameters, enabling the
multicomponent system in question to be explored and directly compared to a suitable reference state in a
single Monte Carlo simulation. We use the method to evaluate the free energies of binary hard-sphere crystals.
Calculations at moderate size ratios
=0.58 and 0.73 are in agreement with previous results, and confirm AB2
and AB13 as stable structures. We also find that the AB
CsCl structure is not entropically stable at the size
ratio and volume where it has been reported experimentally, and therefore that those observations cannot be
explained by packing effects alone.
DOI: 10.1103/PhysRevE.76.066703 PACS number
s: 02.70.c, 05.10.Ln, 65.40.Gr
I. INTRODUCTION
The most efficient packing of hard spheres is one of the
most easily posed questions in physics. For monodisperse
spheres, it forms the basis of Kepler’s conjecture 1,2, and
the solution is any one of many degenerate close-packing
arrangements. In physical systems such as opals and colloi-
dal suspensions, the problem is generalized to calculation of
the highest-entropy structure of a given density. Away from
the degenerate close-packing limit, face-centered cubic turns
out to be the most stable structure, and its exact entropy
advantage has been accurately calculated 3. A further gen-
eralization, which we consider here, is the packing of spheres
of different radii: the so-called binary hard-sphere system.
In order to calculate stable packings for binary hard
spheres, one must evaluate the free energies of candidate
phases as a function of density, composition, and size ratio,
so that the entire phase diagram can be constructed by mini-
mization of the total free energy.
A common approach to free energy calculation is to use
thermodynamic integration. This involves the construction of
a reversible path between the structure of interest and some
reference system. By performing a simulation or set of simu-
lations that traverse that path, one can estimate the free en-
ergy difference between the reference and target systems by
integrating the derivative of the free energy along that path.
For example, the integration technique of Frenkel and Ladd
4 uses an Einstein crystal as a reference system, and has
been successfully applied to a wide range of target systems,
including binary hard-sphere solids 4–7. However, the in-
tegration process introduces an intrinsic source of systematic
error. Small systematic errors can accumulate at each step in
the integration, and integrating across the singularity at a
phase transition causes problems. Of course, careful applica-
tion can minimize these effects, but it remains preferable to
have a technique where systematic errors are more easily
quantifiable. The most accurate free energy differences in
hard-sphere systems are obtained from lattice-switch Monte
Carlo simulation 3, which allows two phases to be sampled
directly during a single Monte Carlo simulation, mapping
between them using a global coordinate transformation. Pro-
vided some biased sampling method ensures that the phase
transformation move is accepted, the free energy difference
between the structures can be determined directly from the
relative probability of finding the simulation in either phase.
Finite-size effects are the only source of systematic errors in
this estimation process. This approach was first used to study
the solid state phase behavior of monodisperse hard spheres
3,8. The method was subsequently developed for different
energy models
the “Hamiltonian switch” and off lattice

“phase switch”. It has been applied to the freezing transi-
tion for hard spheres in Ref. 9,10 and also to soft potentials
in Ref. 11–13.
Here we investigate whether this methodology can be ap-
plied successfully to multicomponent systems by generaliz-
ing the lattice-switch approach to examine the structural
phase behavior of binary hard spheres. The binary hard-
sphere system is well studied, and therefore the validity of
the generalized lattice-switch method can be verified by di-
rect comparison with the results from previous experimental
and theoretical work.
Stable phases of the binary hard-sphere system are char-
acterized primarily by the diameter ratio
=B /A, where
the large and small particles are designated as A and B, re-
spectively. At high values of
, from 1.0 down to 0.92, the
system is expected to form a substitutionally disordered crys-
tal, where A and B particles are randomly distributed on a fcc
lattice 14,15. At slightly lower values, in the range 0.92

0.85, this mixed crystal becomes metastable, and the
solid phase separates into fcc crystals of A and B particles
14,15.
Below this, for 0.85
0.5, three different lattices

shown in Fig. 1 have been observed experimentally: AB2,
AB13, and AB
CsCl. AB13 and AB2 have been observed in
Brazilian gem opals
colloidal silica crystals 16,17, and
later in
poly-methylmethacylate
PMMA colloidal disper-
sions, both of which closely resemble hard spheres. In
PMMA, AB13 was first seen in 18 for
0.61, followed
by AB2 and AB13 at
0.58 in 19. These two structures
were determined to be unstable above
0.62 in Ref. 20
and below
0.52 in Ref. 21. More recent experiments*G.J.Ackland@ed.ac.uk
PHYSICAL REVIEW E 76, 066703
2007
1539-3755/2007/76
6/066703
10 ©2007 The American Physical Society066703-1

have found AB2 to be stable for 0.60
0.425 and AB13
stable over 0.62
0.485 22.
These complex lattices have also been seen in systems
that bear little resemblance to hard spheres, such as rare gas
solids 23, nanocrystals 24, charge-stabilized colloids 25,
and block copolymer micelles 26. At first, it is surprising
that these structures, particularly the highly complex AB13
superlattice, should form spontaneously in such a wide range
of systems. However, it appears that entropic effects of hard-
sphere packing are sufficient to stabilize these structures.
Previously, the lattice-switch approach has been used for
the free energy difference between two candidate structures
for a single species. Two-species phase behavior is more
complicated, as the stable state of a binary system may con-
sist of coexisting phases of different compositions, at equal
pressures and chemical potentials. Therefore, instead of just
comparing the free energy of pairs of candidate phases, we
must evaluate the absolute free energies of all competing
phases so that the overall phase diagram can be constructed.
The basic approach has been presented in Ref. 27 by Bar-
tlett, who considered only the fluid and pure
face-centered
cubic A or B crystal phases. We shall follow the work of
Eldridge et al. 5–7,20 and Cottin and Monson 28,29 and
also include AB2, AB13, and AB
CsCl as candidate two-
component crystals. In the Eldridge et al. 5 work, the sta-
bility of AB2 and AB13 for binary hard spheres was deter-
mined using the thermodynamic integration method of
Frenkel and Ladd 4, with further details of the AB2 calcu-
lations in Ref. 6 and AB13 in Ref. 7. These results show
general agreement with the experimental systems, but sug-
gest that AB2 is kinetically repressed at the AB2 1:2 stoichi-
ometry, as it appears only at higher concentrations of B
particles 22.
Theoretical studies of AB2, AB13, and AB
CsCl have
also been carried out via cell theory 28,29, and density
functional theory
DFT 30–32. For AB2 and AB13, the
results are generally consistent with both the integration
method and experimental results, although some minor dif-
ferences remain. In particular, the range of stability of AB13
is computed via cell theory to be 0.61
0.54 29,
whereas the integration method finds AB13 to be stable over
the range 0.626
0.474 7. Although the latter work
does not consider the possibility that AB2 is more stable than
AB13 over this range of
, it is notable that the range of
stability determined by integration methods is consistent
with the range of observation of AB13 found in the experi-
mental work presented in 22.
The case of the AB
CsCl structure is somewhat less
clear. The close-packed density of this structure peaks at

=0.732, and so may compete with the slightly denser fcc
structure, and indeed metastable CsCl has been observed ex-
perimentally at
=0.736 in 33. The stability of this struc-
ture has not been ascertained by integration methods: it was
not considered to be a valid candidate structure in Ref. 20.
The cell theory work in Ref. 29 did not find CsCl to be
stable for any
, but DFT results have shown metastability
31, and even stability 32. However, the work in Ref. 32
did not consider full phase separation into coexisting A and B
fcc crystals; they authors consider only a fluid coexisting
with a crystal phase made up of randomly distributed A and
B particles, but rich in A, and the fluid composition appears
to be fixed at 1:1. Interestingly, the AB
CsCl structure has
been stabilized experimentally by using oppositely charged
colloids 25,34–36. Moreover, bcc is the lowest-energy ar-
rangement for charged monodisperse particles, so even if the
AB
CsCl structure is only metastable in general, it is pos-
sible that a small amount of charge transfer between particles
or between particles and fluid is enough to stabilize the
structure.
For even smaller sphere size ratios

0.44, the situa-
tion becomes more complicated as additional structures start
to appear e.g., AB
rocksalt 37 and because, at lower vol-
ume fractions
or sufficiently low
, the small spheres are
free to move from site to site within the crystal of larger
spheres. Such systems are not well suited to examination via
the lattice-switch approach, as the algorithm would have to
be modified heavily in order to cope with the free movement
of small spheres. At smaller size ratios still, it is probably
FIG. 1.
Color online The three binary structures considered in this work, with the A particles shown as pale spheres, B particles shown
as darker spheres, and the simulation cell shown as a translucent box.
a The AB13 structure consists of a simple cubic lattice of large
spheres with icosahedral clusters of 13 B particles in the body center of each cell. The 112-particle unit cell contains eight B clusters arranged
in alternating orientations, where the clusters are rotated by 90° between adjacent subcells.
b In AB2, hexagonal close-packed layers of A
particles are stacked directly on top of each other, with a honeycomb pattern of B particles places in the interstitial sites between the A layers.

c The AB
CsCl structure, consisting of two interlaced simple cubic lattices, with the B particles at the body center of the A cube.
A. N. JACKSON AND G. J. ACKLAND PHYSICAL REVIEW E 76, 066703
2007
066703-2

better to model the presence of the B particles indirectly, via
a depletion potential 38–41.
Given the range of results that exist for binary hard
spheres, the system provides an excellent testing ground for
the extension of the lattice-switch approach to multicompo-
nent systems. As our primary concern is simply to extend the
lattice-switch technique to the binary system, it makes sense
to choose a particular value of
to examine, instead of try-
ing to explore a whole range of diameter ratios before the
correctness of the approach has been determined. We have
chosen to examine the stability of AB2 and AB13 at

=0.58, as there has been a large amount of theoretical and
experimental work published for that particular diameter ra-
tio. Furthermore, due to the ambiguity concerning the stabil-
ity of the AB
CsCl structure, we shall also use the lattice-
switch approach to examine the stability of AB
CsCl
relative to the fluid phase and phase-separated solid state at

=0.73.
II. FORMULATION
We consider a system of N particles, of spatial coordinates
r
 and diameters , confined within a volume V, and sub-
ject to periodic boundary conditions. The configurational en-
ergy for this system of hard spheres has the form
E
r
 =

0 if rij  ij, ∀ i, j ,
 otherwise,

1
where rij=
ri
−rj

and ij=
1
2
i+ j. Each microstate has a
Boltzmann weight exp
−E, and so the total configurational
weight associated with this system can be written as

N,V =

i


V
dri

ij


rij − ij ,
2
where 
x1
0 for x0
0, and the product on ij
extends over all particle pairs. The associated
dimension-
less entropy density is
s
N,V 
S
N,V
Nk

1
N
ln
N,V ,
3
where S is the entropy and k is the Boltzmann constant.
This formulation describes the total entropy, integrating
over all possible positions and so over all possible phases. In
order to compare the statistical weights of the configurations
associated with individual
candidate phases, we must for-
mulate a constraint that identifies a configuration as “belong-
ing to” a given phase. To do this, we decompose the particle
position coordinates into a sum of “lattice” and “displace-
ment” vectors:
ri


= Ri
A

+ ui


,
4
where R
ARi
A
, i=1
, . . . ,N forms a set of fixed vectors
associated with the crystal structure labeled A. This set is
defined as the orthodox crystallographic lattice convolved
with the orthodox basis, but will be referred to here as sim-
ply the lattice vectors. We are thus able to initialize a simu-
lation within a particular phase, using the appropriate set of
lattice vectors
R
A and diameters
A for that phase,
just by setting u
i=0
∀ i. We can now use a simple spatial
criterion to decide whether a particular configuration belongs
to the same phase as the simulation progresses. Periodically,
throughout the simulation, the particles are checked to deter-
mine whether, after having taken any center-of-mass drift
into account, u
i is comparable to the interparticle spacing

for more details, see Sec. III C. This can happen only if the
particles have been rearranged and the structure has been
modified, and so this provides a robust criterion for ensuring
that a set of microstates belong to the same macrostate. We
note that this restriction of the configurational integral to
pure perfect phases erases contributions from equilibrium
concentrations of lattice defects, for example lattice vacan-
cies; however, for the systems in question such defects are
extremely rare and make no significant contribution to the
total free energy 42.
By performing this coordinate decomposition Eq.
4,
we have effectively changed our stochastic variables from
being r
 to u
 ,A, operating under a modified configura-
tion energy function cf. Eq.
1 of the form
E
u
,A =

0 if rij  ij, ∀ i, j ,
 otherwise,

5
where the lattice label A dictates the choice of lattice vectors
via rij
u
 ,A=

R
i
A+u
i−
R
j
A+u
j
and diameters via
ij
A= 12
i
A+ j
A
. The configuration weight associated
with a particular structure can now be written as

N,V,A =

i


A
dri



ij

rij − ij ,
6
where A signifies integration subject to the single-phase
constraint. The associated entropy density becomes
s
N,V,A 
1
N
ln
N,V,A .
7
The lattice switch itself is a Monte Carlo move that at-
tempts to switch the structure label, e.g., A→B, thus instan-
taneously exchanging one complete set of lattice sites and
diameters for another,
R
A,A → R
B,B ,
8
while keeping the displacements u
 fixed. The details of
how we construct the lattice switch and how to implement a
simulation capable of switching between two different
phases are given in Sec. III below. Once the two structures

A ,B have been explored, the entropy difference between
them can be determined as
sAB  s
N,V,A − s
N,V,B =
1
N
ln RAB
N,V ,
9
where
LATTICE-SWITCH MONTE CARLO FOR BINARY HARD-… PHYSICAL REVIEW E 76, 066703
2007
066703-3

RAB
N,V =

N,V,A

N,V,B
=
P
A
N,V
P
B
N,V
.
10
Here P
A
N ,V is the joint probability that a system free to
visit both structures will be found in the A structure.
III. IMPLEMENTATION AND METHODOLOGY
A. Defining the structures and the switch
To evaluate the free energy of a particular binary crystal,
we need to perform a lattice switch simulation capable of
visiting the target phase A and some reference phase B. As
an accurate semiempirical expression for the free energy of a
monodisperse hard-sphere fcc crystal has been determined
by Alder et al. 43,44, we choose this structure to be our
usual reference state. A lattice switch which kept the diam-
eters of the particles fixed
such that AB, would
require either a single simulation box containing coexisting
crystals of A and B particles
and thus associated interfacial
effects, making for a poor reference system or two separate
boxes corresponding to separate A and B crystals. This latter
approach can be made to work, but the complexity associ-
ated with adding a second simulation cell is unnecessary
when the two cells are being used to simulate two identical
structures that differ only in scale
by
. A simpler option is
to allow the diameters to change during the switch, turning
all B particles in the mixed phase into particles of diameter
A, and placing all these particles onto the lattice sites of a
single shared fcc lattice. This creates a switch between a
binary crystal and a single crystal, and ensures that the num-
ber of degrees of freedom in each system is the same. Fur-
thermore, as the free energy difference between fcc and hcp
structures is accurately known 8, we are free to use either
close-packed structure as a reference crystal, which affords
us a little more flexibility when it comes to choosing system
sizes. The system sizes used are shown in Table I.
Once the crystal lattices have been chosen, one must
choose a site mapping between crystals A and B. In the
previous lattice-switch work for fcc and hcp structures more
efficient switches, which preserved close-packed planes,
were used. However, for binary to single-lattice switches,
there are no such obvious similarities; therefore, we mapped
sites between crystals at random.
To predict the equilibrium phase behavior, we must first
determine the Gibbs free energy of all the candidate phases,
as a function of the pressure P
=PA
3 /kT and the composi-
tion X=NB /
NA+NB, where NA and NB are the numbers of A
and B particles, respectively. Naturally, we would seek to
measure the Gibbs free energy directly via a lattice-switch
simulation carried out in the constant pressure ensemble.
However, when attempting to switch between the pure and
binary crystals using a single simulation cell, the fact that the
equilibrium cell size is very different for these two structures
means that the switch becomes ineffective. The lattice switch
updates only the lattice and the diameters, not the cell pa-
rameters, and therefore performing the switch from the ref-
erence lattice involves attempting to transform an equilib-
rium density pure crystal into a loosely packed binary crystal

as reducing the diameter of so many particles vastly in-
creases the available volume. The conjugate move requires
a transformation from an equilibrium density binary crystal
to a pure crystal at such a high density that it is likely that all
particles overlap. To facilitate either switching move, the
volume of the simulation cell would have to be dilated so far
that the system would rapidly melt, and indeed such a switch
cannot be made to work in practice. It is possible, at least in
principle, to overcome this limitation by adding an explicit
volume dilation to the switch, so that the cell dimensions, the
lattice, and the diameters are all changed simultaneously.
However, adding a volume transformation factor affects the
probability of sampling the two phases in a nontrivial way.
For example, a dilation to a large cell will almost always be
accepted, whereas the compression to a small cell will be
difficult to achieve. This will bias the measured free energy
difference by a factor which could be determined by careful
application of the mapping transformation matrix approach
outlined in Ref. 8, but this introduces a number of compu-
TABLE I. Structures and sizes used in this work. For fcc, hcp, and AB2 structures, the unit cells were
oriented so that the close-packed layers were lying in the x-y plane. The cubic cells of AB13 and CsCl were
made commensurate with the simulation cell. The run times are shown for determining the equation of state

EOS, the weight function
WF, and the free energy difference
DF, in units of Monte Carlo sweeps

MCSs. Numbers in square brackets indicate the largest number of blocks the run was split into for block
analysis. For all structures and sizes, 20 000 MCS runs were used to equilibrate the constant pressure runs,
while the constant volume simulations needed just 5000 MCSs for equilibration.
N
Target crystal Reference crystal
Run time
MCSs
EOS
Run time
MCSs
WF
Run time
MCSs
DFStructure Dimensions Structure Dimensions
112 AB13 111 hcp 474 5105 106 40 4106 8
896 AB13 222 fcc 8814 5105 6106 240 6106 10
192 AB2 444 fcc 846 106 106 40 4106 80
648 AB2 666 fcc 899 106 5106 200 2106 200
54 CsCl 333 fcc 633 5105 106 40 4106 40
128 CsCl 444 hcp 844 5105 2106 80 5106 20
432 CsCl 666 fcc 869 5105 3106 120 8106 40
A. N. JACKSON AND G. J. ACKLAND PHYSICAL REVIEW E 76, 066703
2007
066703-4

tational overheads, and the slow switching rate mitigates
against gathering good statistics.
All these issues can be avoided in the constant volume
ensemble by keeping the density fixed when switching be-
tween the reference and target structures, thus ensuring that
the displacements are reasonable in both phases, because the
distances between particles are comparable in both phases
under these conditions.
B. Simulation methodology
To estimate the Gibbs free energy, three simulations were
carried out for each state point. The first was a single-phase
constant pressure
NPT simulation, used to estimate the
equilibrium density and cell parameters
c /a ratio for each
binary crystal, at each pressure. Then, two lattice-switch
NVT simulations were performed at the volume fractions and
c /a ratios determined from the NPT simulations. The lattice
switch attempts to transform between the target crystal
at
the target density and c /a ratio and a pure fcc crystal
at the
same density, but with c /a=1.0. The first run is used to
evaluate the multicanonical weight function and the second
uses this function to determine the free energy difference
between the crystals
see Sec. III C for more details.
Having determined the cell proportions, the density, and
the Helmholtz free energy difference between the two struc-
tures

F for a range of pressures, we can then build the
absolute free energy of the binary crystal by adding in the
free energy of the pure phase. To determine this, we integrate
the Alder equation of state 43–45,
PV
NkT
=
3
V
− 1
+ 2.566 + 0.55
V
− 1 − 1.19
V
− 12
+ 5.95
V
− 13,
11
with respect to the reduced volume V
=V /V0
where V0 is
the close-packed volume N3 /2. This yields the Helmholtz
free energy:
Fex
Alder
NkT
= − 3 ln

V
− 1
V

+ 5.124 ln V
− 20.78V
+
19.04
2
V
2
−
5.95
3
V
3 + C4 + 1,
12
where, following 45, the constant of integration C4 is set to
15.05. The Gibbs free energy can be constructed from the
Helmholtz form via
Gex = Fex
Alder +
F + PV − 1 − ln

PV
NkT
,
13
where the last two terms arise from the difference between
the Gibbs and Helmholtz free energies for an ideal gas. To
compute the phase diagrams based on these free energies, we
follow the approach of Bartlett 18 and assume that the par-
ticles are immiscible in the solid phases, and that the binary
fluid behaviour can be described accurately by the Mansoori
equation 46. At each pressure, the densities of the candidate
phases are determined from the equations of state. From
these densities, the Gibbs free energies and chemical poten-
tials of the competing phases can be determined for each
constant pressure line. The condition of equal chemical po-
tentials is then applied along the constant pressure line,
where the coexistence between the crystal phases and the
Mansoori fluid is estimated from the chemical potentials of
each species in the fluid via
see 29
A
f + nB
f
=
1 + n G
ABn ,
14
where A
f and B
f are the chemical potentials determined
from the Mansoori fluid for the A and B particles, respec-
tively, and G
ABn is the Gibbs free energy of the solid phase
that has n B particles for each A particle.
C. Multicanonical biased Monte Carlo procedure
In general, a lattice-switch move from a zero-overlap mi-
crostate in one phase is unlikely to map onto a zero-overlap
microstate in the other. Therefore, the simulation must moni-
tor not only which phase is currently being explored, but also
how close the simulation is to being able to perform the
lattice switch. To find a suitable order parameter, we first let
M
u
 ,A denote the number of overlapping pairs of par-
ticles associated with the displacements u
, within the struc-
ture A. We can then define the overlap order parameter
M
u
  M
u
,B − M
u
,A ,
15
where M is necessarily 0
0 for realizable configura-
tions of the A
B structure. To ensure that the simulation
can reach the so-called gateway states
M=0, we apply the
multicanonical Monte Carlo method 47 as a way of biasing
the simulation in a controlled fashion. To this end, we aug-
ment the system energy function such that
E
u
,A = E
u
,A + 
M
u
 ,
16
where 
M, M=0, ±1, ±2, . . ., constitute a set of multica-
nonical weights 47. Given that a suitable weight function
can be determined
see below, the canonical distribution
P
M
N ,V can be estimated from the measured multica-
nonical distribution P
M
N ,V ,  with the identification
RAB
N,V =

M0
P
M
N,V

M0
P
M
N,V
=

M0
P
M
N,V,e
M

M0
P
M
N,V,e
M
.

17
Here, the exponential reweighting maps the multicanonical
distribution onto the canonical, unfolding the bias associated
with the weights.
The multicanonical weights are determined from nonequi-
librium calculations, which are fast and give a fair approxi-
mation to the equilibrium situation. We conduct a series of
separate Monte Carlo simulations, where the particle dis-
placements
and therefore M are reset to zero at the start of
each. As these simulations will each relax to the equilibrium
M values, this forces the overall run to explore the full range
of M. Throughout this sequence of simulation, the program
monitors and accumulates the macrostate transition probabil-
LATTICE-SWITCH MONTE CARLO FOR BINARY HARD-… PHYSICAL REVIEW E 76, 066703
2007
066703-5

ity matrix, as described in Ref. 48, and at the end of the
whole run, an estimate of the weight function is determined
from the transitions probabilities by a simple shooting
method 49. This weight function is then used for the sec-
ond NVT run, during which the time spent in each phase is
recorded
with the biasing taken into account, and automati-
cally block-analyzed to produce an estimate of the free en-
ergy difference. The actual timings varied depending on the
target system and the number of particles, and are shown in
Table I.
The simulations operated under the usual Metropolis al-
gorithm 50, using random numbers generated by the
Mersenne Twister algorithm 51, as implemented in RNG-
PACK, version 1.1a 52. For each move, a sphere is chosen at
random, and a trial displacement is generated by taking the
displacement of that particle ui
and adding a random three-
dimensional vector
ui
drawn uniformly from a cube of
fixed size. The update in the displacement is accepted ac-
cording to the usual Metropolis prescription:
pa
u
 → u
 +
ui

= min„1,exp− E
u
+
ui
,A − E
u
,A… .
18
Therefore, moves that would lead to overlaps in the cur-
rent phase are always rejected, and moves that would alter
the number of overlaps in the conjugate phase are accepted
or rejected based on the multicanonical weights. The magni-
tude of
ui
was automatically adjusted to yield an accep-
tance rate of about one-third, as this has been shown previ-
ously to produce an acceptable rate of decorrelation in M
49. When M=0, the lattice-switch move is attempted once
per Monte Carlo sweep
MCS
i.e., once every N particle
moves on average, and is always accepted as there is no
energy cost or multicanonical weight difference associated
with this move.
Due to the short-range nature of the hard-sphere interac-
tion, the calculation of E
u
 ,A does not need to extend
over all N2 particle pairs. To make the simulation more effi-
cient, each particle checks only for interactions with its
neighbors, up to a maximum interaction distance of 1.5 unit
cells
i.e., 1.5 times the A-A separation. There are two lists
of neighbors for each particle, one for each phase. As men-
tioned in Sec. II, we require a criterion to determine whether
structure is stable, and the simulation achieves this by peri-
odically scanning the system for particles that have moved
further than the near neighbor distance away from their site

after having taken any center of mass diffusion into ac-
count. If this is observed to happen at any point during any
simulation at a given pressure, then the structure is consid-
ered unstable, and no longer considered as a candidate struc-
ture at that pressure.
To establish the equations of state for these systems, a
number of constant pressure simulations were also per-
formed. In this case, the simulation cell parameters
the
lengths of the sides of the cell, and therefore the volume can
fluctuate during the simulation. Such moves are accepted
with probability
pa
V → V

= min„1,exp− E
u



,A − E
u
,A
− P

V + N ln
V


/V… ,
19
where V


is the volume associated with the trial move. Of
course, this cell-size dilation can change E globally, and so
an expensive energy calculation must be carried out for each
trial volume move. Volume moves were attempted at ran-
dom, once per MCS on average. The maximum size of the
random volume dilation is automatically adjusted so that the
acceptance rate for volume moves is about one-half, which
produces an acceptable rate of exploration of the volume
parameter space 49.
IV. RESULTS
A. Pressure curves and stability
The equations of state, as measured by Monte Carlo simu-
lation in the constant pressure ensemble, are shown in Fig. 2.
The agreement with the data presented in Ref. 7 indicates
that the molecular dynamics simulation presented there used
the smaller
N=112 system size when computing the equa-
tion of state. Also, the close agreement between the Alder
equation of state 43 and the simulated fcc crystal is clear.
However, we note that while the Alder form Eq.
11 pro-
vides a reliable estimate of the equation of state, the indi-
vidual terms are difficult to reproduce. Despite having deter-
mined the density of the fcc structure at each pressure to
within at least 0.1%, attempting to fit an equation of the
Alder form to these data led to a poorly determined fit that
reproduced the original terms only to within one significant
figure. It is therefore clear that a much larger range of pres-
sures must be used to determine Alder-style equations of
state for these structures, and given that quite a modest range
of pressures was explored in the original work 44, it is not





























































0.5 0.55 0.6 0.65 0.7 0.75
φ
0
25
50
75
100
125
150
175
200
P*
Alder FCC [1964]
Eldridge AB13 [1990]
FCC N=648




AB13 N=112

AB13 N=896
AB2 N=192
AB2 N=648
CsCl N=54
CsCl N=432
FIG. 2. Equations of state determined in the constant pressure
ensemble, measuring the density as a function of the pressure. Er-
rors in the density estimates are smaller than the symbol size. The
Alder equation of state for the pure crystal 43 is shown as a solid
line, and the equation of state for AB13 as observed in Ref. 7 is
shown as a solid gray line.
A. N. JACKSON AND G. J. ACKLAND PHYSICAL REVIEW E 76, 066703
2007
066703-6

excellent agreement with previous theoretical results. We
also investigated the experimentally reported CsCl phase at
size ratio
=0.73, and found it to be unstable compared to
phase-segregated fcc or the liquid phase. However, the en-
ergy differences are small, and it is possible that either
charge effects or polydispersity may be responsible for the
experimental observation of this structure.
Only moderate system sizes and run times were required
to accurately estimate the free energies, although the lack of
technical details concerning the run times and statistical un-
certainties of the integration method results in the literature
makes a direct comparison of efficiencies impossible. Unlike
integration methods, the only errors in the lattice-switch MC
estimate of the relative free energy between the reference
and target systems come from finite-size effects and sam-
pling errors, both of which are easy to quantify. However, the
dependence upon the Alder crystal
in order to estimate the
absolute free energy is perhaps unfortunate, as it is unclear
how accurate the parameters of this equation of state are.
Furthermore, and in common with the other work in this
field, the calculation of the phase boundaries relies on the
assumption of immiscibility as embodied by Eq.
14 and
on the Mansoori equation of state for the binary fluid. The
errors associated with these last two assumptions are difficult
to estimate.
However, our results are in good agreement with the lit-
erature, including those that do not depend on the Alder for-
mulation
due to using an Einstein crystal as a reference
state 5. To understand why, we note that from Fig. 7 it is
clear that a change of the order of 1NkT is required to sig-
nificantly alter the pressure at which the transitions occur, as
the gradients of these curves are of order 1. The statistical
errors in our free energy estimates are around 0.001NkT and
the errors arising from the empirical terms of the Alder equa-
tion of state are of order 0.01NkT. Therefore, the positions of
the phase boundaries are dominated by the leading term in
the Alder equation of state and are rather insensitive to the
parameters of the other terms.
Nevertheless, removing these dependencies is a desirable
aim and an interesting topic for future research. For example,
building a lattice-switch mapping to an analytically soluble
reference crystal would allow the use of the Alder equation
of state to be avoided. This Hamiltonian-switch approach
could map from an Einstein crystal or single-occupancy cell
model to the target potential and structure, and would allow
absolute free energies to be measured directly. Also, to ex-
amine how reasonable the assumption of immiscibility is, the
lattice-switch presented here could be used to directly mea-
sure the free energy costs associated with defects and faults.
In the simplest case, one could switch between a pure fcc
crystal of A particles and a fcc crystal where one of the A
particles has been replaced with a B, and thus estimate the
free energy cost of introducing this point defect. Direct esti-
mation of the free energy of the binary fluid may also be
possible by combining the idea of using an ideal reference
state with the techniques drawn from the fluid-solid phase
switch for hard spheres 10.
Other future work includes lowering the size ratio, al-
though as mentioned above, the mobility of the smaller
spheres will present challenges. In this case it may be pos-
sible to use tethers, similar to those employed in the fluid-
solid phase-switch work 10. Also, the AB2 and AB13 struc-
tures have been observed experimentally for noble gas
mixtures 23, and this could be investigated by generalizing
our approach in the same manner as the extension to soft
potentials of the original lattice switch work 11. Given a
suitable reference structure, these soft-potential simulations
should work just as efficiently and accurately as for the bi-
nary hard-sphere system.
ACKNOWLEDGMENTS
This work was supported by EPSRC under Grants No.
GR/S10377/01 and No. GR/T11753/01. A.N.J. would like to
thank Mike Cates for helpful comments and discussions.
1 J. Keplar, Strena seu de Nive Sexangula
G. Tampach, Frank-
furt, 1611.
2 T. C. Hales and S. P. Ferguson, The Kepler Conjecture
1998
Details of Keplar’s conjecture and its solution can be found at
http://www.math.pitt.edu/~thales/kepler98/.
3 A. D. Bruce, N. B. Wilding, and G. J. Ackland, Phys. Rev.
Lett. 79, 3002
1997.
4 D. Frenkel and A. J. C. Ladd, J. Chem. Phys. 81, 3188
1984.
5 M. D. Eldridge, P. A. Madden, and D. Frenkel, Nature
Lon-
don 365, 35
1993.
6 M. Eldridge, P. Madden, and D. Frenkel, Mol. Phys. 80, 987

1993.
7 M. Eldridge, P. Madden, and D. Frenkel, Mol. Phys. 79, 105

1993.
8 A. D. Bruce, A. N. Jackson, G. J. Ackland, and N. B. Wilding,
Phys. Rev. E 61, 906
2000.
9 N. B. Wilding and A. D. Bruce, Phys. Rev. Lett. 85, 5138
2000.
10 N. Wilding, Comput. Phys. Commun. 146, 99
2002.
11 A. N. Jackson, A. D. Bruce, and G. J. Ackland, Phys. Rev. E
65, 036710
2002.
12 Jeffrey Errington, J. Chem. Phys. 120, 3130
2004.
13 G. McNeil-Watson and N. Wilding, J. Chem. Phys. 124,
064504
2006.
14 J. L. Barrat, M. Baus, and J. P. Hansen, Phys. Rev. Lett. 56,
1063
1986.
15 J. Barrat, M. Baus, and J. Hansen, J. Phys. C 20, 1413
1987.
16 J. V. Sanders, Philos. Mag. A 42, 705
1980.
17 M. J. Murray and J. V. Sanders, Philos. Mag. A 42, 721

1980.
18 P. Bartlett, R. H. Ottewill, and P. N. Pusey, J. Chem. Phys. 93,
1299
1990.
19 P. Pusey, J. Phys.: Condens. Matter 6, A29
1994.
20 M. D. Eldridge, P. A. Madden, P. N. Pusey, and P. Bartlett,
Mol. Phys. 84, 395
1995.
21 N. Hunt, R. Jardine, and P. Bartlett, Phys. Rev. E 62, 900
LATTICE-SWITCH MONTE CARLO FOR BINARY HARD-… PHYSICAL REVIEW E 76, 066703
2007
066703-9

2000.
22 A. Schofield, P. Pusey, and P. Radcliffe, Phys. Rev. E 72,
031407
2005.
23 M. Layer, A. Netsch, M. Heitz, J. Meier, and S. Hunklinger,
Phys. Rev. B 73, 184116
2006.
24 F. X. Redl, K.-S. Cho, C. B. Murray, and S. O’Brien, Nature

London 423, 968
2003.
25 C. Royall, M. Leunissen, A. Hynninen, M. Dijkstra, and A.
van Blaaderen, J. Chem. Phys. 124, 244706
2006.
26 Sayeed Abbas and Timothy P. Lodge, Phys. Rev. Lett. 97,
097803
2006.
27 P. Bartlett, J. Phys.: Condens. Matter 2, 4979
1990.
28 X. Cottin and P. A. Monson, Int. J. Thermophys. 16, 733

1995.
29 X. Cottin and P. A. Monson, J. Chem. Phys. 102, 3354
1995.
30 H. Xu, J. Phys.: Condens. Matter 4 L663
1992.
31 A. Denton, Phys. Rev. A 42, 7312
1990.
32 S. Smithline, J. Chem. Phys. 86, 6486
1987.
33 A. Schofield, Phys. Rev. E 64, 051403
2001.
34 M. Leunissen, C. Christova, A. Hynninen, C. Royall, A.
Campbell, A. Imhof, M. Dijkstra, R. van Roij, and A. van
Blaaderen, Nature
London 437, 235
2005.
35 P. Bartlett and A. Campbell, Phys. Rev. Lett. 95, 128302

2005.
36 A. Hynninen, M. Leunissen, A. van Blaaderen, and M. Dijk-
stra, Phys. Rev. Lett. 96, 018303
2006.
37 E. Trizac, M. Eldridge, and P. Madden, Mol. Phys. 90, 675

1997.
38 Ph. Germain and S. Amokrane, Phys. Rev. E 65, 031109

2002.
39 Marjolein Dijkstra, René van Roij, and Robert Evans, Phys.
Rev. E 59, 5744
1999.
40 A. Ayadim and S. Amokrane,Phys. Rev. E 74, 021106
2006.
41 E. Velasco, G. Navascues, and L. Mederos, Phys. Rev. E 60,
3158
1999.
42 Richard Bowles and Robin Speedy, Mol. Phys. 83, 113

1994.
43 B. Alder, J. Chem. Phys. 40, 2724
1964.
44 B. J. Alder, W. G. Hoover, and D. A. Young, J. Chem. Phys.
49, 3688
1968.
45 D. Young, J. Chem. Phys. 70, 473
1979.
46 G. Mansoori, J. Chem. Phys. 54, 1523
1971.
47 B. A. Berg, Fields Inst. Commun. 26, 1
1999.
48 G. R. Smith and A. D. Bruce, Phys. Rev. E 53, 6530
1996.
49 A. N. Jackson, Ph.D. thesis University of Edinburgh 2001.
50 N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. H. Teller,
and E. Teller, J. Chem. Phys. 21, 1087
1953.
51 Makoto Matsumoto and Takuji Nishimura, ACM Trans.
Model. Comput. Simul. 8, 3,
1998.
52 Paul Houle, computer code RNGPACK 1.1A, http://
www.honeylocust.com/rngpack/
A. N. JACKSON AND G. J. ACKLAND PHYSICAL REVIEW E 76, 066703
2007
066703-10