Computer simulation study of free energy barriers in crystal nucleation
J. S. van Duijneveldt
Vun ‘t ~ioj~Luborumy, UttiwrsirJ?
of
litrecht, Pudttaluun 8, 3.584 CH Utrecht, The Netherlands
D. Frenkel
FOM Institutejbr .4 tontic and M&x&w Phy,Gx, Kruisluun 407, 1098 SJSJ,Qmsrerdam, The Nethulands
(Rcccivcd 2 July 199 1; accepted 29 November 199 1)
We show how relatively standard Monte Carlo techniques can be used to probe the free-energy
barrier that separates the crystalline phase from the supercooled liquid. As an ilIustration, we
apply our approach to a system of soft, repulsive spheres [ LI( Y) = E(C/Y) 12]. This system is
known to have a stable Face-centered-cubic (fee) crystal structure up to the melting
temprrat.ure. However, in our simulations, we find that there is a surprisingly low free-energy
barrier for the formation of body-centered-cubic (bee) crystallites from the melt. In contrast,
there appears to be no ‘easy’ path from the melt to the (stable) fee phase. These observations
shed new light on the results of previous simulations that studied the dynamics of crystal
nucleation in the r- ” 1 system. We argue that the techniques developed in this paper can be
used to gain insight in the process of homogeneous nucleation under conditions where direct,
dynamical simulations are inconclusive or prohibitively expensive.
I. INTRODUCTION
Computer simulation studies have contributed consid-
erably to our understanding of the freezing transition. In
1357, Alder and Wainwright’ and Wood and Jacobson”
showed that a system of purely repulsive hard spheres can
undergo a first-order transition from the fluid to the crystal-
line state.
Some
10 years later, Hoover and Ree showed how
computer simulations can be used to locate the freezing
point of an arbitrary ‘atomic’ fluid.” Subsequently, simula-
tions have been used to determine both the freezing point
and the structure of the solid at freezing for a wide range of
interatomic potentials (see, e.g., Refs. 66). In particular,
Hoover et aZ.? investigated the stability of different crystal-
line phases for a class of model systems with repulsive pair
interactions of the form c(r) = e(rr/r)“. This study showed
that for ~27, the stable crystal structure at melting is face
centered cubic (fee), whereas at lower values of n, the solid
melts from the body-centered-cubic (bee) phase. In other
words, ‘hard’ repulsive interactions favor fee solids, while
bee occurs for ‘soft’ interactions.
In the mid-scventics, computer simulations were first.
used to study the ctyrzarCcsofcrystallization, rather than the
thermodynantics. Starting with the work of Mandell et al.,’
a large number of numerical simulations of the homoge-
neous nucleation process have been rcported.d.“-‘” For a re-
view, see! e.g., Refs. 25 and 26.
Computer studies of the dynamics of homogeneous
crystallization arc much less straightforward than the nu-
merical determination of the fluid-solid coexistence point.
R&w, we list some of the factors
that
make computer-simu-
lations studies of crystallization time consuming and, in
many casts, difficult to intcrprct.
Let us first briefly recall the usual approach to study
crystallization in a computer ‘experiment. In such a molec-
ular dynamics (MD 1 Qmulation, the system is first prepared
in a stable, well-equilibrated fluid phase. The temperature of
this fluid is then rapidly quenched to a value well below the
freezing point ( Tf) (often as low as
0.5T,,
or even lower).
The system is then allowed to evolve. During this time evolu-
tion one monitors a number of dynamical observables (e.g.,
the pressure, or the structure factor) that can be used to
detect the onset of homogeneous nucleation.
This procedure has a number of practical drawbacks.
First of all, the system size appears to have a large etfect on
the rate of nucleation. The size of the critical nucleus and the
time of formation of a critical nucleus are smaller in a small
system than in a large system.‘3,27 Recently, Swope and An-
dersen have reported a simulation study of homogeneous
nucleation in a very large system ( 10’ particles). This study
indicates that system-size effects become unimportant only
for model systems containing more than IO3 particles. Such a
system size is a factor lo-50 larger than what is typically
needed to study equilibrium properties of dense fluids. As a
consequence, nucleation studies tend to be computationally
expensive.
This problem is compounded by the fact that it may take
a long time (and therefore a lot of computer time) before a
nucleation event is observed, unless the fluid is strongly su-
percooled. How much supercooling is needed to achieve a
reasonable rate of nucleation (on the time scale of a comput-
er ‘experiment’), depends strongly on the harshness of the
interatomic potential. Mountain et nl.‘” have performed
simulation studies of crystal nucleation in two soft-sphere
systems [L)(T) .-Y-“, with n = 7 and II = 121 and in a Len-
nard-Jones system. They observed that nucleation and
growth take place more easily as the interaction between
particles becomes softer. A similar conclusion was reached
by Robbins et ~1.” who studied crystal nucleation in a
Yukawa fluid. The strong dependence of the nucleation rate
on the harshness of the short-range repulsion between atoms
can be understood intuitively by noting that soft spheres will
slip more easily into a new position than densely packed hard
Jo Chem. Phys. 96 (6),15 March 1992
0021-9606/92/ 064655- 14$06.00
0 1992 American Institute of Physics
4655
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
4656
J. S. van Duijneveldt and D, Frenkel: Free energy barriers to nucleation
spheres. Experimentally, the same correlation between nu-
cleation rate and the softness of the intermolecular potential
has been observed in a study of the formation of colloidal
crystals.”
In order to see homogeneous crystal nucleation in a
computer simulation, the initial fluid phase must be super-
cooled to a temperature that is far below the equilibrium
freezing temperature. This is true for fluids with a ‘soft’ in-
tcrtnolecular potential, but even more so for the more harsh-
ly repulsive systems (i.e., those systems that tend to form
stable fee crystals). Unfortunately, it is not known a priori if
crystal nucleation at such extreme supercooling proceeds in
the same way as crystallization close to the freezing tempera-
ture. This fact severely limits the use of MD studies to inves-
tigate whether nucleation proceeds directly into the thermo-
dynamically stable phase or that a metastable solid phase is
formed first, in accordance with Ostwald’s ‘step rule’.‘” The
latter rule states that the crystal structure that is nucleated
from the melt does not correspond to the most stable phase,
but to the phase that is nearest (in free energy
j
to the fluid
phase. A closely related rule is the one due Stranski and
Totomanow,“” who argue that the crystal phase that forms
first is the one with the lowest free energy barrier for nuclea-
tion. One of the aims of the early work of Mandell
et
al.‘? was
precisely to test whether Ostwald’s rule applied to the homo-
geneous crystallization of the Lennard-Jones fluid. In Ref. 9
evidence was presented that indicated the formation of bee
nuclei in a small, strongly supercooled Lennard-Jones (LJ)
system. However, all other studies of the same system found
evidence for the formation of fee nuclei,‘0,8.‘4 In particular,
the recent numerical study by Swope and Andersen16 of
crystal nucleation in a system of 10’ LJ particles, shows con-
vincingly that, although both bee and fee nuclei can form
under conditions of strong supercooling, only the fee nuclei
(i.e., the
ones
that correspond to the thermodynamically sta-
ble phase) grow to form larger crystal&s. In contrast, nu-
merical studies of homogeneous nucleation in (considerably
smaller) r
‘“-‘soft-sphere’ systems’831” indicate that, al-
though fee is the thermodynamically stable crystal phase,
bee crystallites tend to form even at large supercooling. This
observation is in qualitative agreement with the theoretical
predictions of -4lexander and McTague”’ who argue that, at
least for small supercooling, nucleation of bee crystallites
should be favored in all simple fluids that have a ‘weak’ first-
order freezing transition. A recent theoretical study by Klein
and Leyvraz”’ supports the idea that a metastable bee phase
should nucleate relatively easily from the melt. Experimen-
tally, nucleation of a met&able bee phase has been observed
in rapidly cooled metal melts.“3
nucleation can therefore be considered as the product of two
terms, namely, ( 1) the probability to find the system at the
top of the free-energy barrier to nucleation and (2) the rate
at which this ‘activated’ state transforms into a stable crys-
tallite. Below, we show how relatively simple computer-sim-
ulation techniques can be used to gain information about the
free-energy barrier to nucleation. We shall not discuss the
second factor, although we stress that the techniques de-
scribed in the present paper could be used to prepare a slight-
ly supercooled fluid at the top of the free-energy barrier. The
standard techniques to compute chemical reaction con-
stants33.“5
could then be used to compute the actual rate of
nucleation.
The central problem in determining the nucleation bar-
rier is that the ‘reaction coordinate’ from fluid to crystal is
not known apriori. Moreover, we do not wish to impose any
specific reaction path (e.g., one leading from the fluid to
only one of all the possible crystal phases). Rather, we wish
to use as our reaction coordinate any order parameter 9 that
is sensitive to the overall degree of crystallinity in the system
but much less sensitive to the differences between the possi-
ble crystal structures. Below we show that it is indeed possi-
ble to define such an ‘unbiased’ reaction coordinate.
The Helmholtz free energy of the system, F, is a function
of this order parameter. In fact, it follows from equilibrium
statistical mechanics3” that
f;‘(a) = constant - k,Tln{P(@)),
(1)
where P(Q) is the probability per unit interval to fmd the
order parameter around a given value of <p. We expect that
above the freezing temperature, P(a) will be strongly
peaked at a low (liquidlike) value of a. Close to coexistence,
P(Q) should develop a double-peaked structure. The second
peak corresponds to a solidlike value of 9. In bet.ween these
two peaks, P( Q,) will be very small, but finite. From Eq. ( 1)
we obtain our estimate for the barrier to nucleation from the
minimum value of P(Q) between the solidlike and the li-
quidlike peak.
Clearly, it would be desirable to have a numerical tech-
nique that allows us to study crystal nucleation at tempera-
tures closer to the coexistence point. The straightforward
MD approach will not work under those conditions because
the rat.e of nucleation depends exponentially on the degree of
supercooling. A mildly supercooled fluid will never crystal-
lize during a molecular dynamics simulation. In the present
paper we do not attempt to solve the full problem of crystal
nucleation close to coexistence. Rather, we exploit the fact
that nucleation is an activated process and that the rate of
As P( a,) is an equilibrium property of the system under
consideration, it can, in principle, be probed both by Monte
Carlo and by molecular dynamics simulations. However, in
an ordinary simulation, the system will spend almost all its
time in either theliquid or the solid state and it will therefore
be impossible to obtain good statistics on P(a) for interme-
diate values of the order parameter. We circumvent this
problem by using a non=Boltzmann Monte Carlo sampling
scheme, namely, the ‘umbrella sampling’ t.echnique of Torrie
and ValleauX7 Using this technique, we are able to measure
F(Q) as a function of the crystallinity G. Moreover, we ob-
tain information about the structure of the solid phase at the
other side of the barrier. Below we show that., at least for the
r
“system, the solid phase that is reached from the fluid by
crossing the lowest free-energy barrier is, in fact, not the
equilibrium fee phase but the bee phase.
The rest of this paper is organized as follows. In Sec. II
we briefly summarize the umbrella sampling technique. In
Sec. III, we discuss our choice for t.he crystallinity order
parameter. The computational details are summarized in
Sec. IV. Section V contains our results for the nucleation
J. Chem. Phys., Vol. 96, No. 6,15 March 1992
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
barriers in the I’m ” system and in Sec. VI we discuss some of
the implications of the present work.
II. UMBRELLA SAMPLING
As indicated in the previous section, special numerical
techniques are required to sample the equilibrium distribu-
tion function P(9). There exist two closely related tech-
niques to measure such distribution functions. One is based
on a dynamical determination of the force conjugate to the
‘reaction coordinate’ +.” Such a calculation requires a se-
ries of MD simulations at different fixed values of @. In our
case, where the order parameter happens to be defined as a
discontinuous function of the atomic coordinates (see Sec.
III), this ‘constraint’ method is less attractive. We have
therefore used the umbrella-sampling t.echnique of Torrie
and Valleau-‘7 to measure P(G). The basic idea behind this
method is to bias the Monte Carlo sampling of configuration
space in such a way that points high on the free-energy bar-
rier are sampled frequently. In order to compute P( ‘P), we
must correct for our sampling bias. Below, we briefly sum-
marize how this is achieved in practice.
Consider a Monte Carla simulation in the canonical
(&IY) ensemble. The average of a quantity A which de-
pends on the atomic coordinates qN is given by
(-4 )rvpT I
s
&Y4(qN)exp( -/W‘(q”))
I
@exp( -/3T(q”)),
(2)
where ?“(q”) is the potential energy function and
/&G I/k, T. The umbrella sampling scheme was devised to
handle situations where important contributions to (A )
come from configurations for which the Boltzmann factor
exp ( - pa “) is small, in which case Eq. (2) yields poor
s.tatistics. The umbrella sampling scheme is based on the fact
that Eq. (2) can be rewritten as
(*4 jNkmr =
s
dqNur ‘(@)A(@)exp( -BW’(qN))dqN)
a
I
&f~ut- *(q”)eap( -fl~Y‘(q”))~(q”),
(3)
where 10 is an, as yet unspecified weight function. Rather
than perform a Monte Carlo sampling of the original (Boltz-
maml) distribution function, we now consider configura-
tional averages obtained by sampling according to the biased
distribution function exp( - fl7’)u-r. We Iabel averages in
the original and biased ‘ensembles’ with the subscripts 0 and
w, respectively,
(,4 ),, = ‘:” /““c’
= (A /w),(w)o.
The right-hand side of Eq. (4) shows that, in order to obtain
good statistics for {A jo,
the biased distribution function
exp ( - p”z ‘O) ICI should overlap with A /w, while w itself must
overlap with esp( - PS’.). This ‘bridging’ nature of w ex-
plains the name ‘umbrella ssmpling.’
(4)
the first window to estimate w for the next (partially over-
lapping) window, and so on. Actually, for our purpose it is
not essential that consecutive windows actually overlap. We
wish to determine the variation of F‘(a) with Q, [see Eq.
( 1) 1. If both F( @) and its derivative are smoothly varying
functions of @, we can proceed as follows. We measure
an@) _
_ k T&In p(a)
aa
H
aa
Finally, we have to decide on a suitable form for U. In
general, this choice will depend on the region in configura-
tion space that contributes most to the desired average (A ).
In the present case (A > corresponds to the probability distri-
bution function P( Cp) of some order parameter @:
P(a) = (a(@ - Qi(q”))).
(5)
We are interested in the free-energy barrier F(Ca> [see Eq.
( 1) 1. The minima of this free energy occur at values of Q,
that correspond to homogeneous phases. When the system
moves from one phase to another, Fwill go through a maxi-
mum. Clearly, we should construct a biasing function uy that
allows us to probe P( @) over the entire free-energy barrier.
The
optimum choice
for w
would be
w = exp( + fiF( a>) ), because in that case the biased distri-
bution function would be flat, i.e., all values of 9 are sampled
with the same frequency. However, this choice of w requires
the final answer of the calculation as input. In practice, we
first perform a simulation without any weighting function.
This yields a local estimate for F( Q,). Next, the simulation is
repeated using the current estimate for F(a) to construct
the biasing function wzexp(@). At first sight, it might
seem advantageous to refine the computation of w in such a
way that allrelevant Q values can be sampled in one run. We
have not attempted to do so. Rather, we have performed a
number of umbrella sampling in (partially overlapping) <P
kindows.’ The reason for doing so is the following. Let us
assume that we sample an interval (P,,, - (P,,i,, =A@ in n
umbrella-sampling simulations. The optimum choice of n is
clearly the one that samples the complete @ interval in the
minimum computing time. In order to estimate this time, let
us assume that the system performs a random walk in Cp
space within the window A+/rt. Associated with the ran-
dom walk in Cp space is a ‘diffusion constant’ D,, . The char-
acteristic time needed to sample one interval A+/n is
(ha/it) ’
7, =
Da .
Clearly, the total time to sample all n windows is
The important point to note is that the computing time de-
creases with increasing n. It would, however, be incorrect to
assume that n should be chosen as large as possible. The
actual equilibration time of a run in one of the d> windows
also depends on the rate at which all coordinates ‘orthogo-
nal’ to Q, are sampled. Let us denote this time by TV. Clearly,
once rl becomes appreciably larger than T,,, the total com-
putation will scale as n >( TV. This suggests that the optimum
choice of n is the one for which ‘I;, ~=7~.
In
our simulations, we used the results of a simulation of
in a finite number of (nonoverlapping) windows. We then fit
J. S. van Duijneveldt and D. Frenkei: Free energy barriers to nucleation
4657
J. Chem. Phys., Vol. 96, No. 6, 15 March 1992
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
an analytic function (e.g., a polynomial in Cp) to these values
of the free-energy derivative. Once we have an analytic
expression for [ c?~(@)/&@], we can integrate it to compute,
for example, the height of the free-energy barrier. Although
this version of umbrella sampling is not common, its molecu-
lar-dynamics equivalent is: In a ‘constrained’ MD simula-
tion at fixed Qt, one measures the conjugate force
[JF(S)/a(a] and uses this information to reconstruct
F(Q).
For simulations on small systems, there is no great ad-
vantage in the use of nonoverlapping windows. However, for
large systems, it can be quite expensive to work with overlup-
ping windows. In contrast, the nonoverlapping window ap-
proach will work equally well for small and for large sys-
tems.
III. ORDER PARAMETERS
In order to compute the free-energy barrier that sepa-
rates the isotropic liquid from the crystalline state, we must
first specify an order parameter that measures the degree of
crystallinity in the system. It is important that this order
parameter be chosen such that it does not favor one crystal
structure over all others. This implies that our order param-
eter should be large for ali crystalline phases and small for
the isotropic fluid. Moreover, the value of the order param-
eter should be insensitive to the orientation of the crystal in
space. And finally, it should be relatively easy to compute.
As we shall see below, these requirements severely limit the
number of useful order parameters.
A three-dimensional crystal is characterized by two dis-
tinct types of order that do not exist in an isotropic fluid,
namely, translational order and orientational order. Trans-
lational order is commonly considered to be most character-
istic feature of a crystal: The particle positions are repeated
periodically in three independent directions. A quantitative
measure of this translational order is the magnitude of the
structure factor, S(k), where the wave vector k is a basis
vector (Bragg vector) in the reciprocal lattice correspond-
ing to the crystal lattice under consideration. In simulations
of a finite, periodic system, k should also be commensurate
with the periodic boundary conditions of the system. It is
easy to see that S(k) does
not
satisfy the criteria for a ‘good
order parameter that we have specified above. First of all,
different lattices have digerent Bragg vectors. Hence, a
choice of k that leads to a large value of S( k) for one kind of
lattice, may yield a small S(k) for another lattice. In addi-
tion, the fact that k must be commensurate with the periodic
box makes that a crystallite that is aligned with the periodic
box, will appear to have a much larger degree of crystallinity
than an identical crystallite that happens to be rotated with
respect to the box edges.
There are, however, other measures of the order in a
system that are rotationally invariant. For example, one can
obtain detailed information about the locai order in a many-
body system by analyzing the Voronoi tessellation of space
(see, e.g., Ref. 18). The Voronoi polyhedron associated with
a given particle is defined as the set of points that are closer to
that particle than to any other particle in the system. Clearly,
the Voronoi polyhedra are convex and fill all space. One way
4658
J. S. van Duijneveldt and D. Frenkel: Free energy barriers to nucleation
to characterize a configuration of N atoms, is to count the
number of triangular, square, pentagonal, etc. faces of the
Voronoi polyhedron. It is customary to define the signature
of a Voronoi polyhedron as a string of numbers that denote,
successively, the number of triangular, square, pentagonal,
etc. faces of the polyhedron. For example, the signature of a
Voronoi polyhedron around a particle in a perfect bee envi-
ronment is (0608): its surface consists of 0 triangles, 6
squares, 0 pentagons, and 8 hexagons. In general, the instan-
taneous surroundings of a particle in a fluid or even solid will
differ from those in a perfect crystal lattice. As a conse-
quence, a given structure will be characterized by a distribu-
tion of Voronoi signatures, rather than a single one. Still, it is
possible to use the Voronoi analysis to distinguish different
crystal structures. However, the Voronoi signature does not
provide a convenient orderparameter to measure ‘crystallin-
ity.’ The reason is that different crystal structures have very
different Voronoi signatures, whereas, in the present study,
we need an order parameter that is large for all crystal struc-
tures, while it should vanish (at least in the thermodynamic
limit) for a disordered phase.
Another, rotationally invariant, measure of the crystal-
linity of a 30 system is provided by the bond order of the
system. Orientational order, or bond order, is always present
in three-dimensional crystals: The positional ordering of the
particles implies that any particle is surrounded by other
particles in certain preferred directions. As we shall see be-
low, there exist quantitative measures of this bond order that
are large for all crystal structures of interest and that are,
moreover, invariant under rotation of the crystallite. For
this reason, we shall use the degree of bond order in the
system under study to quantify the degree of crystallinity in
our simulations. As we shall show below, bond order is rea-
sonably easy to calculate, and permits a neat differentiation
to be made between isotropic and ordered states.
The bond-order parameters that we have used were in-
troduced by Steinhardt et ~1.‘~ and have been used by other
authors to analyze the structure of a nucleating system.20*‘”
The analysis starts with a definition of the neighbors of a
particle (for instance, using the Voronoi construction). The
vectors rjoining neighbors are called bonds. Associated with
every bond is a set of numbers
CL (r) = Yfm (Q(r),&-) 1,
(6)
where Y, CS,$? are spherical harmonics and 8(r) and #(r)
are the polar and azimuthal angles of vector r with respect to
an arbitrary reference frame. Only even-l spherical harmon-
ics are considered, because the permutation of a pair of
(identical) particles ought not affect the bond-order param-
eter. The Q,,,, ( r ) defined in Eq. ( 6 ) are still local order pa-
rameters. We now define a global orientational order param-
eter Qlm as
where the sum runs over all ni, bonds in the system. These
e, still depend on the choice of reference frame. However,
from the &,
rotationally invariant combinations can be
constructed:
?%,, -$x
Q!,,, 01,
b
J. Chem. Phys., Vol. 96, No. 6,15 March 1992
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
W[I 2
1>1,,l,l>,n13
(
I 1
1
ml 4
m3 *Qbn,Qlm,Dlm,. (9)
>
811) * m> “5. “$
0
Q, and jV, are called second-order and third-order invar-
iants, respectively. The term in brackets in Eq, (9) is a
Wigner-3jsymbol (see, e.g., Kef. 39). The order parameters
Q, and P-ti, still depend on the definition of neighbors that is
used. In general, this is not a particularly serious problem. In
“2%~ event, it is possible to define a reduced order parmeter
IV
thdt
is hardly sensitive to the precise definition of the
nearest neighbors of‘s particle:
(10)
6; zs Wj(z
lQbn ,,)?
These g7! provide a convenient quantitative measure for the
prevalent crystal st.ructure of the system under study. In Ta-
ble I, some low-Z order parameters are given for simple clus-
ter geometries. Of these only the icosahedral geometry can-
not be repeated periodically to fill all space. For symmetry
re<asons the first nonzero values occur for I = 4 in clusters
with cubic symmetry and for I = 6 in clusters with icosahe-
dral symmetry. Inspection of Table I shows that knowledge
of all four order parameters listed allows us to distinguish the
most common structural units of simple atotnic systems. All
order parameters listed in Table I vanish in the isotropic
fluid phase, at least in the thermodynamic limit. As can be
seen from Table I, the order parameter Q6 is of the same
order of magnitude for all crystal structures of interest. This
makes rZ, less useful to distinguish different crystal struc-
tures, but very useful to act as a generic measure of crystal-
linity. For this mason, we have selected Q6 to play the role of
& crystalline order parameter in all our simulations. The
‘reaction coordinate’ from isotropic liquid to crystal there-
fore corresponds to a path of increasing Q6. The other order
parameters that were discussed in the present section have
been used to analyze the nature of the crystalline structures
formed in our simulations.
IV. COMPUTER SIMULATIONS
A. Introduction
As an illustration of the techniques described in the pre-
vious sections, we have performed a series of computer simu-
TABLE 1. I3o~d oricntationnl order parameters [see Eqs. (8) and (1011 for
a number of simple cluster geometries. fee: face-centered-cubic structure,
hcp: hrxagomil close-packed structure, hcc: body-centered-cubic structure,
and SC: simple cubic structure.
Geometry
>
Q4
Q,
%;
Icasahedt-al
0
0.663 32
0
-
0.169 754
fee
0.190
94
0.574 52
~= 0.159
317 -
0.013 161
hcp
0,OY 7
22
0.484 76
0.134
097 -
0,012 442
bee
0.036
37
0.510 69
0.159
317
0.013 161
SC
0.763
76
0.353 55
0.159
317
0.013 161
C, liquid)
0
0
0
0
J. S. van Duijneveldt and D. Frenkel: Free energy barriers to nucleation
4659
lations to measure the variation of the free energy of a simple
atomic system with the ‘crystallinity’ order parameter Q6.
We selected a model about which much is known from other
simulations, namely, a system of particles interacting
through an r “-soft-sphere potential,
u(r) = E(c7/r)‘2.
(11)
The model parameters E and (T can be used to construct di-
mensionless expressions for the density (p* =pd), the tem-
perature ( T * - k, T/E) and the pressure (P * = Pd./e). In
fact, the thermodynamic state of this model system can be
completely specified by a single dimensionless parameter>
namely, the reduced density
/jqg
(d-/k, T)
i/4 =p*(T*) l/4.
(12)
However, in order to facilitate comparison with other model
systems we prefer to specify T * and P * separately.
In view of the exploratory nature of the present work,
most calculations were performed on quite small periodic
systems containing 125 or 128 particles. The former system
size is a ‘magic number’ for a bee crystal, while the latter
number is optimal for an fee crystal. We cannot hope to
extract quantitative estimates on nucleation barriers from
simulations of such small systems. However, as simulations
of such small systems are relatively cheap, we can easily
study a number of different supercoolings and different ini-
tial conditions. This allows us to discuss most ofthe practical
aspects of the simulation technique. In our simulations of
these small systems, we found that the free-energy barrier
separating the fluid from the (stable) fee phase was, invaria-
bly, much higher than the barrier between the fluid and a bee
crystal. In order to test whether this trend is an artifact of the
small system size, we performed one simulation for a larger
system (N= 1000) at 20% supercooling, This simulation
confirmed that it is relatively easy to go from the fluid to the
bee phase, yet it also showed that the magnitude of the free-
energy barrier is quite sensitive to the system size. However,
the simulation of the lOOO-particle system also demonstrates
that the techniques that we employ are not limited to small
systems. Hence, the fact that relatively large systems must be
studied in order to arrive at a quantitative estimate of the
nucleation barrier does not constitute a serious problem.
B. Computational details
All simulations were carried out at constant pressure,
using the Mont.e Carlo sc.heme described in Ref. 40. The
reason to carry out simulations at constant pressure, rather
than at constant volume is that former condition corre-
sponds more closely to the experimental situation. This is of
particular importance near coexistence. Using a constant
pressure simulation it is possible, at least in principle, to
transform one phase, at constant temperature and pressure,
into the other phase with which it coexists. This is not possi-
ble in constant NVTsimulations. If we measure the distribu-
tion function P(Q, ) at constant pressure, rather than at con-
stant volume, we obtain information about the Gibbs
free-energy barrier between the solid and the liquid phase.
Eq. ( 1) is consequently replaced by
J. Chem. Phys., Vol. 96, No. 6,15 March 1992
Downloaded 11 Oct 2004 to 145.18.129.130. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
