Download Comparison of SiO2 Molecular Dynamics Simulations with Experiments and more Exercises Dynamics in PDF only on Docsity!
American Mineralogist, Volume 73, pages (^) 941-955, 1988
Molecular dynamicssimulationsof SiO2 melt and glass:Ionic and covalentmodels
J. D. Kusrcxlo* A. C. Lls.lcl
Department of Geology and Geophysics,Yale University, New Haven, Connecticut 065 I l, U.S.A.
Ansrru.cr Computer simulations using the molecular (^) dynamics (MD) techniquehave (^) been carried out on SiO, melt and glass (^) ranging from 300 to 9000 K and 0 to 175 kbar. The MD simulations used two sets of interatomic potentials-a simple ionic model with Born- Mayer repulsion terms and a central-force-fie.ldcovalent potential-both of which were derived from quantum mechanical calculations on molecular clusters. Thermodynamic, structural, and diffirsion data were obtained from these model systemsand compared to experimental values (^) wherever possible. The MD resultsindicate that earlier X-ray diffraction data may need to be re-interpreted with respect to the SiOSi angle distribution. BondJength and bond-angle responsesto pressureand temperature changescompare favorably with experimental and theoretical studies on the d-qvartz structure. Ring-distribution analysesshow that planar rings con- taining three silica tetrahedra are present in the simulated glass.Predictions of second- order thermodynamic properties (^) of the system (C, a, B) reveal inadequaciesin the ionic approximation when applied to vitreous silica. Arrhenius plots of ln D vs. temperature and pressurebetween 2000 and 6000 K and between 0 and 200 kbar show that at 6000 K, D increaseswith pressurefrom 1.3 x (^) l0- at 40 kbar to 6.1 x^ l0 (^5) cm2ls at 100 kbar, then decreasesto 1.5 x (^) l0-5 at 175 kbar. Detailed analysisof the diffusion mechanismat pressureslessthan 100 kbar indicates that it is defect-controlled with high correlations between the percentagesof 3- and S-fold silicons, nonbridging oxygens,and the difusion rate.
Ixrnonucrrox In order to understanddiffusion, nucleation, and crys- tal-growth kinetics in molten systems in an integrated way, an atomistic picture of the melt structure and its dynamics is necessary.Spectroscopicstudies of silicate glasseshave increased (^) our knowledge of the structures that dominate these melts, but do not give much infor- mation about reaction mechanisms.With the molecular dynamics (MD) method, it has becomepossibleto follow the motion of individual ions in a theoretical system (^) and to determine both melt structure and reaction mecha- nisms in the molten state over a wide range of pressure and temperature. An underlying concept behind olrr approach to these problems is that interatomic potbntials in condensed phases are dominated by the same short-range forces present (^) in isolated molecules (Nowton and Gibbs, 1980; Newton et al., 1980;Gibbs, 1982).Thus, it is possibleto translate data from quantum mechanical (^) calculations on molecular clusters,such as H4SiO4and HuSirOr,to inter- action potentials between Si and O in silicate minerals and melts. Using this approach,it will be possibleto put the understandingof melt dynamics into a firm theoret-
- Present address:Geophysical Laboratory, Camegie Institu- tion of Washington,2801 Upton Street,N.W., Washington,D.C.
20008,u.s.A.
0003404x/88/091H941$02.00 94r
ical framework that is simple and flexible enough to en- compassmany systemsof geochemicalinterest. In this paper,two potential-energymodelsderived from quantum mechanical calculations are applied to SiO, melts and glassesand used to extract thermodynamic, structur- al, and diffusion data. Shortcomings of using the ionic approximation to simulate covalent systemsare exposed, but the utility of the MD techniquein kinetic studiesalso becomes obvious, becausediffusion mechanisms for Si and O can be displayed using computer graphics.Prelim- inary simulations with a covalent potential are also com- pared to the more extensive ionic results in an effort to increasethe accuracyof model calculations.Inclusion of these short-range and many-body forces should make simulations more realistic for silica-rich melts.
Mor.ncuun DYNAMTcs
The MD calculation technique usesa central box con- taining at least a few hundred ions to simulate a melt or glass.A periodic boundary condition repeatsthis cell on all sidesto createan infinite volume and remove surface. efects. Within the central box, eachion interacts via short- range forces with all other ions up to a radius of 5 to 6 A (except for the longer-rangeCoulomb force). In this study, the central box was subdivided into many smaller boxes (usually 64), a'ndions were assignedto one ofthe
942 KUBICKI AND LASAGA:MOLECULARDYNAMICS^ SIMULATIONSOF SiO,
small boxes. This technique reducesby several-fold the time involved in carrying summations over both short- range and Coulomb forces. The force on each ion is ob- tained by differentiating^ the computed^ interatomic poten- tial energyfor every ion:
F,: (^) > -(dLr/dr,)(dru/du,)^ (l)
where d :^ force on ion i in the u direction (u:^ x, y, z), @u:^ potential energybetween ions i andJ, ru :^ the dis- tance betweenions i and7, and du :^ the x, y, or z com- ponent of ion i. The sum in Equation I is over all T ions including those in the periodic boxes. For the ionic model, Coulomb potential energies(A) are calculated at eachtime step by application ofthe Ewald technique to achieve efrcient convergence(Sangsterand Dixon. 1976):
tained from the force on the ion by assuming^ that ions obey Newtonian mechanics(F :^ rua). The value of Al in thesesimulationswas 1.0 x^ l0-''^ s at 4000 and 6000 K and 2.5 x^ 10-15s at 2000 K. All simulations were run on a cubic cell containing 324 ions (108 SiO, units). Runs were completed with initial temperaturesranging from 300 to 9000 K.
PorpNrru,s usED rN THE MD sTMULATIoNS
The most important factor in MD models is the ac- curacy of the potential equation that is chosento simulate the forces^ betweenions. One of the aims of this paper is to systematizethe development of potential functions for usein MD simulations.Recently,Iasagaand Gibbs (1987) reported ab initio calculations ofthe Si-O potential and applied their results to the prediction of structural and elastic properties for the silica polymorphs d.-qvartz, a-cristobalite, and stishovite (the high-pressure poly- morph). We propose^ to use these potentials, which have a firm theoretical basis and which have been tested on well-known mineral structures,^ in the MD^ simulations. In this fashion, we hope to remove some^ uncertainty from the interpretation of the MD results. Most earlier work on silicates has used simple ionic potentials of the type
v,,:Z'Zi4' t J (^) r.. * A,,exp(-r,,/p,,), (4)
. U
where (^) 2,, A,,, and (^) t4 could be varied to improve the com- parison with experimental properties.The ionic potential usedin this paper^ is a potential ofthe type given in Equa- tion 4 but derived in Lasagaand Gibbs (1987) using ab initio results.The set of constantsusedfor the ionic mod- el (Table l) are based on the usual formal chargeZ"': +4.0, used in previous MD papers(e.g.,Soules,1979; Angell et al., 1982). In this study, modeling of glassesand melts of the SiO, system was basednot only on a simple ionic potential as in Equation 4 but also on a central-force-field covalent potential. The covalent parameters^ (Table 2) were de- rived from a Morse potential fit to ab initio^ quantum mechanical potential-energysurfaces(Lasagaand Gibbs, 1987).Basedon theseapproximations for the Si-O inter- action in the HoSiOoand HuSi'O, molecules,calculations of the a-quartz structure using the wMIN program (Bus- ing, l98l) gave good results for both the ionic and co- valent potential models (Lasagaand Gibbs, 1987).How- ever, only the covalent potential yielded a good equation of state and elastic constants. The Lasagaand Gibbs data and methodology^ are used heret however,the form ofthe covalent^ potential is slight- ly modified for use in the MD^ simulations. The initial covalent potential chosen^ for the MD^ work was a com- bination of the Urey-Bradley force field and the flexible Morse potential:
-:;i [rurco"r,nt
.;>E?,s(o[sP@'
]
4=e)'
where
s(c): f, Z,exp(-iGr,,)
and the index i runs over all N ions in the cell, G is a reciprocal spacevector, S(G) is the structure factor, 4 is an arbitrary constant^ that can be chosento speedup con- vergence(0.tS A-'l was used here), Z,is the ionic charge of speciesi and V.r*^ is the volume of the periodic box. The sum over reciprocal vectors in Equation 2 was cat- ried out using the transformation introduced by Sangster and Dixon (1976),^ which decreasesthe time involved for eachcalculation by a factor of 14 the number of ions. The potential in Equation 2 is also differentiated as in Equa- tion I to add the Coulomb contribution to the total force on each ion. Once the instantaneousforce on each ion has been de- termined, its position and velocity are calculated for the next time step by using the algorithm of Schofield(1973):
r(t+ at):'f',i^lfli'- o(t- at)t^t, (3a) v ( t + A t ) : v ( t )
- t/ul2a(t^ +^ At)
- 5a(t) -^ a(t -^ At)l\t, (3b)
where r(t), v(l), and a(t) are the position, velocity, and acceleration of an ion at tr-te l. Accelerations are ob-
944 KUBICKI^ AND LASAGA: MOLECULAR^ DYNAMICS^ SIMULATIONS^ OF SiO'
F B - I N I T I OE N E R G YS U R F F C E
I----J---t.oo
t.
I.
" 1 a n (^) 13 0. r 4 0. 15 0. 16 0. 1 7 0. 1 8 0.
SiiOSiANGLE
Fig. l. (a) Quantum mechanicalpotential-energysurfacefor HuSirOr.The x axis is the SiOSi bond angle^ in degrees,and the,, axis is the SiO bond length in A. Eneryies are in kcaVmol. (b) Potential-energy surface^ for H5Si2O7calculated from Morse potential equation. (c) Energy differencesbetweenFigs. la and lb test the ability ofthe Morse potential to reproducequantum mechanical results.
7 0
o 2 o
t o l.
o
ID
where
1... \ : I^ > l/,,"o#.' The temperature can also be set according to Equation 8 by scalingthe velocities. Thus, if the desired value of the temperature is ?"oand the computed value according to Equation 8 is Z, the velocities are scaledby vi"* :^ v?'d{To/T)^. (9)
Scalingis the MD equivalent of immersing the systemin an isothermal heat bath. In this study, the ions are initially positioned in a quartz- like structure that has been distorted to fit a cubic cell and then given a Boltzmann distribution of velocities^ for a chosen temperature. Since runs simulate real times of I to 2 ps (picoseconds),the initial temperature used to randomize the system is set at 10000 K (Angell et al.,
- so that equilibrium is rapidly attained. Cooling is then done on successiveruns in 1000-K increments. The temperaturesfor each run were set by scalingpar- ticle velocities for 250 to 500 time steps.Another 1000 to 2000 steps were then calculated without temperature scaling. Cooling to the next temperature was then begun by scaling the velocities of the last time step from the previous simulation. Equilibration at each temperature was complete during scaling after 250 to 500 time steps as evidencedby the establishmentof a well-defined ther- modynamic averagefrom the nonscaledtemperaturecal- culations. Numerical errors were checkedby the conser- vation oftotal energywithin the system (Figs. 4 and 5).
Both the potential energy (Eqs. 4 and 6) and kinetic energy (Eq. 7) will fluctuate around their averagevalues, which should represent^ the bulk thermodynamic value for a macroscopic system.^ Figure 4 shows^ these^ fluctua- tions in kinetic energy,as^ representedby temperature,to be within 50/oof the average.According to statistical^ me- chanics theory (Hill,^ 1962), the standard deviation in temperature, LT/T^,', for a system with N particles is given by
AT/T: o x^ l/(3N)h. (10)
If N :^ 324, as in our system, Equation 10 predicts that the standard deviation, o, should be approximalely 60/o, consistentwith the results^ in Figure 4. In these MD^ simulations, the total energy (,8) of the systemis constant after scalinghas^ been^ removed.^ There- fore, fluctuations in the potential energy,iD,^ and the ki- netic energy,KE, must be coupled by the relation E -- @
- KE. Figure 5 is a plot of the total energyas a function of time. Note that ,E varies by 0.10/odue to numerical error, truncation of the power seriesin the Schofield al- gorithm, and approximation of the Ewald summation. We have assumedthat this small error does not signifi- cantly affect our results. Prediction of pressuresin the system is problematic becausethe instantaneous^ pressurefluctuations within the system can be dramatic (e.g.,instantaneouspressuresmay rangefrom -100^ to +100 kbar). Nonetheless,pressure estimatesare^ made at eachtime stepusing^ the virial theo- rem (Woodcock,1975)
KUBICKI AND LASAGA:MOLECULARDYNAMICS SIMULATIONSOF SiO,
SiOSi ANGLE
DI FFERENCEII1RP
t. f ,
c
" ' l > i 1 3 0. 1 4 0 " 1 5 0. 1 6 0. 1 7 0. 1 8 0.
SiOSiANGLE
Fig. l-Continued.
6 0
o
z
o
(^6) l. o (t
P : - @ A / d n r
: Nkr/v+ 1nn(j (^) 0,."r), (n) \ i : r I
where F, is the force acting on ion l, r, is the position vector, @A/dnr is the derivative of the Helmholtz free energy of the system with volume at constant tempera- ture, and N and V are the number of ions and volume of the central box, respectively. Values ofthe pressureare obtained by averagingover
many time steps(N",.o=^ 100). Figure 6 gives^ results of a typical run. Fluctuations of +10 kbar occur, but Figure 7 showsthat a time-averaged pressureconvergesto a well- defined value. The simulations using the ionic potential, however, yield pressuresthat are substantially different from those predicted from the equation of state of SiO, glassat this temperatureand density. This deviation (Ta- ble 4) stemsfrom the differencebetweenexperiment and theory in the unit cell predicted with this potential in calculations on quartz (Lasaga and Gibbs, 1987). The simulations using the covalent potential predict pressures
CFLCULFTEDENERGYSURFRCE
-----':oo
---l:-:=--
,-,-_--
/-
t---------- o-.-=.----'-
\ \
\ \
--
---;=---o.o (^) o'o]
KUBICKI AND LASAGA: MOLECULAR^ DYNAMICS^ SIMULATIONS^ OF SiO'^941
!I
t l.L u F UI oo
UE ao u^6 Eo-
x oc U o o u l E c
S I L I C A < T > = 2 8 3 K < P > = 2 k b^ S I L I C A < T > = 2 8 3 K < P > = Z k b
- o
- 0 aoo.o
Fig. 4. Variation in averagetemperatureover 100 time-step intervals for a low-temperature, low-pressuresilica glass^ using the ionic model. All scaledtime stepsof equilibration have been excluded.
p: -Q/n@v/aP)r.^ (14)
Becausethesesimulations usea microcanonicalensemble (constant N, V, and.E), several runs with various vol- umes are used, and (dV/An is calculatedfrom the change in calculatedaveragepressurewith volume. Table 6 con- tains estimatesfor d and B of a SiO, melt from our sim- ulations using the ionic model and experimentalnumbers for SiO, glass. Prediction of second-orderthermodynamic quantities is a stringent test of MD models, so the error in these predictions is encouragingly small. However, the more realistic covalent potentials are necessaryfor these^ sim- ulations if accurate^ thermodvnamic data are to be ob- tained.
Srnucrunn Glass and melt structures can be analyzed^ using the pair correlation function
S I L I C A T = 2 8 3 K < P > = 2 k b
T O T A LE N E R G Y
P O T E N T I A L E N E R G Y
.aOO. (^) O 5OO. O T I M E S T E P S
Fig. 6. Extreme calculatedpressurevariations are presenteven when averagedover 100 time stepsfor the simulation in Figs.^4 and 5. Temperaturefluctuationsdo not control the pressurevari- ations in this casebecausethe NkT/Vterm of the virial theorern (Eq. l1) is very small at 300 K. The pressurepredicted at this temperatureand density from the SiOt glassequation of stateis 0 kbar, and the average calculated pressure^ in this simulation is 2 + 5 kbar.
g,i(r) = (l / 4r pr'z)ldN(r)/drl (l 5)
where N(r) is the number of ions of typeT within a sphere of radius ," around a selectedion of type I, and p is the bulk density of the ions of type 7. Note that as r '^ @, short-range order and structure are lost, and N(r)^ ' (4r/3)(pr3). Therefore, g, '^ I as / + oo. To calculategu for these simulations, each ion is selectedover 100 or more time steps,and the number of atoms of type 7 in thin spherical shells with radii of r and r^ +^ Lr^ (Lr^ : 0.025 A) around these ions li.e., dN(r)/drl^ are counted out to l0 A. ttre correlation functions are averagedover the number of time stepsand plotted to show the short- range melt structure. Figure 9 gives the results for the ionic model at the experimental density and 300 K as well as the curve determined from Mozzi and Warren's (1969) experimental data. Our fit to the experimentally derived correlation function is very good. The peaksdis- agreeby 0.1 A or less^ for primary coordination and by
TnaLe4. Cell dimensions,densities,and averagecalculated pressuresfor MD runsat 6000K
oE I
C)E Uz U
- o T I M E (^) S T E P S
2so.0 (^) 5OO.O T I M E S T E P S
Cell length (r) (A)
Density (pl (g/cm")
Averagecalculated pressure(A (kbar)
Fig. 5. Fluctuations oftotal and potential eneryiesofthe sys- tem over 1000 time stepsafter an initial equilibration period of 500 time stepsusing velocity scaling.
O r o o 0. 0
6 1 7
40 41 4 1 45 57 60 1 0 0 175
SILICAT= 283K P= 2kb
948 KUBICKI AND LASAGA: MOLECULAR DYNAMICS SIMULATIONS OF SiO'
tt L o -os,
uE o^ f, ulo trIL c'
F o F
o L .o^ o -r
L =
ao UE IL
r o o. (^) 200. 2oo'o S T E P S ,, "?o'
orrar, 600'o 8oo'o
Fig.7. Overallrunningaverageofthe pressureplottedin Fig.
- Convergenceto a well-definedvalueindicatesthat the calcu- lation is stableand that more realisticmodel potentialsmay allow accuratepredictionsof equationsof statefor glassesand melts.
lessthan 0.5 A for secondarycoordination. This discrep- ancy is especiallynoticeable in the Si-Si correlation, and it is causedby the disagreementbetweenthe averageSiOSi angle determined experimentally (144")^ and the MD-cal- culated averageSiOSi angle of these simulations (160'). Coordination numbers (CN) are obtained by integrat- ing to the first minimum, r,, in the g(r) curve.
cN:4zrpj"" ^rrn nr. ( l 6 )
Table 7 comparesthe values of first and secondcoordi- nation numbers with experimental values (^) found by Moz- zi and Warren (1969) using X-ray diffraction. For pri- mary coordination, the values in Table 7 are in good agreement with one another. The reason for the large errors in the secondarycoordination numbers is unclear, but the breadths ofthe peaksin the correlation functions makes defining secondarycoordination numbers ambig- uous. Coordination numbers may also be obtained by the location of plateaus^ in a plot of N(r) vs. r (Fig. l0); the plateaus define the average number of ions clustered around each type ofion at a given distance.Primary co- ordination plateaus are very distinct, but there is diffi-
TABLE5. Comparisonof MD-derivedCuand experimentalCp for SiO.at hightemperatures
. From Robie et al. (1978) except the 2000 K value, which has been extrapolatedfrom their data.
T I M E Fig. 8. Plot of runningaveragefor pressurecomponentsin the x, y, and z directionsshowsconvergencetoward the same valuefor all threecomponentsindicatingisostaticequilibrium. Note that total averagepressureis slightlyhigherthan compo- nentsbecausethe (NkT/n^ term of Eq. I 1 is significantat the temperatureof this simulation(i.e.,1500K).
culty in defining secondarycoordination becauseof the diffuse nature ofthe secondarycoordination spherein a melt or glass. The average,primary polyhedral units are accurately reproduced by our models (Tables 7 and 8). However, the ionic approximation results in errors in the distribu- tion of SiO bond lengths and OSiO angles(Fig. I l). The centers of the peaks^ for SiO and OSiO distributions are in the proper positions, but the widths ofthese peaksare broader than the widths of the experimentally derived peaks.For example, our SiO bonds vary by t0.2 A and OSiO tetrahedral anglesby + 10",^ whereasX-ray studies show a +0. I A variance in SiO bonds (Mozzi and War- ren, 1969) and esn (Gaskell and Tarrant, 1980) experi- ments indicale a +7" distribution around the averagetet- rahedral angle.Wide distributions are causedby the lack of directionality in the ionic potentials. Also, high-tem- perature structuresmay have been locked in by the very rapid cooling rates of MD simulations. The covalent model corrects this wide distribution for the SiO bond lengths, but the angle distribution remains unchanged(Fig. l lb). The discrepanciesfor the value of the SiOSi anglehave been a matter of controversy for some timp. The original work of Mozzi and Warren (1969) was "corrected" by DaSilva et al. (1974) to give a value of 153.0'. Subse- quently, however, Coombs et al. (1985) have vindicated the original work of Mozzi and Warren and re-calculated
TaeLe6. lsothermalcompressibilities(B)andthermal-expansion coefficients(a)for SiO,
B (ionic MD)^ B (experimental) (bar-1) (bar')
o (MD ionic) @(experimental) ( K ' ) ( K ' )
Temperature (K)
Density (g/cm")
Cv (MD ionic) Cp(experimental)- lJ(mol K)l^ U/(mol K)l 300 1000 2000
< T > = 1 5 0 0 K < P > = 2 8 k b
T A L A V G P R E S S U R E
- 9 x 1 0 6 2. 4 5 x^ 1 0^6 3. 4 x 1 0 - n 0. 5 x 1 0 - ,
C A L C U L A T E D
C O V A L E N T
I O N I C
KUBICKI AND LASAGA:
< T > = 3 O O K < P > = 2 k b
- 2 1. 4 1. 6 l. a s i o B o N Do t s t n r u c e t i >
O - S i - O D I S T R I B U T I O N
1 0 0. 0 1 2 0 .o A N G L E ( d e g )
Notei The deleeuw et al. (1985)values are from a @ntinuous ranclom network model(CRN).The 5T06-31 G*^ valuesare^ from highlevel^ quantum mechanicalcalculationson isolated rings of three silica tetrahedra with nonbridgingoxygens saturated with H atoms' Analytic values are from optimizedring structuresin isolatedmolecularclustersbasedon equations fit to OM potential-energysurfaces.
model of Bell and Dean (1971) with the important ex- ception of the existenceof 3-membered Si rings in the MD models. These planar, highly strained structures^ have been used to explain the D, line in the Raman spectrum of SiO, glass(Galeener,1982;McMillan et al', 1984).If this assignmentis correct, the MD simulations may be a more accurate reproduction^ of glass structure than the ball-and-stick models. The ionic, high-pressurerun (Ta- ble 10) showsthat aI V/Vo:^ 0.6, 3- and 4-memberedSi rings begin to dominate the structure of the glass.This finding supportsthe assignmentofDr (Sharmaet al., l98l) and D, to small rings on the basis of the intensity in- creasesof these two peaks^ with pressure(Hemley et al., 1986).The changefrom 6-membered^ Si rings to 4-mem- bered Si rings may also have some important effectson melt diffusion, as will be discussedlater. Given the above errors in the reproduction of the SiO, system, we have studied the changesof structure with temperature and pressure that are easily obtainable with MD. Compressibilities of individual bond types and an- gles are found from runs at different densities^ so that a "polyhedral" approachto studying the melt or glassmay be adopted(Hazen^ and Finger, 1982).Knowledgeof these factors is of great importance to the study of diffirsion, nucleation, and crystal growth under geologicconditions where higher temperaturesand pressurespredominate. Figure l2 showsthe relative changesin three^ important structural featureswith pressurefor the ionic model SiO bond length, OSiO angle, and SiOSi angle. Clearly, the
lengthsfrom X-ray data (Mozzi and Warren, 1969)are^ narrower than the calculateddistributions, but the covalent model signif- icantly corrects^ the bondJength distribution as compared with the ionic model.
MOLECULAR DYNAMICS^ SIMULATIONS^ OF SiO,
TneLe9. Bond lengthsand anglesof rings^ in vitreoussilica SiO bond OS|O angle^ S|OSiangle Bulk average 160.6(ionicMD) 159.3(covalentMD) 152.7(deLeeuwet al., 1985) 144.0 (Mozzi and Wanen, 1969) 152.0(DaSilvaet al., 1974) Rings of three silica tetrahedra 108.0 140.4(ionicMD) 89.8 129.5(covalentMD)
6 3 1
6 1 3
6 1 7
Rings of four silica tetrahedra
132.0(deLeeuwet al.,1985) 133.s(5106-31G) 136.5(analytic)
159.1(ionicMD) 154.4 (covalentMD) 146.4(deleeuw et al., 1985) 152.9(analytic)
1 4 0. O
8 0. 0 (^) l O O. O 1 2 0. 0 1 4 0 .o A N G L E ( d e g )
1 6 0 .o 1 8 0. 0
Fig. I l. Plots of distribution of SiO bond lengths,OSiO tet- rahedral angles,and SiOSi intertetrahedralangleswith the ionic model results (solid lines) compared to the covalent model re- sults (dashed lines). Actual distributions of tetrahedral angles from Esnexperiments(Gaskelland Tarrant, 1980)and SiO bond
C O V A L E N
I O N I C
S i _ O - S i D I S T R I B U T I O N
C O V A L E N T
KUBICKI AND LASAGA: MOLECULAR DYNAMICS SIMULATIONS OF SiO, 9 5 1
TaeLe10. Vitreoussilicaring statistics
Atomsin pnmarynng Ringsize (%l (^) Model Reference
/Vote.'MD:^ moleculardynamics; CRN :^ continuous random network; ILP :- p : 3. 6 1 7 9 / c m 3 -^ ionic, low-pressure; IHP :^ ionic, high-pressure.
dominant effect ofincreasing pressureis the contraction of the SiOSi angle. As found in static calculations for quartz (Lasagaand Gibbs, 1987), most of the compress- ibility of Si{}, is dependent upon this parameter. (^) Tetra- hedral anglesand bond distancesare much lesscompress- ible and do not undergo drastic changeup to 100 kbar (Hazen and Finger, 1982), where Si changesto 6-fold coordination in stishovite. The tatSito r6tSitransition was observedin these (^) simulations along with a consistentin-
- (^) ao. 90. too. PRESSURE (kbor.) Fig. 12. Percentchangesin averagebond lengthsand angles with pressureshowing the relative flexibility ofthe SiOSi angle over SiO bonds and OSiO anglesin the ionic model.
U) I zo!. N l.
o.
-11.s ;j
t;" -ro.s _1o.o
Fig. 13. This strongcorrelationbetweenln D" (at 6000 K and 40 to 100 kbar) and 5-coordinatedSi suggeststhat the anomalousincreasein diffirsionrate with pressuremay be a defect-controlledprocess.The linearcorrelationcoefficienthere is 0.99.
creasein the percentageof tstSi(Fig. l3). Contraction of the SiOSi angle should be expected in this simulation, becausesimilar effectsare seenexperimentally in quartz under pressure.Levien et al. (1980) have shown that the SiOSi angle decreasesby 6.50/oand the SiO bond de- creasesby 0.260/owhen quartz is subjectedto a pressure of 60 kbar. This is a changeof -0.16"/kbar and -7^ x l0-s A"/kbar (^) in quartz; our simulations give -0.09'lkbar and -18^ x^ l0-s futbar^ over the same pressurerange. Considering that the latter values are for a melt at 6000 K and the former are for crystalline quarlz, this agree- ment is startling and may be an indication of the overall dominance of bulk properties such as compressibility by individual coordination polyhedra.
DrrrusroN MD simulations allow calculationof self-diffusioncoef- ficients by the Einstein equation
c. o
- 4
- 0 1 9. 0
s0.
28.O
6 1. 0
i o
- 9
- 9
M D C R N I L P Covalent I H P - M D C R N I L P Covalent I H P - M D C R N I L P Covalent I H P - M D C R N I L P Covalent I H P - M D C R N I L P Covalent I H P - I L P Covalent I H P -
Soules(1979) Bell and Dean(1971) this study this study this study Soules (1979) Bell and Dean(1971) this study this study this study Soules(1979) Bell and Dean(1971) this study this study this study Soules(1979) Bell and Dean(1971) this study this study this study Soules(1979) Bell and Dean(1971) this study this study this study this study this study this study
where z :^ time interval during diffirsion, lr,(r) -^ r,(0)1'z is the squareddisplacementof an ion of type i during the time interval 7, and N is the number of silicons or oxy- gensin the cubic cell. "(.. .)" in this caserepresentsav- eraging over many different initial times. Equation l predicts a straight line for a plot of (r2) vs. time, and a slope that is related to the diffusion coemcientD,. Given the time scale of MD^ simulations I ps (picosecond): 1000time steps:75 h of CPU time on a vAx tr/lsof, diffusion coemcients are usually determined at geologi- cally unreasonabletemperatures(e.g.,6000 K). However, extrapolation of high-temperaturevalues to more realis- tic conditions may be possible if the structure does not change drastically between the high- and low-tempera- ture ends and there are no changesin diffusion mecha- nism acrossthe extrapolated temperaturerange.
":#(i
r,,r")-,(o)r,)
r x * z =^ O. 9 2 5
< T > = 6 0 0 0 K
\ " !\ qs'o t \ - \
\ \
,roa,\ \ -^
--. .-
S I L I C O N # 9 1 < T > = 6 0 0 0 K < P > = 6 0 k b
aoo.o looo.o I I M E S T E P S Fig. 17. Displacementsquaredof Sio*ion after scalinghas beenremovedfor a simulationat T :6000 K andP :^60 kbar. The rapid risefrom I :^ 550to / :^600 marksa diftrsionevent duringwhich time the atom is also3-coordinated(seeFigs. and l9).
points, the activation eneryiesare 102.3kJ/mol for Si and 9l.l kJ/mol for O. These values are low compared with experimental valuesof 20G400 kJ/mol. Both atom types are governed by similar energeticswith the O moving more rapidly and diffusing more readily becauseof its lower massand becauseit is not as tightly enclosedwithin its primary coordination sphereas Si.
Activation volumes are obtained from the slope of the ln D vs. pressureline (Fig. 16). Becausethe rate ofdif- fusion increaseswith pressureup to 100 kbar, theseAZ" valuesare negative,in agreementwith experimentalfind- ings for other polymerized silicate melts (Kushiro, 1983). The activation volumes from the MD^ results are -4. cm3/mol for Si and -4.53^ cm3/mol^ for O. The changein the diffusion rate for O with pressure is greater than for Si, which hints that the transition-state complex for dif- fusion is more compact for O than for Si, and both have denser structures than normal melt polyhedra. Thus, al- though temperatureand pressurehave opposite effectson the type of defectsformed, the nature of SiO, melts causes both to increasethe diffusion rate. Above 100 kbar, the diffusion coefficients decrease.because Si has become 6-coordinated and is more tightly bound within its pri- mary coordination shell. One of the biggestadvantagesthat the MD technique enjoys is a knowledge^ of the positions of all individual ions as a function of time. Thus, reaction mechanisms may be checkeddirectly with MD, and the ionic behavior should be as realistic as the potential used to model it. Figure l7 follows the displacementof a Si ion throughout a segmentof a simulation. At time step 578, there is a large (> I A') displacement of this Si ion, and it begins to vibrate around an averageposition of 1.3^ A'compared to the averageof 0.8 A' before this step. This diffusion event coincideswith step 578 in Figure 18, which shows that the Si ion is bonded to three oxygen atoms at this point. The structure of this complex is nearly trigonal
KUBICKI AND LASAGA: MOLECULAR (^) DYNAMICS SIMULATIONS OF SiO,
N
N
L
- o (^) 600. o I 200. o
trl(J z
a o a z o gl
o
- O (^) 7 5 0. O T I M E S T E P S
8 5 0. O 9 5 0. O
Fig. 18. Distancesofnearest-neighbor O atoms to the Si atom ofFig. 17 plotted vs. time. The averageSiO bond is l.-. (^) A in this case,and the maximum distancefor a bond is set at 2.0 A in thesesimulations. Note that the loss of no. I coincideswith the retumofno.2toanormalbridgrngOposition.Also,no.5wasabridgingOonanearbytetrahedronthatformeda3-membered Si ring with the central Si atom before the loss ofno. I occurred.
< T > = 6 0 0 0 K < P > = 6 0 k b
954 KUBICKI AND LASAGA: MOLECULAR DYNAMICS SIMULATIONS^ OF SiO,
Fig. 19. Moleculargraphics(onrrr) of step578 in Figs. and 18.Theformationof this trigonalplanarcomplexcoincides with a changein the averagepositionofthe centralSi atom in the targetedtetrahedra.Bond distancesare givenin A.
planar (Fig. l9); the,r coordinatesofthe Si and all three bonded oxygensare within 10.5 A ofeach other, and the OSiO anglesare approximately 120".Within this frame- work, it becomespossibleto visualize^ Si diffusion normal to the equatorial oxygenstoward the upper right ofFigure 1 9. At time step 700 in Figure 18, an O transfer begins betweentetrahedra.As the new O becomesbonded to the central Si, a 5-fold complex is formed that is a distorted square pyramid (Fig. 20) corresponding to step 775 in Figure 18. A nearly planar ring with three silica tetrahe- dra is also formed at this time using the 5-fold Si. This highly strained complex brings another Si within 3 A of
Fig. 20. Distorted squarepyramidal complex from step 780 of Fig. 18 showing a mechanism of O transfer from one tetra- hedron to another. No. I is the O leaving the original tetrahe- dron, and no. 5 is the new atom replacing it. All five of the oxygens bonded to the central Si are bridging oxygens, but two of the surrounding tetrahedra have been eliminated from^ the drawing for clarity.
x -^ Z N B r O
- l i r. 0 (^) - 1 3. o - l?.! - t l. o - l o. o - 9. o - a. o
Fig.2l. Percentageof 3-coordinatedSi andnonbridgingoxy- gensareplottedvs. ln D (at V :^ Voand2000to 6000K) here to illustratethe breakdownof polymerizationand loweringof viscositywith increasingtemperature.Correlationsbetweenthese defectsandln D for Si andO areboth^ 0.99,whichis suggestive of a defect-controlledmechanismfor difusion.
the central Si atom and completes the oxygen transfer (Fig. l8). Within another l0 time steps,one of the orig- inally bonded O atoms leavesthe squarepyramidal com- plex, and the central Si returns to a tetrahedral geometry without the 3-membered Si ring. This mechanism involves a short-lived 5-coordinated Si that the pressure^ derivatives of the diffusion coem- cients indicate may be a defect^ responsiblefor the anom- alous pressure^ dependenceof diffusion. Other workers (Brawer, I 98 I ) have seenthe effectof 5-fold coordination on diffusion in similar systems(BeFr),^ and our data sup- port this conclusion, becausethere is a strong correlation betweenthe percentageof 5-fold Si and ln D in our data (Fig. l3). In addition, the increaseof ln D with temper- ature is correlatedwith 3-coordinatedSi and nonbridging oxygenspecies(Fig. 2l). The role of 3- and 4-membered Si rings in difusion at high pressuresis not clear but well worth investigating. If these strained speciesdo play an active role in melt diffusion, they could explain the anomalous pressurede- pendenceof diffirsion coefflcientsfor network-forming ions while at the sametime slowing the diffusion rates^ of net- work-modifying cations. Computer animation of these simulations may be necessaryto study many diffirsion eventsand determine which structuresare^ responsiblefor negative activation volumes in silicic melts and which are incidental. At higher temperatures (i.e., T^ >^ 6000 K), different mechanisms apparently become the controlling factors for diffirsion becausethe ln D vs. (l/7") plot (Fig. 15) has a sharp increasein slope at this point and 3-fold Si and nonbridging O become abundant. Above 100 kbar, the 3-fold and 5-fold mechanisms^ are no longer applicable either, becausethere is a large decreasein the diffusion coefficientat this point and the transition to a melt dom- inated by 6-coordinated Si has taken place.
o 2 ,".
E
F .o.
(t) o I 20.
z
o
!Q ro.
