Transcript Determining the Crystal-Melt Interfacial via Molecular

Determining the Crystal-Melt Interfacial Free
Energy via Molecular-Dynamics Simulation
Brian B. Laird
Department of Chemistry
University of Kansas
Lawrence, KS 66045, USA
Ruslan L. Davidchack
Department of Mathematics
University of Leicester
Leicester LE1 7RH, UK
Funding: NSF, KCASC
Prestissimo Workshop 2004
Observation #1
Not all important problems addressed with
MD simulation are biological.
In this work we describe application of
molecular simulation to an important problem
in materials science/metallurgy
Molecular
Interactions
Material
Properties
Problem: Can one calculate the free energy of a
crystal-melt interface using MD simulation?
Crystal-melt interfacial free energy, gcl
– The work required to form a unit area of interface
between a crystal and its own melt
Solid
Liquid
Why is the Interfacial Free Energy Important?
gcl is a primary controlling parameter governing the
kinetics and morphology of crystal growth from the melt
Example: Dendritic Growth
Experiment
Model
The nature of dendritic growth from the melt
is highly sensitive to the orientation
dependence (anisotropy) in gcl.
Data from simulation can be used in
continuum models of dendrite growth
kinetics (Phase-field modeling)
Growth of dendrites in succinonitrile
(NASA microgravity program)
Why is the Interfacial Free Energy Important?
gcl is a primary controlling parameter governing the
kinetics and morphology of crystal growth from the melt
Example: Crystal nucleation
The rate of homogeneous crystal nucleation
from under-cooled melts is highly dependent
on the magnitude of gcl
knucleation  AeBg cl / kT
Nucleation often occurs not to the
thermodynamically most stable bulk crystal
phase, but to the one with the lowest kinetic
barrier (i.e. lowest gcl) - Ostwald step rule
Observation of nucleation in colloidal
crystals: Weitz group (Harvard)
Why are simulations necessary here?
Direct experimental determination of gcl is difficult and
few measurements exist
Direct measurements: (contact angle)
•Only a handful of materials: water,
cadmium, bismuth, pivalic acid,
succinonitrile
Indirect measurements: (nucleation)
• Primary source of data for g.Accurate only
to 10-20%
Simulations needed to determine phenomenology
Indirect experimental measurement of gcl
(Nucleation data and Turnbull’s rule)
gcl can be estimated from nucleation
rates: typically accurate to 10-20%
Turnbull (1950) estimated gcl from
nucleation rate data and found the
following empirical rule:
gcl r2/3 CT DfusH
where the Turnbull Coefficient CT is
~.45 for metals and ~ 0.32 for nonmetals
Can we understand the molecular
origin of Turnbull’s Rule???
Calculating Free Energies via simulation:
Why is Free Energy hard to measure?
•
Unlike energy, entropy (& free
Energy) is not the average of some
mechanical variable, but is a function
of the entire trajectory (or phase
space)
 (H({ p,q}) d
E  E({ p,q}); S  k ln
•

Free energy, F = E - TS, calculation
generally require a series of
simulations slowly transforming the
system from a reference state to the
state of interest
H( p,q)  H 0 ( p,q)  H1 ( p,q)
F( H)  F( H 0 ) 

1
0
H1 d
Thermodynamic Integration
Frenkel Analogy
Energy  Depth of lake
Entropy  Volume of Lake
Observation #2
In molecular simulation there are almost always
two (or more) very different methods for
calculating any given quantity
Calculating a quantity of interest using more than one method is an
important check on the efficacy of our methods
Cleaving
Method
Fluctuation
Method
gcl
Fluctuation Method
•
•
•
•
Method due to Hoyt, Asta & Karma, PRL 86
5530, (2000)
h(x) = height of an interface in a quasi twodimensional slab
If q = angle between the average interface
normal and its instantaneous value, then the
stiffness of the interface can be defined
d 2g
gˆ  g  2
dq
The stiffness can be found from the
fluctuations in h(x)

•
h(x)
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
kB T
 h(k) 
ˆ k2
bW g
2
Advantage:
• Anisotropy in g precisely measured
•

Disadvantages:
• Large systems (N = 40,000 - 100,000)
• Low precision in g itself
W
b
Cleaving methods for calculating interfacial free energy
•Cleaving Potentials: Broughton & Gilmer (1986)
- Lennard-Jones system
•Cleaving Walls: Davidchack & Laird (2000) - HS, LJ, Inverse-Power potentials
•Advantages: very precise for g, relatively small sizes (N = 7000 - 10,000)
•Disadvantages: Anisotropy measurement not as precise as in fluctuation method
crystal
cry + stal
Calculate g directly by
cleaving and rearranging bulk phases to give
interface of interest
THERMODYNAMIC
INTEGRATION
liquid
cry
uid
liq + uid
cry uid
How is the cleaving done?
We employ special cleaving walls made of particles
similar to those in the system •Choose a dividing surface – particles
to the left of the surface are labeled
“1”, those to the right are labeled “2”
•Wall labeled “1” interacts only with
particles of type “1”, same for “2”
•As walls “1” and “2” are moved in
opposite directions toward one
another the two halves of the system
are separated
•If separation done slowly enough the
cleaving is reversible
•Work/Area to cleave measured be
integrating the pressure on the walls
as a function of wall position.
Observation #3
Physical reality is overrated in molecular
simulation
In calculating free energies via simulation, we
only care that the initial and final states are
“physical”, we can do (almost) anything we want
in between
Cleaving methods for calculating interfacial free energy
Solid
Liquid
Step 1
Step 2
A B
B A
Step 3
A A
Start with separate solid and liquid
systems equilibrated at coexistence
conditions: T, rc, rf
B B
Step 4
Steps 1, 2: Insert suitably chosen
“cleaving” potentials into the solid and
liquid systems
Step 3: Juxtapose the solid and liquid
systems by rearranging the boundary
conditions while maintaining the
cleaving potentials
Step 4: Remove the cleaving potentials
from the combined system
gAw1 + w4 + w3 + w4
Systematic error: hysteresis
The main source of uncertainty in the obtained results is the presence of
a hysteresis loop at the fluid ordering transition in Step 2
Reducing Hysteresis
•longer runs
•improve cleaving wall
design
•cleave fluid at lower
density (adds an extra
step)
Our approach: a systematic study of the
effect of inter-atomic potential on gcl
• Simplest potential - Hard spheres
• Effect of Attraction - Lennard Jones
• Effect of range of repulsive potential
- inverse power potentials
First Study: The Simplest System
Hard Spheres
Why hard spheres?
Hard Sphere Model
Typical Simple Material
The freezing transition of simple liquids can be well
described using a hard-sphere model
The Hard-Sphere Crystal
Face-Centered Cubic (FCC)
Simulation Details for Hard-Spheres
•Hard-sphere MD algorithmically exact: Chain of exactly resolved elastic
collisions
•Rappaport’s cell method: dramatically speeds up collision detection
•(100), (110) and (111) interfaces studied
•N ~ 10,000 particles
•Phase coexistence independent of T: rs31.037(crystal); 0.939 (fluid)
• gg1 kBT/s2
Results for hard-spheres
Davidchack & Laird, PRL 85 4751 (2000)]
g1000.592(7)kT/s2
g1100.571(6)kT/s2
g1110.557(7)kT/s2
How do these numbers fit in with other estimates?
•From Nucleation Experiments on silica microspheres
0.54 to 0.55 kT/s2 (Marr&Gast 94, Palberg 99)
•From Density-Functional Theory:
predictions range from 0.28 - 2.00 kT/s2
(WDA of Curtin & Ashcroft gives 0.62 kT/s2)
• From Simulation: Frenkel nucleation simulations: 0.61 kT/s2
Can Turnbull’s rule be explained with a hardsphere model?
For hard-spheres, Turnbull’s reduced interfacial free energy
scales linearly with the melting temperature
g = 0.57 (0.55) kTm/s2
grs-2/3 = 0.55 (0.53)Tm since rs 1.037s-3
If a hard-sphere model holds one would predict
that grs-2/3 = C Tm with C ~ 0.5-0.6
Hard-Sphere Model for FCC forming metals
grs-2/3 ~ 0.5 kTm
Turnbull’s rule follows since
DfusS = DfusH/Tm
is nearly constant for these
metals
Continuous Potentials
•Lennard-Jones
•Inverse-Power Potentials
Differences with Hard-spheres
w3 is non-zero
More care must be taken in construction of cleaving wall
potential
need to use NVT simulation to maintain coexistence
temperature throughout simulation (e.g., use Nose´-Hoover
or Nose´-Poincare methods)
Observation #4
In precise simulation work it is important to
always be aware of the damage done to statistical
mechanical averages by the discretization
Free energy simulations involving phase
equilibrium require highly precise simulations and
discretization error in averages can be important
Need: a detailed statistical mechanics of numerical
algorithms
Example of the effect of discretization error
in Nose´-Poincare´ MD
•In Nose´NVT dynamics, constant T is maintained by adding
two new variables to the Hamiltonian
p˜ i2
2
HN  
 V (q) 
 gkTln s
