ATOMIC-SCALE DYNAMIC PROCESSES
IN THE BRITTLE FRACTURE OF SILICA
THOMAS P. SWILER
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
Thomas P. Swiler
This dissertation is dedicated to my parents, Thomas Russell Swiler and Elaine
Palda Swiler, who were always forgiving of my early "fracture experiments." Until we
can fulfill our destructive tendencies through simulation, we must rely on forgiveness.
I would like to acknowledge the support of my advisor, Joseph H. Simmons,
throughout this work. He has continually supported a molecular dynamics effort even
though it has not been supported by external sources. The opportunity to explore the
implications of this simulation work without having to provide specific answers to
questions posed by sponsors has made work in this area interesting and instructive.
I would like to thank John J. Mecholsky, Jr. for enlightening me to current
research in brittle fracture. His observations and interpretations of fractal fracture
surfaces have provided me with the background to make the link between computer
simulated fracture and experimentally observed fracture.
I would like to thank Adrian C. Wright for performing comparisons of our pair
correlation function data derived from molecular dynamic simulations to his radial
distribution data derived from neutron diffraction experiments.
Finally, I would like to thank my family and friends who have provided support
for me during the up and downs of writing this dissertation.
TABLE OF CONTENTS
ACKNOWLEDGEMENTS .............. ................... .. iv
LIST OF TABLES ..................................... ix
LIST OF FIGURES .................. .................... x
KEY TO SYMBOLS ............ .............. ........... xv
Abstract .............. .................. ............ xxiii
1 INTRODUCTION ............... .......... ........ 1
W hat is Fracture .......................... ......... 1
Fracture as a Transformation .............. ... ............ 2
Molecular Dynamics ................... ............ 4
Statement of Problem ............. .................... 7
2 BACKGROUND ............. ............ .......... 8
Brittle Fracture ............... ................... 9
Characteristics of Brittle Fracture ........... ..... ..... 9
The Ultimate Strength of Materials ........... ..... .... 9
The Effect of Flaws ........... ...... ......... . 11
Regimes of Brittle Fracture ............. .... ....... 12
Environmental Effects ................ . .......... 12
Crack Propagation Velocities ............. .... ....... 14
Reading the Fracture Surface .......... ...... ..... 14
Silica ............................... ........ 16
Bonding in Silica ...................... ........ 16
Structure of Silica ................... ...... .... 16
Phases of Silica .............................. 18
The Crystal Structure of Cristobalite ............... ... 19
The Ring Structure of Silica ............ ..... ....... 20
Elastic Properties of Vitreous Silica . . . . . . ..... 23
The Fracture Strength of Silica . . . . . . . . .. 24
Stress Corrosion of Silica . . . . . ..... ........ 25
Molecular Dynamics ................... ....... ... . . 28
Survey of Applications .................... .... .. . 28
Survey of Potential Functions . . . . . . . ..... . ... 31
Potentials for Silica . . . . . ........ .. .. . . 34
Applications of MD in Models of Silica . . . . . ..... 37
Comparisons of Simulated Structure to Experiment . . . ... 39
MD Simulations of Fracture . . . . . ... . 40
Chaos in Fracture: Determinism and Fractal Surfaces .......... . 43
Investigations of Chaotic Processes . . . . . . . ..... 44
Predictability of Brittle Fracture . . . . . . . . 46
Fractal Geometry ................... ........... 47
Applications of Fractal Geometry to Fracture . . . . .... 49
The Fractal-like Nature of Fracture Surfaces . . . . .... 51
Applicability of Fractal Analyses to MD Simulated Systems . . 52
3 PROCEDURE .................. ......... ......... 56
Using Molecular Dynamics to Model Fracture . . . . . . . . 56
Potential Functions ........... .... ... .......... 58
The Bom-Mayer-Huggins Potential . . . . . . . ..... 60
Truncation or Screening Functions . . . . . . . ..... 61
The Soules Potential Function . . . . . . . . .. 63
The Feuston-Garofalini Potential Function . . . . . ..... 65
Simulation Techniques ........ . . . . . . . . ...... 68
Integration Algorithms ... . . . . . . . ....... 68
Verlet algorithm ................... ...... 68
Fifth-order Gear predictor-corrector algorithm . . ... 69
The Time Step ........... ...... .. ......... 69
Limiting Particle Interactions . . . . . . . 70
Application of Periodic Boundary Conditions . . . . . ... 72
Handling Particle Interactions . . . . . . . . 74
Observing the System . . . . . ......... . . .. 76
System thermodynamics . . . . . . . . 76
System structure .................. ........ 77
Force Distributions . . . . . . . .. 79
Damped Motion Molecular Dynamics . . . . . . . ... 80
Assigning Initial Velocities . . . . . . . . .. 80
Analytical Techniques . . . . . . ....... . . . ..... 81
Graphic Representations . . . . . . . . .. 82
Free Volume Sphere Analysis . . . . . . . .. 82
Ring Statistics ...................... .......... 84
Structural Defects ................. ......... 86
Composition Profiles ........................... 86
Mechanical Properties .......................... 87
Program Details ................ . . . . . . . .... 87
Overview .................................. 87
The Structure of Computer Code . . . . . . ...... 88
The Choice of a Programming Language . . . . . ..... 92
The Simulation Program . . . . .......... ...... 95
The Graphic Analysis Program . . . . . . . .... . 96
Other Analysis Programs . . . . . ....... . . . 97
Experimental Techniques ................... . . . . 97
Initial Conditions ................... ........... 97
Application of Uniaxial Strain . . . . . . . ... . 99
Application of Shear Strain ........... . . . . 101
Experiments and Analyses . . . . . ........ . . . . ..... 103
Phases Simulated ................... . . . . 103
Strain Rates Simulated . . . . . . . . . ...... 103
Comparison of Potential Functions .. . . . . . . . 103
Thermal Effects ............................. 104
Studies of Determinism . . . . .......... ........ 104
Studies of Nonlinear Viscosity . . . . . . . 105
4 RESULTS AND DISCUSSION ................... ...... 106
The Models ............................. ........ 106
Static (Structural) Characteristics . . . . . . . 107
Densities and pressure tensors . . . . . . . ... 107
Low-order pair-correlation functions . . . . . .... 108
High-order pair correlation functions . . . . . ... 112
Discussion of higher order peaks . . . . . . ... 113
Bond angle distributions . . . . . . . . . . 116
Comparison of structural distributions to experiment . .. 117
Ring statistics ........ .................. 121
Volume distributions . . . . . . . . ... 124
Defects ......... ..................... 125
Dynamic Properties ............... ............. 126
Melting temperature . . . . . . ... ....... 126
Mechanical properties . . . . . . . . ... 129
Review of Models .................. ......... 131
Fracture by Application of Uniaxial Strain . . . . . . . . ... 132
Overview ................. ........ .......... 132
Low Cristobalite ................... ........... 133
Stress-strain curves . . . . ......... . . 134
Energetics of fracture . . . . . . . . . 142
Fracture surfaces . . . . . . . . 149
Vitreous Silica. ...... .. ........... 154
Stress-strain curves . . . . . . . 154
Energetics of fracture ....................... 160
Structural changes during strain-to-fracture . . . 164
Viscous Flow ........... .............. ........... 174
Determinism in Fracture ....................... ... .. .. 179
Reproducibility and Accuracy of Computer Simulations . . . 179
Variation in Anneal Times ..... . . . . . . . . 180
Round-off Errors ....... .. . . . . . . . . . 182
"Butterfly Effect" ........... ... ............ .. 183
Sensitivity on Simulation Artifacts . . . . . . ..... 185
5 ANALYSIS ......... ................... ........ .. 187
M odels of Silica .................. . .. ........... 187
Stress and Strain ............... ...... ............. 190
Structural Changes During Strain-to-Failure . . . ... . . . 190
Shear Flow .......... .. . .... . . .......... 193
Initiation of Fracture . . . . . .. ... . . . . ...... 194
Crack Propagation .......... . . . . . . . . . 198
Simulation of Stochastic Crack Growth . . . . . . . . .. 199
The Completion of Fracture . . . . . ........ . . . . 202
The Structure of the Fracture Surface . . . . ... . . . ... 203
When is a Bond Broken ..... ........ ......... .... 203
Surface Reconstruction . . . . ........ . . . . 204
Determinism in Fracture ...... . . . . . . . . . . 207
6 CONCLUSIONS AND FUTURE WORK .......... . . . .. . 209
REFERENCE LIST ................ ...... .......... .. 211
BIOGRAPHICAL SKETCH ........ . . . . . . . . . 219
LIST OF TABLES
2-1 Characteristics of the Phases of Silica . . . . . . . . . 19
3-1 Parameters Used in the Soules Force and Potential Functions ...... . 63
3-2 Parameters used in the Feuston-Garofalini Potential Function . . ... 67
4-1 Coordination of species in simulations using various potentials and
methods of formation as found by the ring statistics program. ...... 125
4-2 Coordination of species in simulations of vitreous silica at 0.2 strain and
various strain rates as found by the ring statistics program. . . ... 173
4-3 Viscosities of vitreous silica as modeled by the Soules potential at various
temperatures and strain rates, in units of 1021 Pas. . . . . .... 175
LIST OF FIGURES
2-1. Schematic showing crack propagation velocities as a function of stress
intensity factor.. ................. . . . . . . . . ...... 13
2-2. The silica tetrahedron is the basic structural unit of silica structures, a)
One silica tetrahedron. b) Two tetrahedra linked by comer sharing. 0 is
the Si-O-Si bond angle.. ................... . . . ... 17
2-3. Chemical reaction proposed by Michalske and Freiman for
environmentally enhanced Si-O-Si bond rupture. The strained crack tip
bond (a) adsorbs a water molecule, (b) is cleaved by dissociative reaction,
and (c) is converted to surface silanol groups.. . . . . . . ... 26
3-1. Plot of the Soules potential function for various interactions. . . ... 64
3-2. Plot of the two-body portion of the Feuston-Garofalini potential function
for various interactions. .................. ... . . . . 66
3-3. The effect of periodic boundary conditions in a simulated system. Note
that the neighborhood of an atom may extend beyond the primary cell. 71
3-4. Schematic showing the application of periodic boundary conditions to a
system under shear strain. ........ . . . . . . . . 74
4-1. Comparison the first peaks of the Si-O pair correlation functions for
vitreous silica as modeled by the Soules and the F-G potential
functions . . . . .. . . .. .. .. . . . . . . . . . 109
4-2. Comparison the first peaks of the 0-0 pair correlation functions for
vitreous silica as modeled by the Soules and the F-G potential
functions ................... ............ ......... 110
4-3. Comparison the first peaks of the Si-Si pair correlation functions for
vitreous silica as modeled by the Soules and the F-G potential
functions ................... ......... ........... 111
4-4. Comparison the higher order peaks of the Si-O pair correlation functions
for vitreous silica as modeled by the Soules and the F-G potential
functions ................... ................... 113
4-5. Contribution of various pair interactions to the potential energy of vitreous
silica as modeled by the Soules potential function. . . . . . ... 114
4-6. Contribution of various pair interactions to the potential energy of vitreous
silica as modeled by the F-G potential function. . . . . . .... 115
4-7. Comparison of O-Si-O bond angle distributions for vitreous silica as
modeled by the Soules and the F-G potential functions. . . . . ... 117
4-8. Comparison of Si-O-Si bond angle distributions for vitreous silica as
modeled by the Soules and the F-G potential functions. . . . . ... 118
4-9. Comparison of radial distribution functions derived from a neutron
diffraction experiment (solid line) and a molecular dynamics simulation
using the Soules potential (dashed line). . . . . . . ...... ...... 119
4-10. Comparison of radial distribution functions derived from a neutron
diffraction experiment (solid line) and a molecular dynamics simulation
using the F-G potential (dashed line). . . . . . . . . . 120
4-11. Comparison of the ring size distributions for vitreous silica as modeled by
the Soules and the F-G potential functions and as formed by crystal melt
and random position techniques. . . . . . ..... ....... 122
4-12. Comparison of the free volume sphere size distributions for vitreous silica
as modeled by the Soules and the F-G potential functions and as formed
by various techniques. ........ . . . . . . . . . 124
4-13. Potential energy as a function of temperature for cristobalite as modeled
by the Soules potential. ......... . . . . . . . ....... 127
4-14. Potential energy as a function of temperature for cristobalite as modeled
by the F-G potential .................. .............. 128
4-15. Equilibrium elastic stress-strain curves for cristobalite and vitreous silica
as modeled by the Soules and the F-G potential functions. . . . ... 129
4-16. Comparison of stress-strain curves at various strain rates for low
cristobalite as modeled by the Soules potential. . . . . . . ... 134
4-17. Comparison of stress-strain curves at various strain rates for low
cristobalite as modeled by the F-G potential. . . . . . . .... 135
4-18. View down the C-axis of low cristobalite under an intermediate tensile
strain applied along the C-axis, as modeled by the Soules potential. In
this view and in all following views, oxygen atoms are represented by the
larger, lighter balls, and silicon atoms are represented by the smaller,
darker bals. ..................... ........ ........... 137
4-19. View down the C-axis of low cristobalite under an intermediate tensile
strain applied along the C-axis, as modeled by the F-G potential. . . 138
4-20. View down the C-axis of low cristobalite under a high tensile strain
applied along the C-axis, as modeled by the F-G potential. . . . ... 139
4-21. View down the C-axis of low cristobalite under a high tensile strain
applied along the C-axis, as modeled by the Soules potential. ...... .140
4-22. Potential energy as a function of strain for low cristobalite as modeled by
the Soules potential function under adiabatic strain at various strain
rates. ................... .............. ........ 144
4-23. Kinetic energy as a function of strain for low cristobalite as modeled by
the Soules potential function under adiabatic strain at various strain
rates. ................... ............. ......... 145
4-24. System temperature as a function of strain for low cristobalite as modeled
by the Soules potential function under adiabatic strain at various strain
rates. ................... ............. ......... 146
4-25. Schematic showing the flow of energy in a brittle material system during
strain-to-failure ........... . . . . . . . . . ..... 147
4-26. Potential energy as a function of strain for low cristobalite as modeled by
the Soules potential function under isothermal strain at various strain
rates. ................... ............ .......... 148
4-27. Normal view of a shallow slice of a surface created during the fracture of
low cristobalite as modeled by the Soules potential under isothermal
conditions at 0.05 ps' strain rate. The failure of the slice to fill the entire
area is a consequence of the roughness of the surface. . . . . ... 149
4-28. Side view of surfaces created during the fracture of low cristobalite as
modeled by the Soules potential under isothermal conditions at 0.05 ps'
strain rate .......... ....... ...... .............. 151
4-29. Side view of surfaces created during the fracture of low cristobalite as
modeled by the Soules potential under isothermal conditions at 0.5 ps'-
strain rate .................. ...... .............. 151
4-30. Relative abundance of oxygen and silicon atoms in slices along the length
of the cell perpendicular to the direction of strain. Values plotted are the
predominance of oxygen over silicon. . . . . . . . . .. 152
4-31. Close-up side view of surface created during the fracture of low
cristobalite as modeled by the Soules potential under isothermal conditions
at 0.5 ps' strain rate ................... .. ........... 153
4-32. Comparison of stress-strain curves at various strain rates for vitreous silica
as modeled by the Soules potential. . . . . . . . . 155
4-33. Comparison of stress-strain curves at various strain rates for vitreous silica
as modeled by the F-G potential. ..... . . . . . . 156
4-34. Distribution of forces on the atoms comprising a vitreous silica sample
under strain at 0.05 ps' strain rate. Various curves represent different
engineering strains. ................. . . . . . 158
4-35. Distribution of forces on the atoms comprising a vitreous silica sample
under strain at 5.0 ps' strain rate. Various curves represent different
engineering strains. ......... . . . . . . . . . . 159
4-36. Potential energy as a function of strain for vitreous silica as modeled by
the Soules potential function under adiabatic strain at various strain
rates. ................................... ..... 161
4-37. System temperature as a function of strain for vitreous silica as modeled
by the Soules potential function under adiabatic strain at various strain
rates. ................... ..................... .. 162
4-38. Potential energy as a function of strain for vitreous silica as modeled by
the Soules potential function under isothermal strain at various strain
rates. ................... .................... .. 163
4-39. Si-O pair correlations in vitreous silica. Comparisons between unstrained
and strained to 0.2 at various strain rates are made. . . . . .... 165
4-40. Si-Si pair correlations in vitreous silica. Comparisons between unstrained
and strained to 0.2 at various strain rates are made. . . . . .... 166
4-41. O-Si-O bond angle distributions in vitreous silica. Comparisons between
unstrained and strained to 0.2 at various strain rates are made. ...... 167
4-42. Si-O-Si bond angle distributions in vitreous silica. Comparisons between
unstrained and strained to 0.2 at various strain rates are made. . . . 168
4-43. Comparison of ring size distributions of vitreous silica strained to 0.2 at
different strain rates. ........... . . . . . . . . 169
4-44. Change in the abundances of rings of various sizes during strain to 0.2 in
vitreous silica at various strain rates. . . . . . . . 170
4-45. Free volume sphere histogram--the integral of the free volume sphere
distribution--for vitreous silica strained to 0.2 strain at various strain
rates. ................... ........... . ........ . 171
4-46. Shear stress as a function of time for a vitreous silica sample as modeled
by the Soules potential function at various temperatures. . . . ... 176
4-47. Shear stress as a function of time for a vitreous silica sample as modeled
by the Soules potential function at various strain rates. . . . . ... 177
4-48. Stress-strain curves for nearly identical systems after various anneal
times. ........... .... ......... ............ 180
4-49. Stress-strain curves for nearly identical systems differing by a loss of
precision resulting from truncation of system conditions at various
times. ................ . ........ ........ .. 183
4-50. Deviations from stress-strain curves resulting from a displacement of one
atom immediately before straining. . . . . . . . . 184
4-51. Deviations from stress-strain curves resulting from a displacement of one
atom before an anneal period prior to straining. . . . . . . ... 185
5-1. Pseudo fracture surface created by simple stochastic crack propagation
simulation ......... .... .................. 201
5-2. Schematic showing the fracture of a bond in vitreous silica. Series shows
unstrained bond, strained bond, and broken bond. . . . . . ... 205
KEY TO SYMBOLS
A 1. A parameter in the Born potential. 2. The cross sectional area covered
by a bond. 3. A parameter in one form of the Lennard-Jones potential.
4. A parameter in the two-body portion of the Stillinger-Weber potential.
A,,, A parameter in the two-body part of the Vessal potential for types u and
A,, 1. A parameter in the Born-Mayer-Huggins interaction potential between
atoms of type u and v. 2. A parameter in the two-body part of the Vessal
potential for types u and v.
A,, A parameter in the three-body part of the Vessal potential corresponding
to the interaction between types u, v and w.
B 1. A parameter in the Born potential. 2. A parameter in one form of the
Lennard-Jones potential. 3. A parameter in the two-body portion of the
Stillinger-Weber potential. 4. A parameter in a stocastic crack simulation.
B(9Oy An angularly dependent function used in the Vessal potential.
BI., A parameter in the two-body part of the Vessal potential for types u and
B,, The strength of the three-body interaction of types u, v and w in the
C 1. A constant the maximum crack velocity equation. 2. A parameter in
a stocastic crack simulation.
C11, A parameter in the two-body part of the Vessal potential for types u and
C,, 1. A parameter in the Born-Mayer-Huggins interaction potential between
atoms of type u and v. 2. A parameter in the two-body part of the Vessal
D 1. Fractal dimension. This is an approximation of the box counting
dimension when applied to real materials. 2. A parameter in a stocastic
D" Fractional part of a fractal dimension.
D,,, A parameter in the two-body part of the Vessal potential for types u and
D,, A parameter in the Bom-Mayer-Huggins interaction potential between
atoms of type u and v.
E,, A parameter in the two-body part of the Vessal potential for types u and
E Young's modulus or elastic modulus.
EK System kinetic energy.
Ep System potential energy.
F A mathematical set of points.
F(r) Force function of r, based on a potential function of r.
F,,, A parameter in the two-body part of the Vessal potential for types u and
F, The interatomic force vector between atoms i and j.
Fqk The interatomic force vector resulting from the interaction between atoms
i, j and k.
G Shear modulus.
H,, The strength of the steric repulsion in the interaction of types u and v in
the Vashishta potential.
Kic Fracture toughness, the stress intensity required to cause type I fracture.
L 1. A linear dimension that is smaller than an edge dimension of the
molecular dynamics cell, used in the truncation function of the Born-
Mayer-Huggins potential based on the Ewald sum. The value of the force
cut-off distance is often used. 2. Longitudinal modulus.
M, A material dependent constant in the equation relating the failure stress of
a material and the radius at which a transition in the fracture surface
occurs, such as mirror to mist, mist to hackle, or mist to crack branching.
The subscriptj corresponds to the particular transition.
N The number of atoms in a simulation.
N, The number of atoms of type i.
N,(F) A function that gives the number of boxes of side length 6 that cover the
P System pressure.
P_ Pressure tensor in the mn coordinates.
Pi-, A parameter in the two-body part of the Vessal potential for types u and
Q,,, A parameter in the two-body part of the Vessal potential for types u and
R, A random number evenly distributed between 0 and 1 used in assigning
random velocities to atoms.
R,,, A parameter in the two-body part of the Vessal potential for types u and
R? The crystallographic R factor, a figure of merit used in the comparison of
pair correlation functions where better fits are represented by lower
Sl, A parameter in the two-body part of the Vessal potential for types u and
T Temperature of the system.
T(r) The truncation function applied to the force function in the Soules
Tp(r,) The experimentally derived pair-correlation function at discrete values of
T,,,(r,) The simulation derived pair-correlation function at discrete values of r.
V Volume of the molecular dynamics cell.
Z, The effective charge of an ion of type u in the Vashishta and Born-Mayer-
a 1. A parameter in the Morse potential. 2. A cut-off distance used in the
ao 1. The dimension of a unit cell along the crystallographic a direction. 2.
A parameter having dimensions of length in the relation between fracture
toughness and surface fractal dimension of for a material.
aY The effect of a defect at location (i,j) in a stocastic crack simulation.
b A parameter in the Born-Mayer-Huggins potential.
c Crack length.
co The dimension of a unit cell along the crystallographic c direction.
c. An integer indicating how many cell lengths a particle should be translated
in the u dimension to bring it back into the primary cell.
d, The deviation caused by a defect at location (i,j) in a stocastic crack
dj The u coordinate difference between atoms i andj.
d' The u coordinate difference between atoms i and j after correction for
periodic boundary conditions.
dimF The box-counting dimension of set F.
f2(r) A dimensionless two-body potential, used in the Stillinger-Weber
fo(r1,r2,0) A dimensionless angularly dependent three-body potential, used in the
gdev(1) A function that returns a random Gaussian deviate with a variance of 1.
i 1. An index indicating an atom. 2. An index representing a location on
a discrete x-y plane in a stocastic crack simulation.
j 1. An index indicating an atom. 2. An index representing a location on
a discrete x-y plane in a stocastic crack simulation.
k An index indicating an atom.
k, Boltzmann's constant.
k,, A parameter in the three-body part of the Vessal potential corresponding
to the three-body spring constant for the interaction between types u, v
I An unspecified constant in the Vashishta potential.
m, The mass of atom i.
m. The mass of an atom of type u.
m An index representing a coordinate in space.
n 1. A parameter in the Born potential. 2. An exponent in the truncation
function used in the Soules force function. 3. An index representing a
coordinate in space.
n, The number of valence electrons for atom type u, used in the Born-
p A parameter in the two-body portion of the Stillinger-Weber potential.
q A parameter in the two-body portion of the Stillinger-Weber potential.
q. The charge on atoms of type u.
r Length or distance, often corresponding to an interatomic distance.
ro The cutoff distance for the three-body interaction in the Vashishta
r4, The decay length of the r" interaction term in the Vashishta potential.
r, The position vector for atom i.
ro, The position vector for atom i at time t.
rv The interatomic distance between atoms i andj.
Radius about the initiating flaw where a transition in the fracture surface
occurs, such as mirror to mist, mist to hackle, or mist to crack branching.
The subscriptj corresponds to the particular transition.
The interaction cutoff distance used in the Soules force function.
The u dimension of the MD cell.
1. An index indicating an atom type. 2. An index representing a
coordinate in space.
1. An index indicating an atom type. 2. An index representing a
coordinate in space.
The velocity of atom i.
The maximum crack propagation velocity in a material.
1. An index indicating an atom type. 2. An index representing a
coordinate in space.
The x coordinate location on a discrete x-y plane in a stocastic crack
The x coordinate distance between atoms i and j.
The u Cartesian coordinate of atom i.
The y coordinate location on a discrete x-y plane in a stocastic crack
The y coordinate distance between atoms i and j.
The z coordinate of a stocastic crack at coordinates (x,y).
The z coordinate distance between atoms i and j.
The length of a time step.
The electronic polarizability of an ion of type u in the Vashishta potential.
1. The Pauling coefficient, used in the Born-Mayer-Huggins potential, for
atom types u and v. 2. A parameter used in the screening function for the
long range interactions in the Feuston-Garofalini potential between atom
types u and v.
6 A small length.
e 1. Strain. In this work we use engineering strain. 2. A parameter used in
Lennard-Jones and Stillinger-Weber potentials to scale the potential
energy. 3. A constant in the truncation term of the Born-Mayer-Huggins
potential function based on the Ewald sum, typically between 0.175 and
Y 1. Surface energy. 2. A parameter in the three-body portion of the
X A parameter specifying the strength of the three-body interaction in the
X,, A parameter that gives the strength of the three-body interaction for atom
types u, v and w in the Feuston-Garofalini potential.
,1- The exponent of the steric repulsion in the interaction of types u and v in
the Vashishta potential.
o0 A parameter in the Morse potential.
42,,(r) Pair interaction potential function of r for atom types u and v. The
subscript 2 indicates that this is a two-body interaction potential. If only
one atom type is present, the subscripts uv are not written.
43.,-(r,, r,, r Three-body interaction potential function of positions r,, r, and r, for atom
types u, v and w. The subscript 3 indicates that this is a three-body
interaction potential. If only one atom type is present, the subscripts uvw
are not written.
O1.w(ri, rIk,' 0,) Angularly dependent three-body potential of r,, r,, and the
included angle, 0k for atom types u, v and w. The subscript 0
indicates that this is a simple angularly dependent three-body
potential. If only one atom type is present, the subscripts uvw are
p 1. Density. 2. Hardness parameter in the Born-Mayer-Huggins potential.
a 1. Stress. 2. A parameter used in Lennard-Jones and Stillinger-Weber
potentials to scale the interaction distances.
ac Critical stress required to cause fracture in a sample with surface flaws.
of Failure stress.
aTh Theoretical cohesive strength.
a, Radius repulsive parameter for atom type u in the Bom-Mayer-Huggins
00 The tetrahedral bond angle, cos-'(-1/3), in the three-body portion of the
0,k The i-j-k bond angle.
0,, The preferred bond angle for the three-body interaction of types u, v and
w in the Vashishta potential.
r7, Shear in the u-v plane: displacement in u dependent on change in v.
v Poisson's ratio.
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
ATOMIC-SCALE DYNAMIC PROCESSES
IN THE BRITTLE FRACTURE OF SILICA
THOMAS P. SWILER
Chairperson: Joseph H. Simmons
Major Department: Materials Science and Engineering
The molecular dynamics (MD) simulation technique was used to study the atomic-
scale dynamics processes that may take place in the fracture of brittle materials. Silica
was chosen as the host material for this study because it occurs in both crystalline and
vitreous phases, it is of great technical interest, and it has been widely simulated using
Two models for silica based on different interatomic potential functions were used
and compared. The Soules two-body potential resulted in structures having better
densities and closer correspondence to neutron diffraction data than did the Feuston-
Garofalini three-body potential, which resulted in fewer topological defects and a
cristobalite structure that did not undergo a symmetry transformation under uniaxial
strain. Both vitreous silica and cristobalite were modeled.
Fracture experiments were performed by applying uniaxial strain to the various
samples at strain rates ranging from 0.05 ps' to 5.0 ps'. The fracture behaviors for all
samples were dependent on the strain rate. Samples strained at higher strain rates
exhibited higher fracture strengths, less localized fracture, and rougher fracture surfaces.
Cristobalite samples were stronger than vitreous silica samples.
Both permanent structural variations and thermal vibrations affected the path of
fracture through the samples. The fracture path was extremely sensitive to initial
conditions, indicating that fracture is a truly chaotic process.
We found that the effects of the strain-rate dependence could be traced to the
"flow" of strain throughout the sample during fracture. In a normal fracture process, the
local strain-rate at the crack tip may increase as a crack propagates, causing fracture to
appear as a high strain-rate process in real material. A fracture mechanism is proposed
that may result in a fracture surface containing self-similar features like those observed
in the fracture of real brittle materials. A simple stochastic crack growth simulation is
performed to show the viability of the proposed mechanism.
What is Fracture
The word fracture originated from the latin wordfrangere, which simply means,
to break. In days of antiquity, brittle objects that were whole simply broke into parts if
they were not strong enough to support the applied load. Of course, objects still break
if they are not strong enough to support the applied load, but we now have the tools to
investigate how they break. Fracture is not considered an instantaneous event by which
one piece becomes two or more pieces, but a transformation process that can be followed
from start to finish. However, as much as fracture has been experienced, the process
is not fully understood.
We propose that a complete understanding of fracture can be obtained only if we
understand what happens in a material under strain to fracture at the atomic scale.
Fracture displays a complexity that cannot be fully modeled using a continuum model for
a material body. Griffith' and Orowan2 proposed criteria for predicting the growth of
cracks in brittle materials using strain-energy release and stress-intensity criteria,
respectively, in a body that is assumed to be composed of a continuum of material. The
formalism put forth by these models is still dominant in fracture mechanics. The stress
intensity criterion has been extended to predict slow crack growth of materials in a
corrosive environment by Michalske and Freiman.3 However, although these models
predict the growth of cracks in an object under an applied load and the role of cracks in
reducing the strength of objects, they cannot predict the complex geometries that the
fracture process creates. Thus, we must advance beyond a continuum model and view
fracture as a process that occurs in a many-body system, instead of as a process that
occurs in a continuum of matter, to understand all the processes that result in the details
observed in fracture.
We seek to model fracture on the atomic scale and then attempt to extrapolate the
behaviors that we observe to the mesoscopic scale.
Fracture as a Transformation
Fracture is a simple process in some respects, and it is nearly infinitely complex
In terms of fracture, materials fall into two gross categories--brittle and ductile.
Brittle materials are those materials that fracture with little deformation before the
initiation of fracture. Ductile materials are those that are capable of undergoing
significant deformation before fracture occurs. Certain materials may be either brittle
or ductile under different conditions.
Within the field of materials science and engineering, the costs associated with
both fracture and corrosion to the national economy on a yearly basis are often cited to
emphasize the negative aspects of these processes. These are the processes that cause
the failure of useful objects. As a result, much of the research in the area of fracture has
been in involved in predicting the strength of materials.
The processes by which the usable strengths of materials in real-world
applications are limited are well understood. The strengths of ductile materials are
limited by the formation and movement of dislocations that cause permanent deformation
at the yield stress. This may cause the parts made from these materials to fail in their
purposes, whether the parts have completely fractured or not. The strengths of brittle
materials are usually limited by flaws that result in regions of significantly higher than
average stresses within the materials, thus causing fracture to be initiated even when the
average stress on the material is less than its cohesive strength.
How do ductile materials fracture when dislocations are pinned, and how do
brittle materials fracture when there are no flaws? Dislocations in engineering ductile
materials are routinely pinned by work hardening and by precipitation hardening so that
processes giving ductility are deactivated and failure occurs by brittle fracture. When
there are no flaws present in brittle materials, the fracture of these materials becomes
dependent on the dynamic aspects of materials.
The importance of the dynamics of failure in ductile materials has long been
accepted. As a result, the dynamics of failure in ductile materials has been extensively
studied. However, the dynamics of failure in brittle materials is not so well understood,
because brittleness has traditionally implied that fracture is complete almost as soon as
it has begun.
In order to study the dynamics of fracture in brittle materials, we require high-
speed investigative techniques on the atomic scale. These techniques do not currently
exist. Thus, in order to study the dynamic processes occurring in brittle fracture, we
perform simulations to examine the behaviors that we cannot see by normal experimental
Molecular dynamics (MD) is a numeric technique that builds on a knowledge of
interatomic potential functions to investigate the behavior of a system of atoms under
various conditions. It was first used by Rahman4 to investigate the motion of Ar atoms
in a liquid, and has since been used in studies ranging from enzymatic reactions to the
fracture of materials.
Molecular dynamics can be compared with another commonly used simulation
technique, Monte Carlo simulations, by examining the path that the system takes while
being simulated. Monte Carlo simulations make random perturbations of a system to
find favorable configurations. It is not necessary that the individual perturbations be
realistic because only the thermodynamics and structure of the final state are important.
In contrast, molecular dynamics simulations allow movement of atoms according to
Newton's equations of motion. Thus, the Monte Carlo simulation technique allows the
simulation of physically realistic structures, but not physically realistic dynamics.
Molecular dynamics, on the other hand, is capable of simulating dynamics, but not
necessarily the most optimal structures. In MD, movement to optimal structures may be
blocked by either inertial barriers, where particles may possess too much kinetic energy
to stay in the deepest potential wells, or potential barriers, where particles may not have
enough kinetic energy to get out of shallow potential wells, within the short time allowed
by the simulation.
Although the limitations of molecular dynamics are numerous, its strengths are
obvious. Molecular dynamics is a technique that can be used to study the atomic
processes and the resulting properties of material systems that are not conducive to a
study by analytic means. A quick inspection of the disordered surfaces of fragments of
a fracture event illustrates the complexities involved in fracture. Thus, molecular
dynamics offers the ability to witness the events that occur during fracture that result in
these surface morphologies.
Processes in noncrystalline solids such as vitreous silica cannot be understood
without such a simulation technique because of the complexities involved. Traditionally,
the properties of glasses have been described by thermodynamics because of an inability
to understand the dynamics involved. Most glasses have an enthalpy function of
temperature that displays a glass transition region. The glass transition region is the
temperature region wherein the enthalpy in a super-cooled liquid ceases to behave as a
liquid and begins to behave as a solid upon cooling. The laws of thermodynamics tell
us that within the glass transition region there must be certain limits on the changes in
enthalpy. For example, since both the enthalpy and the specific heat of a liquid are
greater than those for a crystal, a linear extrapolation of the enthalpy function for a liquid
to lower temperatures may cross the enthalpy function for the crystalline phase.
However, because of the excess configurational potential energy of a liquid as compared
to a crystal, the enthalpy of a liquid should never be lower than the enthalpy of a crystal
at the same temperature. This is the Kauzman Paradox. In a sense, it guarantees that
the liquid of the composition must go through a glass transition if it does not crystallize
first. After a material passes through a glass transition upon cooling, its structure
becomes configurationally arrested, and both glass and crystal phases have nearly the
same specific heats. At this point the glass and crystal enthalpy curves become parallel,
so that the glass curve will not cross the crystal curve, therefore avoiding the paradoxical
situation. Simulation is important in the investigation of these transformations because
one can pinpoint the types of liquid-like behaviors possible at different temperatures,
rather than lumping behaviors into a featureless normal distribution. Thus,
nonequilibrium molecular dynamics has become an important technique in the study of
It is in this backdrop that simulations of fracture in a brittle material, a glass,
vitreous silica are performed. Molecular dynamics not only gives us the ability to
simulate a dynamics event such as fracture, it provides us also with a way to obtain the
structure, in the case of a glass, so that these simulations can be performed on facsimiles
of real materials. Limited as it may be, MD gives us the ability to study the events that
we cannot view directly.
Statement of Problem
Vitreous silica is a material in which the intrinsic fracture behavior has become
important. Silica optical fibers are routinely made to near theoretical strength by
eliminating surface and bulk flaws. Despite this control over the state of flaws in this
material, variations in the strength of samples still exist. The Griffith criterion, when
applied to the variation in sample strengths, predicts flaw sizes on the order of a few
angstroms. This work seeks to find the source of this variation in strength in this
Additionally, the dynamics of fracture, although not directly observable in fracture
experiments, leave a signature on the resulting fracture surfaces. Reading this signature
is the subject of the science of fractography. Fractography has been successful in
locating the flaw that initiated fracture and in finding the state of stress in the material
when fracture was initiated. However, the dynamic processes that result in the exact
path that fracture takes after being initiated have remained elusive. The simulation
technique that we use should be able to shed some light on these processes.
In this work, we use molecular dynamics simulations to study the atomic-scale
fracture processes in silica. Therefore, we will review the important aspects of this
study. Since we hope to generalize the results of our brittle fracture simulations to brittle
fracture in general, we will start with a discussion of fracture. Then we will examine
the various phases of silica that are at the center of our particular study. We will then
discuss the technique of molecular dynamics, including the limitations of this technique
as well as the extent of applications. Various simulation models for silica will then be
reviewed. We will then describe how fracture is simulated and studied using molecular
dynamics. Finally, we will identify the aspects of fracture that a dynamical technique
such as MD may allow us to study. Because of the various fields that are drawn
together, a smooth progression between topics will not be attempted in the background
Characteristics of Brittle Fracture
Brittle fracture is the mechanical failure of a material without plastic deformation.
Brittle fracture is not limited to traditionally brittle materials. Metals or polymer
materials that are considered to be plastic or ductile can be brittle under certain
conditions of temperature and strain rate. Even inorganic glasses that are the epitome
of brittle solids at room temperatures behave as viscous liquids at temperatures above the
glass transition and are not brittle at low strain rates.
Experimentally, brittle fracture is generally a sudden process, so the dynamics of
this process are not well known. Observations made during the fracture process by
Dickinson et al.5 show that brittle fracture is accompanied by an ejection of particles and
light. To date, measurements of these fractoemitted particles remain one of the few
methods for characterizing the dynamic processes leading to brittle fracture.
The Ultimate Strength of Materials
The ultimate strength of a material without pre-existing surface flaws that fails by
brittle fracture can be calculated based on a knowledge of its structure and the properties
of its interatomic bonds. This value is known as the cohesive strength of the material,
and it can be calculated for a material given the Young's modulus, the surface energies,
the atomic structure, and an assumption of the form of the interatomic potential.267
A cohesive strength calculation may use one of many pair potential functions.
Two possible potentials are the Morse potential and the Born potential. The Morse
potential is applicable to covalent systems and has the form 4(r) =4e 'r-r) -2e -(r-
where 42(r) is the pair potential energy function of r, r is the interatomic distance, ro is
the equilibrium interatomic distance, and &, and a are constants that can be found from
materials properties. The Born potential is applicable to ionic systems and has the form
02(r) =A +B where 4 and r are as described above, and A, B, and n are constants that
can be found from materials properties.
Cohesive strengths are derived from these pair potential functions by calculating
the derivatives of the functions and finding an expression for the maximum stress
encountered in breaking the chemical bonds. Stress is defined as .= 1 d (r) where a
is the stress, and A is the cross-sectional area covered by the bond. The maximum stress
is then found where i(r) is a maximum, or when d 2(r) is zero. Thus, the
theoretical cohesive strength is found by calculating the value of r for which d ,2(r)
is zero, then substituting this value of r into the expression for 1 d0(r)
The resulting expressions for the cohesive strength for the two potentials are, for
the Morse potential, % Yy and, for the Born potential, E[ 2 r) where y
is the surface energy, E is the Young's modulus, and n=-rE-
For vitreous silica, with E = 7.3 x 01O Pa, y = 4.56 J/m2, and r. = 1.62 x 10-10
m, the theoretical strength given by the Morse potential is 23 GPa. This compares
favorably with the highest strength of 14.3 GPa measured for treated silica fibers in
liquid nitrogen by Proctor et al.,' but is about two orders of magnitude greater than the
measured strength for bulk silica or untreated silica fibers.
The Effect of Flaws
Discrepancies between theoretical strengths and measured strengths are common
in brittle materials. This occurs because the presence of flaws in materials result in local
strains that are much greater than the average strain within the material under stress.
Plastic deformation does not readily occur in brittle materials, so the stresses in these
areas of locally high strain cannot be relieved. These regions of high stress cause failure
to occur in a material even when the average stress throughout is much less than the
theoretical cohesive strength. Because of this, flaws may be considered as stress
concentrators in the fracture of brittle materials.
Griffith' calculated the tensile stress, a,, which could be supported by a sample
of a material containing a crack of a given length, c. His calculation was based on the
assumption that a crack will propagate under an applied stress when the rate of decrease
in elastic strain energy during crack propagation is greater than the rate of increase in
surface energy resulting from the formation of new surfaces. Thus, for a sheet with an
edge crack of length c, the critical stress is given by' -_ 2Ey Similar expressions
have been derived for various types of crack geometries.
Regimes of Brittle Fracture
The limiting processes in the brittle fracture of glasses in reactive environments
are different in the three regimes categorized by Wiederhorn.' Fracture processes of
regimes I and II are limited by environmental effects and occur at stresses, or stress
intensities, well below the theoretical strength of the material, and crack growth
velocities are affected by species in the chemical environment. In type III fracture, the
material becomes unstable under the applied stress without regard to environmental
effects and the high crack growth velocities that result are limited only by the mechanical
properties of the material.
Stress-induced corrosion or slow crack growth corresponding to types I and II
fracture are important in limiting the strengths of materials in many applications.
1 r--. T -T
0.1 In moist nitrogen
In dry nitrogen
0.2 0.3 0.4 0.5 0.6 0.7 0.6
Stress Intensity, MPa'm-0.5
Figure 2-1. Schematic showing crack propagation velocities as a function of stress
intensity factor. Adapted from Michalske."
Michalske" studied the fracture surfaces created by slow crack growth to better
understand the mechanisms involved. In this type of fracture, environmental species
diffuse to the crack tip, weakening the chemical bonds at the tip and causing crack
growth at stress intensities less than o, but greater than some minimum stress. The
velocity of the crack tip is then controlled by different mechanisms in the two regions.
Since the reactivity of the material to corrosion at the crack tip is dependent on the local
stress, the rate of reaction and therefore the crack growth velocity is dependent on the
stress intensity at the crack tip in type I fracture. The rate of diffusion of the
environmental species to the crack tip controls the crack growth velocity and is nearly
independent of the stress intensity in type II fracture. Thus, the crack growth velocity
in region II has a lesser dependence on stress intensity than in region 1. A schematic
showing this behavior is given in Figure 2-1.
Crack Prooaeation Velocities
The maximum velocity of crack fronts in type III fracture is limited by the speed
at which information of the approaching crack front may travel through the system. In
a mechanical system, this speed is the speed of sound. Different methods have been
applied to calculate to maximum crack propagation speed in real materials, giving similar
expressions of the form2 v=C E where C is a constant between 0.4 and 0.6, and p
is the density of the material. Maximum crack velocities are found by experiment to be
consistently less than, but on the same order as, the values given by the above
Reading the Fracture Surface
Because brittle fracture is a sudden process that does not easily lend itself to real-
time study, much of the analysis of brittle fracture is performed by fractography, the
study of fracture surfaces. Since we seek to describe the processes occurring during
fracture, we need to understand what fracture surfaces can tell us about these processes.
Through fractography, one can determine the location of the flaw that initiated fracture,
determine the size of the initiating flaw, and determine the stress at which the sample
fractured by examining the patterns on the fracture surface.
Brittle fracture surfaces have four distinct regions--the smooth mirror region near
the initiating flaw, the rougher mist region beyond the mirror region followed by the
coarse hackle region, and finally the crack branching region extending beyond the hackle
to the limits of the fractured sample. By locating the smoother regions of a fracture
surface, one approaches the area of the initiating flaw, and by advancing to the center
of the smoothest region, the mirror region, one can pinpoint this flaw.
Fracture surfaces can provide information on the strain state of the material before
fracture as well as the location of the initiating flaw. The transitions between the mirror
and mist regions, the mist and hackle regions, and the hackle and crack-branching
regions are distinct. It has been found that the radii of these regions are related to the
failure stress by the relationship13 oaJ=M. where a, is the failure stress, r. is the
transition radius corresponding to either the mirror-mist, the mist-hackle, or the hackle-
crack-branching transitions, and M, is a material dependent constant for each of the
Recently efforts have been made to describe the aspects of the fracture surface
that have traditionally eluded a meaningful quantitative description. These aspects are
the seemingly random fluctuations that exist within all regions of fracture on varying
scales. Fractal geometry has been used to describe the statistical nature of these
fluctuations, as will be discussed later in this review.
Bonding in Silica
Silica is the chemical compound SiO2. Oxygen is a group VI element and silicon
is a group IV element. Oxygen is more electronegative than is silicon, but the difference
in electronegativities is not great enough to make the bonding in SiO, purely ionic, so
the bonding has significant covalent nature. Because silicon is a group IV element, it
requires four bonds to complete its valence if covalent bonding is assumed, or it has four
valence electrons to donate if ionic bonding is assumed. Oxygen is a group VI element,
so it requires two bonds to complete its valence if covalent bonding is assumed, or it will
gain two electrons to complete its valence if ionic bonding is assumed.
Structure of Silica
The basic structural unit of most condensed silica phases is the silica tetrahedron
as shown in Figure 2-2.14 In each silica tetrahedron, there is one silicon atom
tetrahedrally coordinated by four oxygen atoms. Each oxygen atom is then shared by
two such tetrahedra to maintain local charge neutrality. Because the silica tetrahedron
is the basic structural unit of silica glass as well as most crystalline phases of silica, even
the glass phase is not without short-range order on the scale of a silica tetrahedron.
Because of this, we refer to this phase as vitreous rather than amorphous.
Figure 2-2. The silica tetrahedron is the basic structural unit of silica structures.
a) One silica tetrahedron. b) Two tetrahedra linked by corer sharing.
0 is the Si-O-Si bond angle. Adapted from Holloway.'4
Phases of Silica
Because silica structures are relatively open, they are relatively free to change
symmetry with temperature. Thus, silica comes in several polymorphic phases. The
major crystalline phases of silica are quartz, the phase thermodynamically stable at
temperatures below 1140 K, tridymite, stable in the temperature range of 1140 K to 1770
K, cristobalite, stable between 1740 K and the melting point of 2000 K, and coesite, the
high pressure phase.'5.'6 In addition, there are further distinctions within the individual
phases. The quartz, tridymite, and cristobalite phases are subdivided into high and low
subphases. The transformations between the subphases are displacive, meaning that there
is no rearrangement of the structural elements necessary for transformation, only
reorientation of the structural units. The transformations between the major phases are
reconstructive, meaning that there are structural rearrangements required to go from one
phase to another. The vitreous phase is the usual result of cooling molten silica to a
solid when no effort is made to induce crystallization. For a summary of the phases of
silica, see the review article by T.D. Coyle."
In crystalline silica, the tetrahedral units are linked in an orderly fashion, resulting
in a structure possessing long-range order. In vitreous silica, the tetrahedra are linked
in an apparent disordered fashion. This results in a structure possessing short-range
order given by the rigid constraints on the size and shapes of the tetrahedral units, which
diminishes rapidly beyond the first coordination distance because of the less ordered
Characteristics of the Phases of Silica
Phase Symmetry Density (g/cm3) Tet Linkage
Low Quartz Hexagonal 2.6 Corner
Tridymite Hexagonal 2.22 Corner
Cristobalite Tetragonal-Cubic 2.33 Corer
Coesite Monoclinic 3.01 Corner
Keatite Tetragonal Corner
Shishovite Tetragonal 4.35 Octahedrally
Silica-W Tetragonal 1.97 Edge
Vitreous NA 2.20 Corner
linkage of the tetrahedral units. The characteristics of the phases of silica are given in
The Crystal Structure of Cristobalite
Cristobalite is of interest to us because its high symmetry, low density structure
is probably the most comparable to the vitreous phase. At room temperatures, both
phases are essentially thermally arrested high temperature phases. Thus cristobalite is
the crystal phase of choice to compare to vitreous silica in studies of the effects of
structural order and disorder in silica.
Although quartz may be the equilibrium phase, the cristobalite phase can be
relatively stable when quenched to low temperatures because the reconstructive
transformations between the major phases is slow. Within the cristobalite polymorph,
two sub-polymorphs are present, high, or beta, and low, or alpha, cristobalite, which
undergo a displacive transformation within the temperature range of 470 540 K. This
displacive transformation results in a symmetry change from the high temperature cubic
symmetry to the lower temperature tetragonal symmetry upon cooling.
The structure of high cristobalite"1 is similar to the diamond structure. The silicon
atoms are in the diamond positions and oxygen atoms are located between each pair of
silicon nearest neighbors. The symmetry of high cristobalite is cubic, and the structure
belongs to the Fd3m (Oh7) space group. The unit cell of high cristobalite has an edge
length of 7.16 A at 560 K and contains eight formula units.
The structure of low cristobalite has the same connectivity as the high cristobalite
structure. Although this phase has tetragonal symmetry with prism dimensions
ao = 4.96 A and co = 6.92 A, the structure is pseudocubic when the cell is rotated by
450 about the c-axis and the sides a. are transformed into diagonals of the base of the
prism, that is, ao = 2a, = 7.02A and co = co 6.92A. In this description, the
structure is a relatively small distortion of the high cristobalite structure. The unit cell
of low cristobalite contains four formula units and belongs to the space group P412,
The Rine Structure of Silica
A useful description of intermediate structure in the phases of silica is the ring-
size distribution. The statistical nature of this measure allows it to be used for vitreous
silica structures as well as crystalline structures. The ring structure plays prominently in
the description of silicate structures because these structures are relatively open as a
result of the low coordination of silicon and oxygen atoms, allowing identification of
individual rings in the network.
A ring in a silica structure is a closed loop of adjoining tetrahedra that are
identified by the number of tetrahedra contained in the ring rather than by the number
of atoms that are bonded to form the ring. Work by Galeener" has shown that different
numbers of tetrahedra in rings give rise to differences in the energies of the structure.
Works by Marians and Hobbs20,21 and Marians and Burdett22 have described the
topological constraints on the ring size distributions in crystalline and vitreous forms of
silica and on vitreous silica surfaces. Most crystalline forms of silica are composed
entirely of six-membered rings, whereas vitreous forms necessarily have some ring sizes
different from six.
The ability to make vitreous silica by a sol-gel route has enabled a detailed study
of the role of rings in the formation of a silica network structure. In sol-gel processing
of silica, silicon alkoxide precursors in an aqueous solution undergo condensation
reactions such as the following, where alkoxides are identified by the formula unit
Si(OR)4 + H20 (RO)3Si-OH + ROH
2 (RO),Si-OH 4 (RO)3Si-O-Si(OR)3 + H20
(RO)3 Si-OH + Si(OR)4 7 (RO)3Si-O-Si(OR)3 + ROH
As the polymerization process continues, a network structure gradually forms. Because
of the thermal and environmental constraints on the system, the resulting gel is not
generally structurally similar to bulk silica. The absence of a high temperature step in
the gelation process prevents the system from exploring different structural configurations
that would result in the system seeking a global metastable state similar to that obtained
by fused silica. Instead, the system structure is the result of a growth process by which
Si(OR)4 molecules react first to form chains, then to form a network structure as
SiO,(OR)(4 groups become linked with an associated loss of water molecules. Because
the interim structure results from the dynamic growth process emanating from reaction
centers, the structure tends to be fractallike--the geometry of the structure grows from
a center point and through a series of branches, attempts to fill space. Silica gels are
transformed into bulk silica structures by the processes of drying, stabilizing, and
densifying. During drying, liquid water and reaction by-products are removed from the
gel. In stabilization, physically adsorbed water is removed. Finally in densification,
chemically adsorbed water is removed, and the structure attains full connectivity. During
the densification step, structural changes occur that elucidate some of the medium-range
structural details of silica.
Raman studies of gels during the densification process tell an interesting story.
Superimposed on a somewhat broad Raman spectrum are narrow "defect" bands at 490
and 608 cm', so called because their presence cannot be explained by continuous random
network models.24 As gels are stabilized by driving off chemically combined water in
the structures and then densified by collapsing the porous gel structure, the intensity of
these defect bands increase, then decrease. This has permitted identification of these
bands with the condensation of dangling silanol groups on the surfaces of the gels. As
these dangling bonds become linked, highly strained three- and four- membered rings are
formed and subsequently destroyed as the gel densifies and the internal surfaces are
eliminated. This shows that the concept of rings corresponds to a real structural
component in silica glass.
Elastic Properties of Vitreous Silica
Mallinder and Proctor" found that vitreous silica at low temperatures displayed
a nonlinear elastic modulus in which the apparent elastic modulus is increased at high
strains. At 77 K, the Young's modulus, E, and the shear modulus, G, as functions of
strain, E, were found to be E=71.9(1 +5.75e) GPa and G=31.5(1+3.06e) GPa,
respectively. The authors attributed this behavior to the opening of Si-O-Si bond angles
in the structure, causing the elongation of silica rings. They reasoned that extension of
the structure could be more easily accomplished by initially opening the Si-O-Si bond
angle than by elongation of Si-O bonds.
The same authors also found that sodium silicate glasses exhibited a nonlinear
elastic modulus in which the apparent elastic modulus decreased with increasing strain
under similar conditions. In this case the authors reasoned that the extra cations were
positioned in the structure so as to prevent the rotation of tetrahedra under low strain.
The function for the elastic modulus given for this material was given as
Thomas and Brown26 studied the elastic properties of fused silica and found that
this material exhibited a similar nonlinear elastic modulus at room temperature. The
authors suggested that this anomaly may play an important role in crack propagation and
Proctor2 examined the effects of a nonlinear elastic modulus on the fracture
behavior of fused silica. The author concluded that the nonlinear elastic effect must be
considered when calculating either failure strains in high strength samples or stress values
from failure strains, but that the nonlinear effects on the calculations of the size of
fundamental Griffith flaws are self-cancelling and generally can be neglected.
Kulawansa et al.2 performed an atomic force microscopy (AFM) study on silica
fracture surfaces. The authors of this study analyzed the details of the fracture surface.
They explained that some of the surface features observed were the result of a crack
blunting mechanism caused by a hardening of the material under high strain, as they
claimed would be consistent with a nonlinear elastic modulus that increases with strain.
The Fracture Strength of Silica
Measurements of the strengths of silica fibers indicate that nearly flaw-free fibers
have been produced. Proctor et al.8 performed strength measurements on as-drawn fibers
and found a narrow distribution of material strengths, indicating that samples were nearly
flaw free. When the samples were fractured in liquid nitrogen, the authors obtained a
maximum strength of 14.3 GPa. Kurkjian and Paek9 performed an analysis of the
tensile strength statistics for silica fibers and found that most of the variability of the
apparent strengths of fibers could be attributed to the variability of fiber diameters along
the gauge length of the tensile test, giving a nearly single-valued strength of 5.75 GPa
with a variation of less than 1%.
In a seminar at the University of Florida Department of Materials Science and
Engineering in the Spring of 1990, Kurkjian noted that measurements suggested that fiber
strengths were not completely single valued and indicated the presence of flaw sizes on
the order of a few angstroms. A scanning tunnelling microscope (STM) was used to
examine the topography of as-drawn fiber surfaces to explore the nature of the possible
Proctor27 interpreted fracture strengths for silica fibers and found that average
stresses were well below the highest measured value. Using elastic modulus data
obtained in earlier work, he obtained values of flaw sizes of 3.6 nm for fused silica
fibers and 7 nm for soda-lime glass. He then made the assumption that the most
prevalent surface flaw is the largest broken ring of silica tetrahedra on the surface of the
fiber and thus deduced average ring sizes of 5 to 7 tetrahedra for fused silica and 8 to
12 tetrahedra for borosilicate glass. The corrected strengths obtained from the later
strength data ranged from 17.9 to 18.7 GPa for fused silica and 6.8 to 10.1 GPa for
Stress Corrosion of Silica
A model for the chemical reactions resulting in the stress corrosion of vitreous
silica was proposed by Michalske and Freiman3 and was further explored by Michalske
so d H^
% O O SO
0o oW 0 0 0
(a) (c) (c)
Figure 2-3. Chemical reaction proposed by Michalske and Freiman for
environmentally enhanced Si-O-Si bond rupture. The strained crack tip
bond (a) adsorbs a water molecule, (b) is cleaved by dissociative reaction,
and (c) is converted to surface silanol groups. Adapted from Michalske
and Bunker.30 In this model, the authors propose that a water molecule aligns itself with
a strained Si-O-Si bond, resulting in a weakening of the bond, and the water molecule
is consumed in a reaction that breaks the bond and creates a pair of silanol groups. By
this model, three properties may be common to all active corrosive species: (1) the
molecule has at least one lone (nonbonding) electron pair (Lewis base), (2) the molecule
possesses a labile proton (Bronsted acid), and (3) the distance between the acid and base
sites conforms with the Si-O bond distance (0.163 nm). A schematic of this process is
given in Figure 2-3 and shows three stages in bond cleavage.
Hench and West''32 studied the reaction paths in the fracture of Si-O bonds in
rings of various sizes using molecular orbital calculations. They found that when
particular bonds of unconstrained rings were strained, the energies required to break the
bonds were significantly lower than the expected energy of the Si-O bond. This behavior
resulted because unconstrained rings of various sizes, each having one strained Si-O
bond, underwent a ring contraction mechanism whereby another tetrahedron within the
ring approached the oxygen atom of the strained bond, five-coordinating the tetrahedron,
and weakening the initially strained Si-O bond. The authors found that 4-membered
rings were the weakest rings, with an activation energy for Si-O bond breakage of 320
kJ-mol', followed by 3-membered rings with an activation energy of 400 kJmol', and
5-membered rings with an activation energy of 430 kJ-mol"'. When water was
introduced into the system, the various sized rings showed different tendencies for
hydrolysis, the critical step in environmentally assisted fracture. In this case, the authors
found that 3-membered rings were the easiest to hydrolyze, with an activation energy of
30 kJImol'l, followed by 4-membered rings, with an activation energy of 120 kJImol'
and 5-membered rings with an activation energy of 160 kJ-mol-. In both the water-free
and water-present cases, the authors concluded that the ring structure of bulk silica
influences the fracture of Si-O bonds, and must be considered when studying this
Because water is an active component in the stress corrosion of silica, and since
it is prevalent throughout the atmosphere, it plays an important role in all fracture
experiments on silica except for those experiments in which it is purposely and carefully
excluded. Understanding the role of water is of vital importance when interpreting the
observations of any fracture experiment. The fracture of silica in the absence of water
is quite different from the fracture of silica in a water-free environment, as demonstrated
by a comparison of strengths measured by Proctor8 (14.1 GPa) in a water-free
environment at 77 K and by Kurkjian and Paek29 (5.75 GPa) in a laboratory environment
at room temperature.
All simulations to be presented in this work pertain solely to a water-free
Survey of Applications
A molecular dynamic (MD) simulation is, as its name implies, a technique that
models the movement (the dynamics) of particles (of molecules). As pointed out by
Hoover3 in an introductory article on the subject, interesting results can be obtained from
even the simplest of systems because of the many-body interactions that may be involved.
Molecular dynamic simulations are not difficult to visualize. Imagine the
collective motion of several balls on a billiard table with a frictionless surface and
perfectly elastic bumpers. Imagine now a computer model of this motion, and you have
envisioned an MD simulation of a two-dimensional hard-sphere system with rather
peculiar boundary conditions. Further, MD simulations are capable of extracting
thermodynamic quantities such as pressure and temperature from the dynamic system,
providing a simulation link between the system of discrete particles in motion and the
thermodynamic quantities that describe that system.
Molecular dynamic simulations were first used by Rahman4 to study the
correlations of atomic motions in liquid argon. MD was first used in simulations of
liquids for various reasons. Clearly, liquids are dominated by particle dynamics,
therefore, their properties can be quite well simulated by MD methods. Since liquids are
fluid, their structures are sufficiently dynamic so that slow nucleation and growth
processes are not required to obtain reasonable simulations of their structures, and
comparatively short MD simulations are capable of simulating reasonable structures.
Finally, despite the dynamic nature of liquid, there are strong interactions among atoms
and transient structures of sufficient interest to warrant a simulation of structure and
properties, as opposed to an ideal gas where there is little interaction between atoms and
The applications for MD simulations have expanded over time. Non-equilibrium
MD is concerned with the dynamics of systems and not in thermodynamic equilibrium,
and is used in simulating glass structures. The structure and dynamics of materials
surfaces have been studied. For example, Heyes et al.3435 and Heyes3" studied the
structure and dynamics of surface and near-surface regions in ionic crystal films.
Broughton and Gilmer37n."39 studied crystal-fluid interfaces of Lennard-Jones metals.
Gilmer, Grabow, and Bakker4 used MD simulations to model epitaxial growth in
Lennard-Jones systems. Briinger et al.41 reported the use of MD simulations in the
determination of crystallographic structures of protein molecules from X-ray and neutron
There has been recent interest in using MD to study the interactions between
surfaces in a field called "nanotribology." Landman and Luedtke42 reported MD
simulation results in studies of the interaction of metal tips with metal substrates.
Harrison et al.43 reported the use of MD in studies of friction between diamond surfaces.
Belak et al." used MD simulations to model atomic scale deformation of surfaces by a
process similar to machining.
The evolution of computer hardware, along with the evolution of computer codes,
is continually increasing the capabilities of MD simulations. MD simulation codes are
capable of being run efficiently on parallel computers with some modification.
Furthermore, researchers have developed specialized computers for performing certain
simulations more quickly and efficiently. An example of such a computer is described
by Bakker et al.45 A result of the development in computer hardware is that systems
composed of tens of millions of particles can now be simulated on parallel computers
consisting of thousands of processors, whereas systems having fewer than one thousand
particles were simulated in the infancy of MD.
Recently, MD simulations have reached popular appeal with the advent of
commercially available MD simulation software packages, such as the Insight package
from BioSym Inc. One can now buy a simulation system that is capable of simulating
a wide range of materials and configurations. It is likely that the use of MD in all types
of materials simulations will increase dramatically since it is no longer necessary to either
write simulation codes or modify existing simulation codes for a wide range of
applications. Despite the advent of commercially available simulation codes, many
researchers will most likely continue using custom-made codes, which are optimized for
the simulation of certain systems or to run on specific computer hardware.
Survey of Potential Functions
The identity of a simulation is strongly dependent on the potential function used,
which in turn is based on the type of bonding in the material being simulated. There are
several common potential functions that are used in MD simulations.
In order to simplify the description of the various potential functions in common
usage, we will use the same means of summation for all potentials, and common symbols
whenever possible to permit comparison between the potentials (see the Key to Symbols
in the Appendix). Summations are not necessarily listed as they would be performed the
most efficiently in the computer code, but are written for logical simplicity. The general
interaction potential for a system will be written as
where N is the number of atoms in the system,
i, j, and k represent individual atoms as specified in the summations,
u, v, and w represent the atom types of i, j, and k, respectively,
2l.(r0,) is the function representing the two body portion of the potential for types
u and v, and
-3,_(r,, r,, r) is the function representing the three body portion of the potential
for types u, v, and w. Since the three body portion of the potential is
most often a function of two bond lengths and a bond angle, it may
usually be equivalently referred to as .ik(rj, rj, OJ). We have written
the summations to work most directly with this form because it is
applicable among all the three-body potential functions that we will
The use of indices u, v, and w to denote atom types is a departure from the usual
notation in the literature. This notation is used to unambiguously differentiate the
constants or functions are dependent on atom types only from the variables that are
dependent on individual atoms. Although this is usually clear from context, there have
been some cases of ambiguities and errors in the literature. For systems of only one
atom type, the subscripts u, v, and w will not be used.
A common pair-potential used for the condensed phases of noble elements is the
Lennard-Jones 6-12 potential. The function for this potential is A2(r)=--_ A+ where
A and B are positive constants and r is the distance between atoms. This potential is
usually written in the dimensionally more appealing form (r)=4ef2(), where
f2(r) = r12 r' is a dimensionless function, and a = (BIA)"6 and e = A2/4B are
constants.46 This potential was also commonly used in the simulation of metals, since
it often results in a close-packed crystal structure similar to the structure of many metals.
Recently, a potential based on the embedded-atom-method (EAM) developed by
Daw and Baskes47 has become a popular method of simulating processes in metals with
empty or filled d bands.4" This potential uses the assumption that atoms are embedded
into a local electron density provided by the remaining atoms of the system. The
electron density function of the system is a many-body interaction problem, so this
potential is a many-body potential. This potential has an advantage in metals because it
is less empirical than is the Lennard-Jones 6-12 potential traditionally used in these
Covalently bonded materials require a potential with an angularly dependent term
to simulate the angularly dependent bonds that result from bonding through molecular
orbitals. One such potential was developed for silicon by Stillinger and Weber.49 An
angular dependence is obtained by adding an appropriate three-body term to a truncated
Lennard-Jones-type potential. The form of the potential is 02(r) =f' for the two-
body part, and (r, r 9~ )e.r' ( k o) for the three-body part. The function f in the
two-body part of the potential is
[A(Br P-r )exp( r
while the function f, in the three-body part of the potential is
Xexp (Y Y coso+ I2,
f0(r,r6k j rq
[0, r >a V r >a,
and the constants e and a serve to make2 and f dimensionless functions.
Ionic materials are generally modeled using a potential that incorporates a
coulombic potential with an exponential repulsive term. The Born-Mayer-Huggins
potential is such a potential. It has the general form
(.) +A=- e brc rj S+D rS
where r> is the inter-ionic distance, q, and q, are the ionic charges, and A,,, b, C,, and
D,, are constants.
Potentials for Silica
There are several commonly used potentials for modeling silica. The first of
these is a variant of the Born-Mayer-Huggins described by Busing,so first used in silica
by Woodcock et al." and later modified by Soules. The second is a variant of the Born-
Mayer-Huggins potential that incorporates a three-body term similar to that used in the
Stillinger-Weber potential. These two potentials will be described in detail in the
procedure section of Chapter 4.
Vashishta et al.52 derived a potential function that included a three-body part and
applied it to various silica polymorphs. The two-body portion of the potential is given
H qq (aZ,2 Z2)
S(r r. 2i
and the three-body portion of the potential is given by
B ,exp (-- -- coos ,- r)
0,._(r,rO = k) r-ro i -roj
0, r> ro V > r
with constants H,, and 7,, representing the strengths and exponents of the steric
repulsion, Z, and a. representing the effective charge and the electronic polarizability
of ions of type u, r4, representing the decay length of the r' interaction term, B,,, ro,
and 0,, representing the strength, the cutoff distance and the preferred bond angle for
the three-body interaction, and I an unspecified constant. The authors claim that this
potential gives the correct Si-O-Si bond angle distribution where earlier potentials failed.
Another potential function that is now commonly used for silica is described by
Vessal et al.53 This potential function has the most adjustable parameters of any
reviewed here, and claims an excellent fit to experimentally measured properties of
various phases of silica. Wright54 found this potential created a vitreous silica structure
having the best fit to neutron diffraction data of any MD simulated structure at the time
of the comparison. The authors also claim agreement with the presumed Si-O-Si bond
angle distribution in real silica. The functional form of this potential is quite complex.
The two-body portion of the potential has long-range and short-range parts. The long-
range Coulombic portion of the potential is calculated using the Ewald55 method applied
to a fully ionic system of silicon and oxygen. The short-range part of the potential is
given by the four-part Buckingham potential
A ..r +B l.u +Cl.r3 +Dr42? +El.r, +Fl.., r,,,< r < r2.
P.u pr+Qiuri2 +Rui +Sr, r2~
-- r3., < < r.
where r, is the interatomic distance between i andj, and most other constants are purely
empirical. The three-body potential for O-Si-O interactions has the form
U(rrk,Ok) = iA (B(GO) )'exp(- r )
where A, is a constant given by A k B(O,) is a function given by
B(O,) =(00-7)2-(,, -I)2, and the constants kk and 00 are the three-body spring constant
and the tetrahedral bond angle, respectively. With various spline approximations, etc.,
this potential has more than 40 parameters. With so many adjustable parameters, it is
not surprising that this potential function can claim the best fit to the experimental
properties of silica, although the simultaneous optimization of all parameters is
Alavi et al." proposed a potential that adjusts the effective charge of silicon and
oxygen ions based on their coordination. This potential was proposed to correct
problems created when network-modifying ions are included in the structure in which
some of the network forming ions are not fully coordinated. The concept of variable
charge may be important in fracture studies as well where coordination of atoms change
during fracture. Tsai et al." have used such a mechanism in the simulation of fracture
in silicon. Thus, this charge-transfer potential may be one to watch in the future.
Applications of MD in Models of Silica
Woodcock et al." first simulated a silica-like structure using a variant of the
Bom-Mayer-Huggins potential as used for ionic systems. These authors acknowledged
that a purely ionic potential was not fully applicable to a system with a fair degree of
covalent character, but they tried the potential since it was available and easy to apply.
Despite the known limitations of this potential, comparisons of simulation derived pair
correlation functions (PCF's) with PCF's deduced from x-ray diffraction data showed
rough agreement and a comparison of system energies showed that the potential energy
calculated from the simulation was 5% less than an experimentally derived Born-Haber
lattice energy for quartz. The authors made observations and performed experiments on
the simulated system. They measured diffusion constants at various temperatures and
pressures, and studied the permanent compaction of silica by the application of a
compressive stress. Mitra et al.58 performed simulations of vitreous silica structures
using a variant of the Bom-Mayer-Huggins potential and reported good agreement to x-
ray and neutron diffraction data. These authors also examined the internal energy as a
function of temperature and reported that the simulated system exhibited a glass transition
at 1500 K.
Soules59 performed a simulation of sodium silicate glasses using a potential similar
to Woodcock et al.5" Again, the structure was reported to be in significant agreement
to x-ray diffraction data. Soules found that alkali and alkaline earth cations clustered
around non-bridging oxygens, as is generally predicted for these structures. Coordination
information showed that the number of non-bridging oxygen atoms was strongly
dependent on the amount of alkali ion present. Ring size distribution data showed
changes with varying alkali content as well, although it was not clear from the limited
extent of the reported ring size distribution what the changes were. Soules and Busbey6
further performed simulations to study diffusion of sodium in alkali silicate glasses and
found that diffusion coefficients obtained from simulation at high temperatures were in
agreement with extrapolated experimentally measured values. They then concluded that
purely ionic potentials for silica were sufficient to model transport properties.
Garofalini and Levine," and Feuston and Garofalini62 examined the structural
changes on surfaces of vitreous silica. These simulations created arbitrary surfaces by
removing periodic boundary conditions between a pair of sides of the MD cell. The
general finding of these simulations was that less-highly charged ions tended to segregate
toward the surfaces. In pure vitreous silica, these species were oxygen ions, while in the
alkali silicates, alkali ions also tended to segregate toward the surfaces.
Garofalini" studied the reactivity of water on vitreous silica surfaces. He found
that water reacts with a silica surface through a series of penta-coordinated silicon
transition geometries in which water molecules are somewhat paired on the surfaces.
Comparisons of Simulated Structure to Experiment
Researchers generally do not like to point out the discrepancies of their models
in the published results of their MD simulations. They often compare radial distribution
functions derived from simulation with those deduced from diffraction experiments, and
if there are the same number of peaks in the distribution function and they occur in the
same general location, they call it a "good" fit. However, based on some minimal
constraints on the structure such as the preferred Si-O bond length for silica, in no case
can radial distributions be too far wrong. Wright"" has attacked the casual comparisons
of RDF's derived from MD simulated structures to the RDF's deduced from
experimental techniques such as neutron diffraction data and has proposed a standard for
comparing the experimental and simulated RDF's. This standard is a variant of the
crystallographic R factor, and is given by R = Tp( Ti (r) By this standard,
the current champion in the fit of simulated structure to experiment is Vessal et al." with
an R7 factor of 9.0% .
MD Simulations of Fracture
The use of molecular dynamics for studying the failure of materials is not new.
Two approaches have been taken in these studies. The first has been to perform fracture
experiments in simplified MD systems in order to perform detailed studies of the
elements of fracture. The second has been to study fracture using MD simulations that
attempt to model real systems.
Paskin et al.,6,"7 Dienes and Paskin," and Sieradzki et al.6 performed simulations
of crack growth in two-dimensional crystalline lattices. These works demonstrated the
applicability of the Griffith criteria in simple 2-D systems modeled by molecular
dynamics, and studied mechanisms of crack growth, including the formation of
dislocations and speed of crack propagation, as functions of applied stress, crack length,
and potential parameters. These works provide an atomistic basis to fracture mechanics,
however, they do not provide details on the fracture of particular materials.
MD simulations of fracture in silicate glasses were first performed by Soules and
Busbey70 on sodium silicate glass samples at elevated temperatures. The periodic
boundary conditions used made this experiment equivalent to the fracture of a fiber with
a diameter of about 50 A, so this study did not model fracture in a bulk material. The
essential characteristic of glass fibers in the real world that gives them interesting
mechanical properties is the absence of surface flaws, resulting in nearly theoretical
strengths, but the simulated "fibers" did not share this characteristic. At the atomic
scale, the surfaces appear rough and would not be considered pristine if the simulated
"fibers" were scaled to the size of real fibers. Thus, behavior similar to fibers is not to
be expected, and is not observed. What is observed is that samples undergo necking as
the simulated surfaces dominate their mechanical properties. Cavitation within the
samples is also observed, indicating that the system sometimes takes this route to absorb
the added volume while maintaining connectivity and minimizing the total energy at each
stage of strain.
Soules71 examined the role of the vitreous structure in the fracture of silica by
comparing fracture in vitreous silica to fracture in high cristobalite silica. He
investigated how the Si-O-Si bond must break during fracture and found that for the bond
to break under an applied strain, the bridging oxygen must vibrate out of its potential
well. Thus, he showed that fracture in silica is an activated process and undergoes a
relaxation. He also found that high cristobalite exhibited a high strength of 70 GPa as
compared to a strength of 24 GPa for vitreous silica.
Kieffer and Angell72 applied an isotropic tensile strain to a silica sample in order
to rupture the structure. They found that the resulting structures had a fractal character
similar to that of aerogels, and that the fractal dimension was dependent on density but
independent of strain rate, indicating that the fractal dimension for silica under a given
set of thermodynamic constraints is unique.
Ochoa and Simmons7 examined the role of strain rate in the fracture behavior of
vitreous silica. They found that within the range of 0.4 4.0 ps', the strength, as
defined by the maximum stress supported by the system during strain to failure,
increased as a function of strain rate from about 45 GPa to 95 GPa. They attributed the
decrease in strength with decreasing strain rate to a relaxation process where atomic
vibrations caused some reorientation of the structure capable of lowering its cohesive
strength. They also found that the effects of strain on the structure were dependent on
strain rate, with the highest strain rates causing an elongation of Si-O bonds and lower
strain rates resulting in an opening of Si-O-Si bond angles in the strained structures prior
to fracture. Finally, they discovered that the maximum stress did not occur at a strain
immediately prior to rupture of the structure, but rather at some point where some
fraction of the bonds were broken internal to the structure.
Ochoa et al.74 examined the role of thermal vibrations in the fracture of vitreous
silica and high cristobalite by simulating the fracture of these materials at strain rates
ranging from 0.2 8.0 ps-'. They used a potential function based on the Born-Mayer-
Huggins potential with two sets of potential parameters to investigate the sensitivity of
results on a specific potential. The authors found that the maximum stress during strain-
to-fracture increased with increasing strain rate for the vitreous silica samples, ranging
from 30 GPa to 65 GPa, but was nearly independent of strain for the cristobalite sample
at 80 GPa. They explained that the decrease in strength of the amorphous structure with
decreasing strain rate was a result of rotation of (SiO)4' tetrahedra during strain-to-
fracture at low strain rates that was not allowed at high strain rates. They found that
voids placed in the samples resulted in a decrease in the strengths of the samples and
altered the mechanism of fracture at low strain rates. Finally, they found that the
differences in the potential parameters did not significantly affect the behaviors observed.
Simmons et al.75 examined the mechanisms required for stabilization of silica
surfaces following fracture. They found that fracture is stabilized by the movement of
dangling silicon atoms back into the structure, leaving predominantly oxygens on the
surface, resulting in a repulsion between the two freshly created fracture surfaces.
Chaos in Fracture: Determinism and Fractal Surfaces
We earlier discussed brittle fracture in silica. That discussion was based on the
principles of fracture mechanics--maximum strengths of static structures, the effects of
flaws on the strength of materials, and the role of environmental species during fracture.
However, although fracture mechanics may serve to identify the causes of fracture, it is
ultimately unable to predict the path of fracture that results in the seemingly random
structure of fracture surfaces. The complexities of fracture surfaces indicate that the
particular dynamics of the fracture process is important in determining the path of
fracture. Although molecular dynamics simulations may be capable of elucidating the
role of dynamic processes during fracture, we cannot analyze every atomic motion in the
fracture of a simulated system. Because we can easily observe the deviation of the
system when we make a small perturbation in the starting conditions, the concept of
chaos with its sensitivity on initial conditions gives us a handle on analyzing the complex
behavior that cannot be analyzed on an atom-by-atom basis. Thus, without looking at
all dynamic interactions in detail, we can find the significance of these interactions in a
In nature, fractal-like shapes often result from chaotic processes. This is seen in
brittle fracture as well. Thus we will also discuss fractal geometry and its appearance
in brittle fracture.
Investigations of Chaotic Processes
"Chaotic" is a descriptive term applied to systems that follow paths which are
strongly dependent on initial conditions. The term, "strongly dependent on initial
conditions" implies that any deviation in initial conditions ultimately results in system
paths as divergent as system paths resulting from large deviations in initial conditions.
The concept of chaos originated from simulation work of weather patterns by
Lorenz.76,77 Lorenz found that small deviations in the starting conditions of a numerically
solved differential equation resulted in a divergence of the system behavior over time.
Although the paths of the systems remained close for some time, they gradually deviated
to the point that they no longer possessed any correlation that would indicate a shared
beginning. The paths taken by the simulated systems deviated as a result of the loss of
precision associated with converting the computer's internal representation of numbers
into a decimal representation of limited precision, then reconverting this representation
into an internal binary representation when re-starting the simulation. Although the
simulation was not necessarily capable of a fully accurate representation of the weather,
it demonstrated that without a total knowledge of starting conditions and external effects,
the prediction of weather for indefinite periods of time would be impossible. The general
implication of Lorenz's work is that the behavior of systems described by coupled
nonlinear equations cannot be predicted within a certain range in the long term. In an
MD simulation where a material described by coupled nonlinear potential functions,
similar behaviors may be expected.
Mathematical descriptions of chaos are particular. One of the tools used in
investigating the chaotic nature of a system is a Poincard map, which is a map of the path
of a system in a space of some or all of its variables. Chaotic systems are often found
to have maps that approach but never achieve a periodic behavior, and such a map was
constructed by Lorenz. Often a projection of these maps creates either a fractal or an
intricate pattern. Another method used in investigating the chaotic nature of a system
is an investigation of an apparently randomly varying quantity. In fractoemission
experiments, Langford et al.78 performed autocorrelation, Fourier transform, and box-
counting dimension analyses on the curves of fractoemission intensity as functions of
time. They found evidence of deterministic chaos from these analyses, and related it to
the fractal dimension of the resulting fracture surfaces as found by the slit-island
technique. We will discuss the implications of a fractional dimension in a real material
later in this review.
Predictability of Brittle Fracture
Fracture occurs when a mechanical system becomes unstable by the application
of strain. What does it mean when a system becomes unstable? A system is unstable
when it is in a state that will undergo a transformation to a favored state under any
perturbation. This means that as a system becomes unstable, even previously
insignificant dynamics within the system may play a role in determining the path to the
favored state. This seems to be the case in brittle fracture where the fracture mechanics
that can quite accurately predict when a system will become unstable is ultimately unable
to predict the precise path to fracture after it has become unstable.
The choices available in predicting the path of brittle fracture are quite limited.
Either brittle fracture either occurs by cleavage, where fracture surfaces approach atomic
smoothness, or it does not. In cleavage, the path to fracture is known or controlled. But
in all other fracture events, the path is not known and is not exactly controlled. For
example, distributions of sample strengths have not been completely explained by known
sample flaws as pointed out by Kurkjian previously. Also, according to a private
discussion with T.A. Michalske in 1991, there appears to be some uncertainty in the
dynamics of slow crack growth in pristine silica under vacuum where there should be no
environmental effects, and thus, no effect of stress on the rate of crack growth. The
question arises--is the inability to completely describe the fracture process the result of
this process being chaotic?
A fractal is a mathematical set of points that is best described by fractal geometry
and has a fractional dimension. Fractal geometry is a branch of mathematics that
encompasses Euclidean geometry as a special case--the case of integral dimensions.
Objects in geometry are idealized and are not realized in the physical world, and are
referred to as sets. Although a set is a mathematical abstraction, physical objects may
approximate sets. Thus, when a surface is described as being fractal, it really means that
the shape of the surface approximates a fractal set, just as when a physical object is
described as spherical, the description implies that the shape of the object approximates
a sphere to certain tolerances. Thus, physical objects may befractal-like, but they can
not be fractal.
One of the most fundamental descriptions of a set is its size. The measurement
of sets is made using mathematically rigorous techniques. In order to measure a set, the
proper dimension for the set must be used. The Hausdorff dimension for a set is the
dimension in which a meaningful measurement of the set can be made. As its definition
is beyond the scope of this discussion, refer to Falconer79 for a complete definition of the
The Hausdorff dimension is the only dimension for which the set may have a
finite non-zero measure. For example, in Euclidean geometry, the measure of an area
is infinite if the measurement is made using a measure of length, is zero if the
measurement is made using the measurement of volume, and can only be non-zero and
finite only if the measurement is made using the same dimension of the set, that is, using
a measurement of area. However, in contrast to examples of sets in Euclidean geometry,
it is not always obvious in what dimension a measurement of a non-Euclidean set should
be made. Thus, although the Hausdorff dimension is not in itself a measure, it is
essential to the measurement process, and it is a useful descriptor for a set.
The Hausdorff dimension is the most fundamental definition of the dimension of
a fractal set, but other definitions are also used. The box counting dimension operates
on the graph of a set, and is defined as dim F lim logN(F) where F is the set that
describes the object, the function N, gives the number of boxes of linear dimension 8
required to cover F.79 In many cases, the box-counting dimension is equivalent to the
Fractal sets are of interest to mathematicians because of their unusual properties,
but why are sets that are best described using fractal geometry of interest to physical
scientists? The answer is that shapes found in nature often approximate the shapes of
fractal sets. Since mathematicians have discovered the generating functions for creating
fractal sets, physical scientists may assume that the physical processes that result in the
formation of such shapes are described by similar functions.
Several authors have related how physical objects approximate fractal shapes.
The "proof" that they give to show that these objects are fractal is that the "fractal
dimension," a limited approximation to the box-counting dimension, is not an integer.
The definition of "fractal dimension" used by Mecholsky et al.80 is given by the relation
Nr = 1, where D is the "fractal dimension," N is the number of boxes of size
proportional to r that cover a set, or the number of ruler lengths proportional to r that
measure the length of a curve. These are related to quantities in the definition of the
box-counting dimension where dimB F -D, N, N, and 5 r. An important distinction
is that no limit is involved in the definition of the "fractal dimension" as described here.
Thus, the quantity D is a useful description of an object that approximates a fractal over
a limited range of r. However, it does not prove that the object is fractal, only that it
possess some properties that approximate those of a fractal object over a limited range
Although the Hausdorff dimension remains the most precise definition of the
dimension of a fractal set, the term "fractal dimension" is often used in the physical
sciences to denote the approximate box-counting dimension as defined above. Because
of this, we will follow the same usage for the remainder of this work. However, one
should bear in mind that this fractal dimension is obtained only over a limited range of
scale, and the validity of extrapolations to different size scales must be tested.
Applications of Fractal Geometry to Fracture
Recently, there has been a great deal of interest in describing fracture surfaces
using fractal geometry. There are several reasons for doing this.
The first is that all statistically significant details of these seemingly random
surfaces can be described with a limited amount of information, because the same surface
features may be found over a wide range of scales, with the atomic scale being the
absolute lowest limit. Thus, the description of fracture surfaces using fractal geometry
provides a possible link between those observations made by fractography and those
made by atomic scale simulation.
The second is that a scalar quantity, the fractal dimension, is obtained for the
description of fracture surfaces. Just as the pair correlation function is a useful
description of a random glass structure, the fractal dimension is a description of the
random structure of a fracture surface. However, a fractal dimension, like a pair
correlation function, does not unambiguously describe the structure. It is merely a
statistic of the structure that may be compared to test the possible similarities of various
structures. Despite its limitations, the fractal dimension can be related to material
fracture properties, as will be shown.
Finally, and probably most important from the standpoint of understanding
fracture, a description of a fracture surface as fractal-like may help in the understanding
of fracture processes. Fractal sets can be created by iterated generator functions. By
describing a fracture surface as fractal-like, one may be able to find the generating
function for the surface, and then elucidate some of the processes which occur during the
progression of fracture. Since fractal-like shapes often imply a chaotic generation
process, do the fractal characteristics of fracture surfaces imply that fracture is chaotic?
Later in this work, we will attempt to show how fracture may proceed to form fractal-
The Fractal-like Nature of Fracture Surfaces
If fracture occurred only by the orderly cleavage of crystal planes in crystals and
it always created minimum energy surfaces in glasses consistent with the applied
mechanical constraints, then the science of fracture would be nearly complete. However,
brittle fracture is not that simple. During the formation of fracture surfaces, the effects
of the dynamic lattice seem to counter the attempts of the material to minimize its surface
energy. Even in the fracture of a highly ordered material such as single-crystal silicon,
Tsai and Mecholsky" found that the resulting fracture surfaces were not smooth, but had
a fractal-like character. Mecholsky et al.0O reported that fracture surfaces of ceramic
materials also exhibit fractal-like characteristics from investigations of alumina and glass
ceramic fracture surfaces. Such surfaces indicate that a dynamic process is at work in
determining the path that fracture follows in brittle materials.
Because the description of fracture surfaces as being fractal-like implies that
fracture surfaces may have detail to nearly the atomic scale, such a description may
provide a link between macroscopic measurements and microscopic or atomic scale
behavior. Thus, the finding that fracture surfaces have a fractal-like character opens new
avenues for study into brittle fracture.
Mecholsky et al.82 reported a relationship between the toughness of a material and
the measured fractal dimension of a fracture surface of that material. They found that
the toughness was related to the fractional part of the fractal dimension by the relation
Kc =E/oD where Kc is the fracture toughness, E is the Young's modulus, D" is the
fractional part of the fractal dimension, and ao is a constant with units of length, which
may be representative of some characteristic length in the material or of the process zone
during fracture. Since surfaces with a higher fractal dimension generally appear more
tortuous, an easy way to remember this relationship is the phrase, "the rougher, the
Applicability of Fractal Analyses to MD Simulated Systems
An important question that needs to be answered is the capability of modeling
fractal-like structures by MD simulations. Quality fractal analysis requires a range of
scale that is not readily available to MD simulations. Can one derive any valid results
from an analysis of fracture surfaces formed in MD simulations.
Kieffer and Angell7 reported on the isotropic fracture of silica. Under these
conditions, a bulk silica sample appeared to attain a structure reminiscent of silica
aerogels. By using a box-counting technique applied to pair-correlation functions, they
found the fractal dimension of such samples to be less than two.
In meetings with John J. Mecholsky and Joseph H. Simmons, we discussed what
would be necessary to perform an analysis of fractal dimension on surfaces created by
uniaxial fracture. Mecholsky used a slit-island technique that is essentially a box-
counting applied to a surface on real fracture surfaces. We looked at Kieffer's paper as
a model to see how a fractal analysis might be accomplished on the fracture surface
created by an MD simulation. Kieffer's method viewed the fractured silica as a mass
fractal, whereas the method used by Mecholsky examined surface fractal characteristics.
In a mass fractal analysis, an object may have a dimension ranging from three,
corresponding to a continuous three-dimensional network solid, through two,
corresponding to a layer structure, to one, corresponding to a linear molecule. All
samples are composed of a finite number of atoms that can be arranged in three-, two-,
or one-dimensional arrays. Since scattering or pair correlation functions treat atoms as
point particles, this analysis has no bias towards finding a set to be three dimensional,
and the dimension found is wholly dependent on the arrangement of particles in three-,
two-, or one-dimensional structures.
Clearly, this analysis was not suited for finding the surface fractal dimension of
a solid as 2 + D' where 0 < D' < 1. Such an analysis would be best performed by
some kind of slit-island analysis. There are two problems associated this type of
technique. The first is that it is not a simple matter to define the fracture surface in a
MD system because the system is logically a system of point particles. The second is
that our samples are too small to perform a fractal analysis, because samples of size on
the order of tens of Angstroms and composed of units on the order of Angstroms could
not have the range of scale required to perform an accurate analysis of fractal dimension.
Thus, although a slit-island analysis could be performed in an MD system, it might not
give an accurate representation of the true fractal character of the surface.
Y.L. Tsai et al." recently conducted a molecular dynamics study of the fracture
of single crystal silicon on the atomic scale. This system was chosen because of previous
measurements of the fractal dimension on fracture surfaces of this material using a
scanning tunneling microscope. Also, this system had already been simulated using
molecular dynamics, and the potential function (Stillinger-Weber) was similar to one that
was already written into code (Feuston-Garofalini). The authors were successful in
performing a slit-island analysis of a fracture surface simulated by MD, and were able
to calculate a fractal dimension, although the range of measurement on which this
dimension was based was necessarily limited.
The technique used by Tsai differs from the technique used by Kieffer in an
important respect. Tsai's method is dependent on the space occupied by the constituent
atoms of a material, while Kieffer's method is dependent only on the spatial distribution
of atom centers. Tsai assumes that the constituent atoms can be described as spheres of
a convenient radius while Kieffer assumes that the constituent atoms are point particles.
This difference in approach is due to the difference between surface and mass fractals.
An important effect that results from this difference is that the method that used by Tsai
is dependent on the shape of the space chosen to envelope the atom center, or simply,
the shape of the atom. Since atoms do not have unambiguous shapes, the possible effects
of the shape chosen on the dimension measured is worrisome. On the other hand, the
creation of a surface by carefully chosen atom shapes may compare well to the surface
that may be measured by a scanning tunneling microscope or an atomic force microscope
if such techniques were truly capable of single atom resolution. A fractal dimension
obtained by this technique would then be comparable to one obtained from a microscopic
analysis on real fracture surfaces. Results of that study show a strong correlation
between fractal dimensions and surface roughness calculated from the MD simulations
and those calculated from experiments for 100, 110, and 111 surfaces.
Because of the problems associated with finding the fractal dimension of a
fracture surface created by MD simulations of fracture, we will not report any such
analysis in this work. We will, however, try to relate the atomic processes we observe
to possible generating functions that may result in the formation of fracture surfaces
having a fractal-like character.
Usine Molecular Dynamics to Model Fracture
In this work, we conduct a series of simulations of brittle fracture in a simple
"model" material for which a wealth of experimental measurements is available.
Although limitations of the modeled interatomic potentials and the spatial and temporal
extent of the simulated systems make it impossible to supply an exact match of the
experimental strength of the material or to calculate exact crack growth velocities, the
simulations were performed to investigate the underlying mechanisms for the trends
observed in fracture.
One of the biggest challenges in modeling fracture on the atomic scale in real
materials, and especially in glasses, is sorting out all the interacting processes that occur
during fracture. Logically, a material may be treated as a large system of interacting
particles. The dynamics of such a system are difficult to predict, and the cause and
effect relationships are difficult to describe analytically. Often the symmetry of such
systems can simplify the analysis, but this approach can not been applied to problems of
fracture where the order of systems is broken. In some cases, the modeling of fracture
in symmetric systems can be treated analytically, as in modeling cleavage fracture where
the system symmetry is dominant over the system dynamics. When a system has little
structural order at the outset, as in the case of a glass, an analytic treatment fails, and
a numerical treatment must be used.
Authors who present results from MD simulations often try to directly relate their
results to experimentally obtained results. An example of this is the direct comparison
of the system pressure as a function strain in MD systems to stress-strain curves of real
materials undergoing fracture. However, molecular dynamics gives a perspective of
materials fracture on a microscopic scale that cannot be directly compared to materials
fracture on a macroscopic scale. Therefore, the results of MD simulations are best used
to build an understanding of the processes that occur on the microscopic scale. With an
understanding of microscopic fracture processes and how such microscopic processes
might effect fracture on the macroscopic, one may be able to predict what the
macroscopic experimentally observed behavior will be, possibly to explain previously
unexplainable behaviors. It is not reasonable to expect that MD simulations will give
results which will directly relate to the macroscopic.
The potential function gives the simulated material its molecular identity.
Boundary conditions, system size, system density (for constant volume simulations), and
temperature serve to describe the equilibrium state of a system. In non-equilibrium
molecular dynamics, the path by which the structure is derived is also significant, and
needs to be studied. The study of fracture using MD then breaks down into a simulation
of the sample structure followed by the simulation of strain application and examination
of the resulting processes.
Modeling of physical systems is possible because many of the interactions within
a physical system can be ignored and fairly accurate predictions can still be made about
certain aspects of the system. Molecular dynamics modeling is possible in our systems,
for example, because the collective behavior of a large group of atoms is relatively
independent of the motion of electrons about individual atoms. In an MD system, the
processes that one intends to model are atomic motions that are controlled by interatomic
potential functions. Because MD approximates electronic bonding structure with
interatomic potentials and ignores electronic states, it cannot be applied to study the
formation of electronic defects in materials, but it has shown great success in studies of
multiple particle dynamics.
Although the fracture of a bond begins with a restructuring of the electronic
configurations forming the bond, the role of nearest neighbors becomes more significant
as the bond is stretched and, as will be shown, dominates the bond fracture process. The
complications caused by environmental species on bonding are eliminated by examining
fracture in a vacuum. In this case, again as will be shown later, MD calculations yield
a bond fracture process similar to the molecular orbital models proposed by West and
Hench." The MD simulations further benefit from the ability to examine properly the
effect of the broken bond on the behavior of the ensemble of neighboring atoms. A
problem not examined here is that of environmentally driven fracture, because of the
major role played by the chemical interaction of environmental species (H' and OH in
water) on the strained bond. In this latter case, ab-initio approaches that examine the
details of the electronic configurations of the bond are essential, as shown by West and
It is recognized that the forces between atoms are the result of the interaction of
electrons about the individual atoms, however, good results can be obtained by assuming
that the electrons are in some sort of stationary configuration, and that interatomic forces
can be modeled from some pseudo-potential. Systems of atoms with spherically
symmetric electron configurations are well modeled by a two-body potential, which
represents purely ionic character. Systems of atoms with non-spherically symmetric
electron configurations are better modeled by a three-body potential. A three body
potential that models a covalent-like electron configuration thus attempts to model a
In the case of silica, the potential functions that are commonly used either model
silica as a purely ionic system (Woodcock), a pseudo-ionic system (Soules), an ionic
system with a pseudo-covalent character (Feuston-Garofalini), or a pseudo ionic-covalent
system (Vessal). All these potential functions have limitations.
For our simulations of fracture, we have picked two potential functions that have
demonstrated good fits to various pieces of experimental data for silica, and yet are quite
simple. These are the truncated Bom-Mayer-Huggins potential developed by Soules83
and the empirical three-body potential for vitreous silica developed by Feuston and
Garofalini4 that includes an angularly dependent term. The Soules potential is purely
pseudo-ionic in nature, while the Feuston-Garofalini potential incorporates a pseudo-
The Born-Mayer-Huggins Potential
The basis of both functions is the Born-Mayer-Huggins potential. This potential
function has the form
r) P. b.exp ) + Cr-' + Dr-',
where u, v signify atom types,
0,, is the potential between atoms of type u and v as a function of r,
r is the separation distance between atomic centers,
q,, q, are charges on atoms of types u and v respectively,
e e, are radius repulsive parameters,
b is a positive constant,
p is a hardness parameter,
C,,, D. are the coefficients of dipole induced dipole and dipole -
induced quadrupole terms, respectively,
.,, is the Pauling coefficient, given by
0 = (1 + -U + -),
where Z,, Z, are the units of charge of type u and v, respectively, and
n., n, are the numbers of valence electrons of u and v, respectively.
The multiple terms become important when the atoms are highly polarizable and are
separated by small distances, but are generally neglected in the simulation of silica glass.
Some have questioned the accuracy of this simplification since oxygen exhibits a high
degree of polarization and the atoms are closely spaced in condensed phases. The
justification for ignoring these terms is that atoms within a solid are generally constrained
to the point that these terms would have minimal effect, and these terms involve many-
body interactions that are computationally difficult. When the multiple terms are
neglected, the potential function simplifies to
42(r) = .exp(-- --)+ .
By collecting terms, this becomes
2, (r)A exp=(-) qiq
A_, -= b expi .
From here, the Soules and the Feuston-Garofalini potential functions take divergent
Truncation or Screening Functions
In a simulation, it is desirable to have a potential of limited range in order to
minimize the number of particle interactions. However, the Coulombic portion of the
Born-Mayer-Huggins potential approaches zero slowly as the reciprocal of distance and
many neighbors would be required for calculations of all potentially significant
interactions. To overcome this problem, a screening or truncation function is used to
limit these interactions. This truncation process simulates the result of a lattice of
alternating charges as used in the calculation of a Madelung constant.
One method of dealing with the slowly convergent Coulomb potential is to
perform an Ewald sum." If one knows the symmetry of a crystal, the electrostatic
potential of ions in the crystal can be more quickly summed by performing the sum in
Fourier space. This is the basis of the Ewald sum approach. From the Ewald sum, a
screening factor is derived that breaks up the Coulombic part of the potential into a local
part that is summed over individual ions, and thus contributes to inter-particle forces, and
a background part that is just added to the overall potential energy and has no effect on
inter-particle forces. Strictly speaking, the results of an Ewald summation are only valid
for the crystal system for which they were derived, however, the functional form of the
screening factor is useful even when the inverse lattice sum has not been performed for
a particular system. Thus, the Ewald convergence parameter may be used in the
summation of the system potential energy using the modified Bom-Mayer-Huggins
2(r) A_ ,exp(] r q. erfc (,
where erfc() is the error function complement function,
e is a constant, typically 0.175 0.35,6 and
L is a linear dimension of the molecular dynamics cell or r,.
Parameters Used in the Soules Force and Potential Functions
Potential Function Parameters
p = 0.290023 A
Asisi = 3300 p ergs
As,-0 = 3000 p ergs
AO-O = 2300 p ergs
qs, = +4e
qo = -2e
where e is the electron charge
Truncation Function Parameters
r. = 5.5 A
The term eL is termed the convergence parameter. It is dependent on an arbitrarily
chosen constant and the size of the molecular dynamics cell, but it is independent of the
type of interaction.
The Soules Potential Function
The Soules force function has been shown to simulate vitreous silica structures
with approximately the correct density and with nearest neighbor distances in significant
agreement with neutron diffraction data. It is the product of the force function derived
from the Born-Mayer-Huggins potential and an empirical truncation function. A direct
differentiation of the Born-Mayer-Huggins potential gives
Fd-, (r) A~ exp )
dr p p r2
-100 1 1 1L
0 1 2 3 4 5 6
Figure 3-1. Plot of the Soules potential function for various interactions.
The truncation factor has the form T(r) = [ (r )]. The Soules force function is
Fo...v(r) = F sv(r) T(r)
S exp(-) + -]. -(
The potential function that corresponds to the Soules truncated Born-Mayer-
Huggins force function is
'Sol = 2 Fs^'oe2uv(r) dr
--A J exp(- 1.,]I-(r
i- exp(-'+ t Q -(rr-[ dr,
which is not analytically integrable. Thus, the potential energy of the system is
calculated using a numeric integration of this expression. The parameters used in the
Soules force function are given in Table 3-1. A plot of the Soules potential function is
given in figure 3-1.
The Feuston-Garofalini Potential Function
The three-body potential developed by Feuston and Garofalini is a sum of the
modified Born-Mayer-Huggins potential an empirical convergence parameter and the
three-body portion of the Stillinger-Weber potential.4 The two-body portion of the
2.(r) = A. exp(- + erfc( )
Notice that the eL term of the convergence parameter has been replaced by /,. This
makes the screening function independent of cell size, but dependent on the type of
200- -- -
', ', ~ Si-OSio' -- .. ..
0 1 2 3 4 5 6
Figure 3-2. Plot of the two-body portion of the Feuston-Garofalini potential function
for various interactions.
interaction, and adds another set of adjustable parameters. The three body portion of the
h exp[ + osOkCcsO < rf,. A k< r
46.U(r,,k(r -) -r2-) Q .) -rr<)
0, r,_ r V Tjk -> rt
The parameters used in this potential function are given in Table 3-2 and a plot of the
two-body portion of this potential is given in figure 3-2. This shows a shallower
Parameters used in the Feuston-Garofalini Potential Function
Two-body Potential Parameters
p = 0.290023 A
As,.s, = 1880 p ergs
As,,o = 2960 p ergs
AO-0 = 720 p ergs
qsi = +4e
qo = -2e
where e is the electron charge
3s,i = 2.29
si-o = 2.6
o-o = 2.34
Nonzero Three-body Potential Parameters
Xo-s-o = 180 p ergs
Xsi-o-si = 3 p ergs
mo-si-o = 2.6 A
Ys,-o-s, = 2.0 A
rso = 3.0 A
rsi-o.sC = 2.6 A
cos 0'o-si-o = -1/3
cos 's,.o-s, = -1/3
potential well in the Si-O interaction and a weaker repulsion in the Si-Si and 0-0
interactions than is obtained for the Soules potential.
Molecular dynamics is a computationally intensive technique in which the
simultaneous integration of the motion of many individual particles is performed over
short time steps.
Integration of functions may be performed in many different ways, but only a
couple of different techniques are commonly used when integrating the equations of
motion of a molecular dynamics system because they are well known to provide
reasonably good accuracy with a minimum amount of computation. These are the Verlet
"leapfrog" algorithm and the Gear Predictor-Corrector algorithm.
Our simulations generally made use the Verlet algorithm,5 given by the vector
r,,,. 2. r,, 2 At + FK ,
mi i (j k) i
where r,,,, is the position of the particle after the time step,
r,, is the position of the particle before the time step,
r,,,, is the position of the particle before the last time step,
t is a particular time,
At is the length of a time step,
m is the mass of the particle,
F5 is the sum of all two-body force components on the particle.
F is the sum of all three-body force components on the particle,
Fifth-order Gear predictor-corrector algorithm
The fifth-order Gear predictor-corrector algorithm, as encoded by J.M. Haile,8
was also tested. This was done to verify that the Verlet algorithm, known for its speed
and stability, but not for its overall precision, gives reasonable results as compared to the
slower but more precise Gear algorithm.
The Time Step
The time step is the basic unit of time in a molecular dynamics simulation. The
time step is an approximation for an infinitesimal amount of time that is valid for the
processes under observation. During this fundamental time step, in the Verlet algorithm,
the acceleration is assumed constant. In the Gear predictor-corrector method,
approximations are developed that move the constant function to a higher order time
derivative. In the selection of the fundamental time-step, it is desirable to use a value
short enough to maintain a constant system energy, yet a value long enough to observe
a physically significant process. The usual test for deciding if a time step is of
sufficiently short length is to perform a simulation over about a thousand steps and check
to see if the system energy remains fairly constant. If not, then the time step is too long.
A time step of 0.5 X 10 s was found to be of acceptable length, and was used in these
simulations. This corresponds to about 1/40 of an atomic vibration period.
Limiting Particle Interactions
The magnitude of the task of calculating particle interactions is reduced by
limiting the spatial extent over which these interactions are considered. Neighbors lists
for all particles that include all neighbors which may possibly interact with a central
particle are compiled. Only these neighbors are then considered when calculating
interactions with the central particle. The neighbors lists are updated periodically by
finding all particles within a distance given by the maximum interaction distance plus the
maximum distance a particle may travel during the period between updates. The period
between updates in these simulations is ten time-steps. The maximum distance, found
by Ochoa,8 that a particle may move in a solid in ten time steps at 0.5 X 10"5 s each is
0.5 A. Thus, the maximum distance for inclusion in a neighbors list from the central
atom is 0.5 A greater than the cutoff distance for interaction between particles, or 6.0
A in the case of a interaction cutoff at 5.5 A.
The Pythagorean theorem is used to find the squares of inter-particle distances.
Finding inter-particle distances requires square root operations that are computationally
time consuming, and are avoided. The force table is tabulated in terms of the square of
inter-particle distances, so inter-particle distances are required only when pair correlation
functions or bond angle distribution functions are accumulated.
S* 0* 0 0* 0
r .. m0 ; .
0000 0 *0 * *0 00
SPr imary cel
Repea t ed second a r y c e lls
Ne ghborho od abou t t om
Figure 3-3. The effect of periodic boundary conditions in a simulated system. Note
that the neighborhood of an atom may extend beyond the primary cell.
Application of Periodic Boundary Conditions
The use of periodic boundary conditions allows one to model a bulk-like material
with a limited number of particles. A schematic showing the effect of periodic boundary
conditions is presented in figure 2-3.
Periodic boundary conditions are applied to the simulated system when finding
inter-particle distances. If the difference in any coordinate is greater than one-half the
corresponding dimension of the cell, then the corresponding image of the particle located
in an adjacent cell will be closer to the central particle than the particle in the primary
cell. The difference in the coordinate is then adjusted to reflect this by either adding or
subtracting the corresponding dimension of the cell to minimize the absolute values of
the difference. It is the resulting difference that is used in the Pythagorean relation to
find the square of the distance between the particles. The algorithm by which this is
performed is as follows: The differences in coordinates are calculated from
d,i= x, x- ,
and are subsequently adjusted with
S= d., s. round(),
and used in the Pythagorean relation
where u = 1, 2, 3 denote the Cartesian dimensions of space,
x,, is the u Cartesian coordinate of atom i,
x,5 is the u Cartesian coordinate of atom j,
di is the u coordinate difference between atoms i and j,
d',, is the u coordinate difference between atoms i andj after correction
for periodic boundary conditions,
s, is the u dimension of the cell,
round(x) is a function that returns the integer nearest to its argument, x,
r, is the distance between atoms i andj.
The application of periodic boundary conditions gets more difficult when the cell
is sheared. Suppose that one puts the cell oriented with the coordinate axes into simple
shear as shown in figure 2-3. The opposite sides of the cell are shifted with respect to
each other. This must be taken into account when applying periodic boundary conditions
across the shifted faces. The periodic boundary condition equations are
c = round( -)
d'l d,, ,, c,) 7.(sC)
d" = d s. round(
S-TAl y Sy
---- --------- 0 1----*------ -* -,
/ / / /
*. e-.. ;*. . "
I M - -/
S 0 / *$ e 0
0 ** ** /
Pr ary ce ll
Logical vin ew of slheared primary cell
Zanlin Pri may celons
Logical v ew of sheared secondary cells
Figure 3-4. Schematic showing the application of periodic boundary conditions to a
system under shear strain.
By this method, full periodic boundary conditions are maintained even when the virtual
shape of the cell changes.
Handling Particle Interactions
Interactions between particles are handled by the following process--interaction
geometries are found, interaction forces are found, components of interaction forces are
summed on all participating particles, and particle positions are updated to reflect the
effect of interactions.
In the case of two-body potentials, the interaction geometry is simply the square
of the distance between pairs of central atoms and neighbors. The interaction force is
found from the force table, which is subsequently broken into components and summed
on each particle.
When the three-body component of the force function is applied, the interaction
geometry becomes the distances between the central atom and two of its coordinating
neighbors of a certain type, and the angle between the two "bonds." Since the number
of interaction variables is large, but the actual number of interactions is rather small, a
force table for a three body potential would be difficult to compile, but fortunately, is
not required. The force is calculated for all three body interactions as they are needed,
and components are summed on each interacting particle.
As a particle moves, it may leave the primary cell. When this occurs, the particle
is translated back into the primary cell by subtracting the appropriate cell dimension and
adjusting for shear, as required. Thus particles that leave through one side of the cell
reappear through the opposite side of the cell, perhaps shifted if the cell has undergone
shear. This allows one to control the dimensions of the cell while permitting the atoms
to move as their environments dictate.
Observing the System
As a simulation progresses, one can collect information to find out what the
system is doing. Any aspect of a simulation can be recorded. A few that are commonly
recorded are instantaneous or average thermodynamic properties, particle positions and
velocities, pair correlation functions, and bond angle distribution functions.
Thermodynamic properties, such as system temperature, energy, and pressure are
calculated. A running average of these values may be kept, so that outputs of either
instantaneous thermodynamic properties or average properties may be made.
The temperature of the system can be found by assuming that the system behaves
as a classical statistical system, with average particle kinetic energies given by the
equipartition theorem, EK = k T. The temperature can then be adjusted by scaling the
velocities of all particles to give the proper average kinetic energy.
The system pressure or stress tensor is calculated using the Virial of Clausius,
which has the form
P +(NknT rI ENrVw(ro])
When converted for use in simulations, it takes the form
P O 3 NN