2
2m s
2Q
• g = Number of degrees of freedom
•NVE dynamics generated by HN , after time transformation dt/dt=s, yields a
canonical (NVT) distribution in the reduced phase space

•In Nose´-Poincare´ (Bond, Leimkuhler & Laird, 1999) the Nose´
Hamiltonian is time transformed to run in real time
H NP  s [ H N  H N (t  0)]
•Can be integrated using the Generalized Leapfrog Algorithm (GLA)
• GLA is symplectic
Example of the effect of discretization error
in Nose´-Poincare´ MD
•If canonical distribution is correctly obtained
then the equipartition theorem holds
T  Tinst
1
pi2


g k B mi
•The difference between T and Tinst is a measure
of the uncertainty in T due to the discretization
Nose´-Poincare´ with GLA this can be
•For
worked out (S. Bond thesis). Similar formulae
for Nose´-Hoover integrators
h2
A NP  A C   A2 C  A C 2 C
kT

1 
1 T –1
1
kT T –1
2 ( p,q)  pT M –1V 
( p M p  gkBT) 2  V T M –1V 
( p M p  gkBT)
12 
2Q
2
Q

Results for the truncated Lennard-Jones system
(Davidchack & Laird, J. Chem. Phys., 118, 7651 (2003)
kT/
Orientation
This work
g ( /s2)
Broughton &
Gilmer*
Morris &
Song
0.617 (tp)
(100)
0.371(3)
0.34(2)
0.369(8)
(110)
0.360(3)
0.36(2)
0.361(8)
(111)
0.347(3)
0.35(2)
0.355(8)
(100)
0.562(6)
N/A
-
(110)
0.543(6)
N/A
-
(111)
0.508(8)
N/A
-
(100)
0.84(1)
N/A
-
(110)
0.80(1)
N/A
-
(111)
0.75(1)
N/A
-
1.0
1.5
*J.Chem.Phys.
84, 5759 (1986).
Note that anisotropy in LJ differs from HS in that the order of
ANISOTROPY in Interfacial Free Energy
Expansion in cubic harmonics (Fehlner & Vosko)
3
  3




3
17
g ( nˆ )  g 0 1 1 ni4   2 ni4  66n12n22 n32  
5  i1
7 
 i1
Hard Spheres
Lennard-Jones System

T* = 0.617
1.0
1.5
g0(/s2)
0.360(2)
0.539(4)
0.798(6)
0.573(6)kT/s2
1
0.093(17)
0.13(3)
0.16(4)
0.09(4)
2
-0.011(4)
-0.022(9)
-0.019(6)
-0.005(11)
(g100 -g110)/g0
0.03(1)
0.035(15)
0.050(7)
0.036(16)
(g100 -g111)/g0
0.07(1)
0.10(2)
0.113(7)
0.061(16)
Inverse Power Potentials
s 
V(r)   
r 
n
•Important reference system for studying the effect of potential range
•For n=∞, we have the hard-sphere system.

•Phase Diagram: For n>7, crystal structure is FCC, for 4<n<7 system freezes
to BCC transforming to FCC at lower temperatures (higher densities)
n>7
fluid
fluid
n<7
fcc
bcc
fcc
density
Inverse Power Potential Scaling
s 
V(r)   
r 
n
•Only one parameter sn
•Excess thermodynamic properties only depend upon
n rs3 (kT/) -3/n = r* T* -3/n
•Phase diagram is one dimensional,
 only depends upon n
•Along coexistence curve:
Pc = P1T(1+3/n); gg1T(1 + 2/n)
Inverse Power Potentials and Turnbull’s rule
From the scaling laws
g  g1 T
(12 / n )
; n (solid) rs Tm
3 / n
 rs   T
3/n
n m
So…
g
˜  g rs2/ 3  g1 Tm12/ n [n (s)Tm3/ n ]2/ 3  g
˜  [g1 n2/ 3 (s)]Tm
as in hard spheres, we see scaling with Tm
And since DHfus TDSfus
g 2 / 3 
1 n,s
˜  
g
 DS 
DH fus  CT ,n DH fus
 fus 
Turnbull’s rule is
exact for inverse
power potentials!
Results for the Inverse-Power Series
Similar to Fe (Asta, et al.) the bcc interface has a lower
free energy
Turnbull Coefficient
Similar to Fe (Asta, et al.) the bcc interface has a lower
free energy
Results for the Inverse-Power Series
(Anisotropy)
Similar to Fe (Asta, et al.) the bcc interface has a anisotropy
Summary
• We have measured the crystal/melt interfacial free energy, g,
for hard-spheres, Lennard-Jones and inverse-power series.
Our simulations can resolve the anisotropy in this quantity.
• Comparison of data from fluctuation method and cleaving
method shows the two methods to be consistent and
complementary
• We show the molecular origin of Turnbull’s rule and give an
alternate formulation
• For the inverse-power series gbcc < gfcc consistent with
fluctuation model calculations and nucleation experiments also bcc is less anisotropic than fcc.