Application of empirical ionic models to SiO2 liquid: Potential model approximations and integration of SiO2 polymorph data

Cosmoch~mca AclaVol.51,pp.

Grochrmmz d C Pergamon Journals

Ltd. 1987. Printed

1209-1218 in U.S.A.

00 16-7037/87/$3.00

+ .oO

Application of empirical ionic models to SiOz liquid: Potential model approximations and integration of SiOz polymorph data ROBERT L. ERIKSON and CHARLES J. HOSTETLER Pacific Northwest Laboratories, Richland, WA 99352, U.S.A. (Received January 2 I, 1986; accepted in revised form February 13, 1987) Abstract-Structural and thermodynamic properties of crystalline SiOl and SiOz liquid have heen examined with Monte Carlo (MC), molecular dynamics (MD), and energy minimization (EM) calculations using several ionic potential models obtained from the literature. The MC and MD methods calculate the same structural and thermodynamic properties for liquids when the same potential model is used. The EWALD (192 1) method of calculating coulomb interactions reproduced most successfully the structure of liquid silica. Approximating the coulomb interaction by eliminating the inverse lattice sum results in predicted bond distances that are too short and an average (Si-0-a) angle of approximately 180”. Introduction of a cut-off in the potential energy function produces irregular tetrahedra and inconsistencies in predicted Si-0 coordination in silica liquid. The system internal energies show that liquid structures derived from random starting configurations can be metastable relative to structures calculated from crystalline starting configurations. The static lattice properties of the polymorphs alpha-quartz, coesite, and stishovite were used to evaluate further the accuracy of different sets of repulsive parameters for the full Ewald ionic model. Most of the models studied reproduced poorly the measured structures and elastic constants of the polymorphs. The major weakness of the ionic model is the unreasonably large Si-0 bond strength (120 X lo-‘* ergs/bond) when formal ionic charges are used. Fractional charge models with a small Si-0 bond strength (30 X IO-‘* ergs/bond) improve the agreement with experimental data. However, further improvement of the ionic model should include reducing the Si-0 bond strength to values in better agreement with published estimates (7 x IO-‘* to 13 X lo-” ergs/bond). By using additional information to constrain the parameterization of the ionic model, such as estimated bond strengths and static properties of the silica polymorphs, a model more representative of the interparticle interactions may be obtained. INTRODLJCIION SIMULATION METHODS using potential energy functions have been applied to study the properties of liquid, glassy and crystalline silicates. The computational methods used include Monte Carlo (MC), molecular dynamics (MD), and energy minimization (EM) calculations which allow direct observation of microstructure and the correlation of thermodynamic and structural properties. In MD and MC techniques, the liquid or glass structure evolves from an initial set of particle positions. The initial configuration of particles is based either on a crystalline analogue of the liquid (glass) or from random particle positions, both of which possess a potential energy different than that of the equilibrium liquid (glass). The EM technique is used to calculate the structure of a crystal with the unit cell geometry and particle positions constrained by the symmetry relations for a particular space group. The cohesive energy is minimized by varying the space group parameters and structural, elastic, and dielectric properties can be calculated for the minimum energy crystal. These three simulation methods, although computationally distinct, ah require an interaction model to describe either the interparticle potentials (MC, EM) or forces (MD). The ionic model has been most often used as an approximation representing interpatticle interactions in crystalline and amorphous silicates. These interactions are expressed in terms of paitwise additive potential (force) functions which give the potential energy

of (or force between) pairs of ions as a function of their separation distance. Using this model, the functions represent the interaction between pairs of non-polarizable, soft-shell ions and contain a long-range coulomb term and repulsive terms arising from the overlap of electronic charge distribution at small separations. The short-range potential between particles i and j is given by: AijTBu

(1)

which is the form given by PAULING (1928), or by: Aij exp( -Bijr)

(2)

which is the Born-Huggins-Mayer form (TOSI and FUMI, 1964). In these expressions Aij and Bij are parameters adjusted to fit the interaction model to some physical properties of the system, and r is the interparticle separation. The coulombic potential has the form: ZiZjh/r

(3)

where e is the unit of electronic charge and Zi is the charge on the ith ion. For a crystalline system with periodic boundary conditions, the total coulomb energy is calculated using the EWALD (1921) method (e.g., TOSI, 1964; p. 112). The parameterization of a particular potential in EM calculations usually involves fitting values for the repulsive constants using experimentally measured structural and physical properties of the crystal. Much 1209

K. 1. Erikson and C. J. Hostetler

1210

work using this technique has been applied to crystalline magnesium silicates important in understanding mantle mineralogy. Models used to simulate structural properties of these silicates have used either formal ionic charges, or fractional charges to approximate the covalency of the Si-0 bond (PRICF and PARKER, 1984; CATLOW rt ~1.. 1982; PARKER c1 al., 1984; PRICE er ul., 1985; MATSL!Iand MATSUMOTO, 1985). Covalency of the Si-0 bond has also been approximated by adding Morse potential terms to the coulombic potential (PRICE and PARKER, 1984: PRICE et al., 1985). In general. the fractional charge rigid-ion models used in EM calculations have been the most successful in reproducing the observed structural and elastic behavior of crystalline silicates. An alternative to applying fitting procedures to empirical data is the modified electron gas approach recently used by POST and BURNHAM(1986) to calculate repulsive parameters for TiOz polymorphs, quartz, forsterite, and diopside. Our major interest is the application of interaction potential models to calculate the structure of silicate liquids. Application of either MC or MD calculations to silicate liquids and glasses using ionic potential models introduces additional difficulties in the theoretical treatment and parameterization of the potential. Several methods have been developed to approximate an extended system from the finite system actually modeled. The model system consists of a cubic box which is expanded through the use of periodic boundary conditions. The nearest or minimum image convention (BRUSHet al., 1966) is used so that interactions between particles are confined to those within the basic box, or with the nearest image of particles in one of of the previous

the adjacent boxes. Using this convention, the expression for the total coulomb energy is modified by removing the sum over direct lattice vectors. This full Ewald expression represents the coulombic potential of a periodic (i.e., Pl symmetry) ionic system where the minimum image convention minimizes the error introduced by the small size of the model system.

However, additional approximations have been introduced in the application of the ionic model to silicate liquids and glasses that make it difficult to judge the adequacy of the ionic model to calculate structures of amorphous silicates. Because the summation over inverse lattice vectors in the full Ewald expression is expensive to compute, some studies (WOODCOCKet ul., 1976; SOULES, 1979; SOULES and BLJSBEY, 1981: GAROFALINI, 1982, 1983: KELSO ~‘1ul., 1983; DOAN. 1984) formed a new potential by eliminating the inverse lattice sum. We refer to this expression as the abbreviated

Ewald approximation

as an approximation dition, contributions

that can be viewed

to a fully ionic potential. In adto the potential energy from any

particles separated by more than half the boxlength (L) have also been neglected in some studies (SOULES.

1979; SOULESand BUSBEY,1981; MATXJI et al., 1982). This approximation is equivalent to a cut-off in the interparticle

potential

functions

at 1./I.

The application of MD and MC techniques to silicate liquids and glasses also requires an adequate parameterization of the model. The empiric-1 method of deriving the repulsive parameters is less rigorous for liquids than for crystalline analogs because of the absence of a long-range periodic structure to constrain the fit of the parameters. Most previous studies involving the structures of amorphous silicate compositions have used molecular dynamics and Ionic potential models parameterized to model one particular aspect of the structure (most often, bond distances and angles). However, the structure of liquids and glasses in these studies often show significant differences in bond distances, angles, coordination numbers, and the ratio of bridging to non-bridging oxygens. These dif:. ferences may be the result of the potential, the approximations used in some calculations, or differences in application of the simulation technique. Therefore, the effects of the approximations mentioned above must be evaluated in order to determine the adequacy of the ionic model to reproduce the structural properties of silicate liquids and glasses. In this study, we investigate the ability of the Ionic model to reproduce accurately liquid and crystalline structures in the Si02 system. We chose to model thy SiOZ system for several reasons. This system has obvious geochemical significance. and accurate neutron diffraction data for vitreous silica and several silica polymorphs are available. Secondly, experimental measurements of the structural. physical, and thermodynamic properties of binary alkali silicate and aluminosilicate melts suggest a regular modification of the framework structure with alkali addition. An interaction model that could describe the changes in properties as a function of composition (without having to reparameterize the framework interactions) would be more useful than models requiring reparameterization for each composition investigated. The effects of starting configurations, temperature histoq. and the form of the pair potentials (and associated approximations) on the calculated structures of SiOz liquid are systematically examined. Calculation of lattice transition energies and elastic properties of silica polymorphs using EM methods are used (with the same parameterizations) to examine the suitability of ionic potentials to model the thermodynamic and structural properties of Si02 crystals. The combination of calculations for both crystalline and amorphous Si02 enables a more complete examination of the suitability of the ionic model to reproduce interparticle interactions in SiO:

PROCEDURE: LIQUID In the MD method, the 3N equations of motion are numerically integrated and the evolution of the system is calculated as a function of time (ALDERand WAINWRIGHT, 1959). MD requires initial particle velocities which are selected at random from the Boltzmann distribution appropriate to the run temperature. In the MC method, system configurations are randomly selected such that the ensemble distribution of system energies satisfies the Boltzmann distribution. Theo-

Ionic models of SiOZ liquid retically, if the model system is in thermal equilibrium, liquid structures calculated using the MC and MD methods will be the same. The structure of liquid silica was predicted at a temperature (T) of 2000K using both the MC and MD techniques. The number of ions (N) in the simulated system was 192 and the molar volume (V) was chosen to correspond to actual (or extrapolated) values for silica liquid at the run temperature (BOCKRIS etal., 1956; RIEBLING, 1966). The MD code used (MDIONS) was originally written by ANASTASIOU and FINCHAM(I 982) but was modified to calculate properties in silicate systems, and the MC code (MCARL) was developed by HOSTETLER(1982).The MC simulations were performed in the canonical (N, V, T) ensemble, and the jump parameter was chosen to keep the configuration acceptance ratio near 0.5. Both (N, V, T) and (N, V, E) MD runs, where E is the total energy of the system, were performed using a timestep of 1 X 1O-‘5sec. Two repulsive parameterizations of the ionic potential model were used to calculate the structure of silica liquid; the model described by SOULES (1979),and the potential model of MATSUIef a/. (1982). Both models use Born-Huggins-Mayer repulsive terms and ionic charges of +4 and -2 for Si and 0, respectively. For the SOULES ( 1979) potential, we examined the abbreviated Ewald treatment with and without the L/2 cut-off approximation as well as the full Ewald treatment of the coulomb term. The full Ewald treatment was used in all calculations using the potential model of MATSUIet a/. ( 1982). Initial configurations were derived from the structure of beta-cristobalite, the only cubic silica polymorph. Three crystal structure refinements for beta-cristobalite were used as starting configurations in the simulations. Bond distances, angles, and the space groups for the three starting configurations are listed in Table I. A random starting configuration was also generated in which no ions overlapped and local charge balance was approximately satisfied. The approach to an equilibrium system energy was evaluated for all simulations by monitoring the change in internal energy with configuration number (MC) or the number of timesteps (MD). This procedure was used for the different starting configurations derived from beta-cristobalite and random ion positions. In addition, two different temperature histories were used in the simulations. In the direct melting (DM) simulations the run temperature was fixed at 2000K. Homogenization and cooling (HC) runs were started at 6000K, and the final configuration ofthe homogenized liquid at 6000K was used as the initial configuration for a run at 2000K. Homogenization of the liquid at 6000K produces disorder arising from the thermal broadening of the nearest neighbor bond distances. The use of several different temperature profiles allows the approach to equilibrium for the liquid to be made from order to disorder (DM runs) and from disorder to order (HC runs). The structure of silica liquid computed in this study is directly compared to experimental diffraction measurements. The total pair correlation function (T(r)) was computed from the predicted individual pair correlation functions using the neutron scattering factors of BACON(1972). An experimental

T(r) function for silica glass was obtained from the neutron diffraction measurements of WRIGHTand SINCLAIR(19%). PROCEDURE:

CRYSTALLINE

POLYMORPHS

Using an energy minimization technique, the static properties for alpha-quartz, coesite, and stishovite were calculated using various parameterizations of the ionic model developed to calculate the properties of amorphous SiO*. The coulomb term was treated using the Ewald transformation with direct space sum and inverse space sum limits of 3 and 5, respectively. In these EM calculations, the lattice parameter was varied to find the minimum energy unit cell size for each potential model while keeping the symmetry of the ion positions and unit cell constant. No thermal corrections to the ion positions

Table 1.

Reflnenent

1211 The space group, bond distances, and bond angles for Scrlstobalitc used as starting configurations

space group

'si-O(X)



Wyckoff (1925)

Fd3m

1.54

180°

Barth (1932)

P213

1.59 - 1.68

180' - 130'

Pewor

Fd3m

1.615 1.607

147.70

(1973)

were made and the bulk properties were calculated using the static lattice equation of state (TOSI, 1964). The calculated

structural and elastic properties correspond to values at zero pressure and temperature. The alpha-quartz, coesite, and stishovite polymorphs were studied because of the availability of crystal structure refinement data and bulk modulus measurements. For alphaquartz and coesite, low-temperature neutron diffraction measurements were available (LAGERet al., 1982; SMYTHet al., 1985, respectively). Because no low temperature structure refinements were available for stishovite, we used the refinement of SINCLAIRand RINGWOOD(1978)

at 298K. The isothermal bulk moduli for the polymorphs were measured at 298K which are the nearest experimental analogue to the bulk moduli of the static crystals.

RESULTS

AND DISCUSSION

Measurements of the radial distribution function (RDF) or total pair correlation function (T(r)) for liquid and vitreous SiOz have been reported by several investigators (MOZZI and WARREN, 1969; WRIGHT and SINCLAIR, 1978; MISAWA et al., 1980; WASEDA, 1980). The average structure measured for Si02 liquid or glass using either X-ray or neutron diffraction techniques is similar although better resolution of the average structure is possible using neutron diffraction (WRIGHT and SINCLAIR, 1978). All previous measurements indicate the average nearest-neighbor bond distances to be in the range 1.60- 1.62A (Si-0), 2.62-2.65A (O-O), and 3.10-3.13A (Si-Si). Average intratetrahedral (0-Si-0) and intertetrahedral (Si-0-Si) angles consistent with these measurements range from 109.8”-I 10.4”, and 146.2”-156.0”, respectively. The next-nearest neighbor coordination numbers (Nij) calculated from the RDF data show Ns,.o = 4, No., = 6, and Nsi_si = 4 suggesting the Si04 tetrahedra is the basic structural unit in the liquid and glass. Because of the increased resolution and wavevector independent nature of the neutron scattering factors for Si and 0, we compare our calculations directly to the neutron measurements of WRIGHT and SINCLAIR (1978). Evaluation

of ionic potential model approximations

The major problem evaluating the adequacy of the ionic model for amorphous silicates based on previous work (Table 2) are approximations introduced to reduce computation time. Earlier investigations of the bulk structure of SiOz and Na*O-SiOz glasses (WOODCOCK et al., 1976; SOULES, 1979) introduced the abbreviated Ewald approximation (i.e., elimination of the inverse space sum in the Coulombic energy expression)

R. L. Erikson and C. J. Hostetler Table

2.

Woodcock

Sumaary

et

of

MC and

EID studies

of

silicate

liquids

al.

and

EHN Abbrev.

(1976)

glasses

HC,

b

HC,

random

HC,

random

HC.

rand”,,,

HC,

random

HC,

cryeraillnr

HC,

random

“C,

random

HC,

rando*

HC,

random

HC.

random

Ewald

Ji BHM Soules

and

Abbrev.

Busbev

(1981)

LIZ

So&

and

Evald

cutoff

Varshneya (19Rl) Angel1

et

(1981,

al.

Na

0-S102

19 3

1982)

i”“s

I iquid

,

glase

nitre

et

SiO2

al.

375,

(1981)

3000

LOW

Mitre

lfquid,

(1982)

Slash

MSO-Al20

Hostetler

PaulinS

-

19 2

sio2,

(1982)

Full

Evald

ions liquid M@-SW* 21s

BHM Full

ions

liquid

L/2

/

Evald cutoff

glass Soules

and

Bnsbey

(1983)

Ne

O-S10

Shielded

3040 i,,t

Garofalinl

BHM

(1982)

Abbrev. a

Garofalini

St02

(1983)

288

al.

et

Kel*0

BHU

(1984)

Evald

,glass

B

K o-s102

BHM

ions

Abbrev.

g1nxs Doan

Evald

Abbrev.

ions

2a 8

(19a3)

ionic

E.&d

a

sio2

BHM

h4a

ions

Abbrev.

glZ%?e

Evald

il

a”d (1984)

Hitra parker

Deslpaeyand Kavamura

aHN

(1984)

Full

Evald

n Dempsey (1984)

et

al.

BHM Full Evald R

-.__ _~.~_I__-._ ;

No

information

given

about

cutoff

No = No

inf”cPation information

given given

about about

starting treatllent

BtIpl -

BorPBuggios-Mayer

repulsive

Pauliog - Pauling repulsive HC - homogenization-cooling D”

-

direct

melting

tempereture

configuration of coulomb

term

term

term temperature sequence

and an L/2 cutoff in the potential energy functions. These approximations have subsequently been used in calculations involving bulk and surface structural differences in glasses (GAROFALINI. 1983; KELSO et ul.. 1983), evaluating vibrational spectra (GAROFALINI, 1982), and calculating point defects in silica glass (DOAN, 1984). On the other hand, these approximations have not been adopted by all investigators. To separate the effects of the repulsive parameterization from the treatment of the coulombic potential on the

sequence (see

(see

text)

text)

calculated structure, we have recalculated the structure of amorphous Si02 using the MD and MC techniques with and without these approximations. In the first set of calculations, the abbreviated Ewald approximation and the repulsive parameters of SoUL_t:s (I 979) were used. The equilibrium structure ofthe liquid at 2000K is compared to the neutron diffraction measurement of vitreous silica in Fig. 1. The calculated structure was obtained from a set of MD and MC runs using two temperature histories (direct melting (DM)

Ionic models of SiOz liquid

r i*.l

FIG. 1. The structure of silica liquid at 2OOOKcalculated using the abbreviated Ewald approximation and the repulsive parameters of SOULES(1979).Solid and dashed curves represent predicted results and the experimental neutron diffraction data (WRIGHT and SINCLAIR,1978) for the total pair correlation function (T(r), barns/A2) of vitreous silica, respectively.

and homogenization-cooling (HC)) and three crystalline starting configurations. Figures 2a and 2b compare the change in internal energy of the system starting from the three crystalline starting configurations. At the beginning of the DM simulations (i.e., configuration number (Fig. 2a) or timestep (Fig. 2b) = 0), the starting configuration from the WYCKOFF (1925) refinement has the lowest internal energy of the three crystal structure refinements. The preference for this idealized beta-cristobalite structure (Si-0-Si = 1SO”)

FIG. 2. Approach to an equilibrium system energy for SiOZ liquid at 2000K predicted in direct melting (DM) runs from three crystalline starting configurations using (a) Monte Carlo and (b) molecular dynamics.

1213

over the PEACOR (1973) structure (Si-0-a = 147.7”) suggests that this potential model will not accurately predict the bond distances and angles for this polymorph. (The inconsistencies of the earlier structure refinements were discussed by FEACOR, 1973). The changes in internal energy of the system during these simulations represent an early period of structural rearrangement during which the internal energy changes rapidly, and a longer period in which the internal energies converge to the same value (approximately - 12525 kJ/mol). These results show the equivalence of the MC and MD techniques in predicting the structure of the liquid at 2000K using the same ionic model. The final state of the liquid was also obtained by ordering the homogenized liquid using an HC temperature sequence (Figs. 3a and 3b). A comparison of the calculated and experimental T(r) functions in Fig. 1 shows that this potential poorly models the average structure of the liquid. The calculated first nearest neighbor distances for (Si-0), (O0), and (Si-Si) are approximately 1.55, 2.54, and 3.09& respectively, and are all shorter than measured values. Although the predicted (and experimental) (Si0) and (O-O) bond distances are consistent with regular tetrahedra ((0-a-0) - 1 lo”), the shorter Si-0 and Si-Si average distances predicted using this potential model are consistent with (Si-0-Si) - 17 1’. The

FIG. 3. Comparison of approach to equilibrium system energies for SiO2 liquid at 2OOOKpredicted in direct melting (DM) and homogenization-cooling (HC) runs using (a) Monte Carlo and (b) molecular dynamics.

1214

K I _. Erikson and C. J. Hostetler

average bond distances and angles of the predicted liquid structure more closely resemble a disordered form of the idealized beta-cristobalite structure rather than that of silica liquid (or glass). Furthermore, the calculated total T(r) function shows considerably more structural detail between 3.5-7.OA that is not present in the experimental measurements. The addition of a cut-off in the potential functions at L/2 was considered a reasonable approximation in previous studies because the change in potential energy caused by the cut-off is only four percent when used in conjunction with an abbreviated Ewald potential (SOULES, 1979). The structure at 2000K using this additional approximation and the repulsive parameters of SOULES (1979) has been recalculated using the DM and HC reversal method adopted in this study (Fig. 4). This structure was reproduced using both the MC and MD methods and represents the minimum energy structure for this potential model. Introducing the 1~ 2 cut-off increases the (Si-0) and (Si-Si) nearest neighbor distances to 1.61 and 3.17& respectively. These average values and the intermediate range structure obtained are in better agreement with experimental measurements than the results computed using this potential without the cut-off. However, the (O-O) nearest neighbor distance is too short (2.59A) and has a wider distribution of distances implying the existence of irregular tetrahedra in silica liquid. The structure in Fig. 4 most closely resembles the structure of silica glass calculated in previous work using this potential model (SOULES, 1979; GAROFALINL 1982. 1983), or using a slightly different repulsive parameterization (WOODCOCK et al., 1976). The prediction of irregular tetrahedra, thought possibly to be caused by lack of directional character in the ionic model (WOODCOCK et ul., 1976), is the result of the L/2 cut-off in the potential functions. Full Ewald ionic model The remaining studies listed In Table 2 have used a full Ewald (or equivalent) ionic potential model to cal-

RG. 4. The structure of silica liquid at 2000K calculated using the abbreviated Ewald and L/2 cut-off approximations, and the repulsive parameters of %HJLES (1979). Solid and dashed curves as in Fig. 1.

,,m,”

I+‘IG.5. Comparison of approach to equihbrium system energies for SiOz liquid at 2000K calculated in direci melting (DM) experiments from crystalline and random startmg con. figurations.

culate the structural properties of vitreous ,tr hquld silicates. Because of the computer intensive nature of these calculations only two repulsive parameterizations (SOULIS, 1979; MATSUI et ui. 1982) used with thi Ewald transformation were reinvestigated using the MD method. Both potentials use Born-Hug&s-Mayer repulsive parameters with formal ionic charges. ‘The SOULES ( 1979) potential, used in previous calculations in this study, was parameterized to match experimental first-nearest neighbor bond distances for vitreous silica. The second potential used is the Cl-c potential oi MATSIII et al. (1982) obtained from energy minimization calculations of the structure of forsterite and applied by MATSUI et al. ( 1982) to forsteritc. enstatite. and silica melts. The structure of silica liquid at 2000K was calculated from crystalline and random initial configurations using the repulsive parameters of SOULES (1979). The Cl-c potential (MAJSLI cv
Ionic models of SiOz liquid

FIG. 6. The structure of silica liquid at 2OOOKcalculated from a crystalline starting configuration using the full Ewald coulomb potential and the repulsive parameters of SOULES ( 1979). Solid and dashed curves as in Fig. 1.

space summation is included in the coulomb calculation. Increasing the magnitude of the O-O repulsion to obtain better agreement with the experimental O0 average bond distance (SOULESand VARSHNEYA, 198 1) is not necessary if the coulomb and cut-off approximations are not used. The differences in the two calculated structures using the Ewald transformation (Figs. 6 and 7) are the result of the repulsive parameterization. Although the repulsive parameterization does not appear to have an effect on the silicon-oxygen coordination (Table 3), small differences between the distribution of bond distances are evident. For the SOULES (1979) potential, the O-O bond distance distribution is broader than the calculation using the C lc potential and the corresponding neutron diffraction peak. The Si-0 and Si-Si peaks for the Cl-c potential model are narrower than the corresponding neutron diffraction peaks, although both models are in better agreement with experimental results without the coulomb and L/2 cut-off approximations. The computed lattice energy for the equilibrium liquid (- 12330 kJ/ mol) using the repulsive parameters of SOULES(1979) is approximately 1500 kJ/mol larger than the energy for the liquid using the C 1-c potential. This is the result of the larger O-O repulsion of this model which causes the broader distribution of O-O distances. The system energies computed in the runs started from the crystalline and random configurations using the SOULES( 1979) repulsive model both reach steady state values (Fig. 5). Several of the studies listed Table 2 have used random ion positions as an initial configuration and allow the system to relax over many time-

sou1es

(1979)

llatsui et (19821

al.

o-o

2.64 1.61 3.17

6.0 4.0 4.0

llO.1°

159.8’

si-0 si-si O-0 si-0 si-si

2.65 1.62 3.18

6.0 4.0 4.0

109.8”

157.90

1215

FIG. 7. The structure of silica liquid at 2000K calculated from a crystalline starting configuration using the full Ewald coulomb potential and the repulsive parameters (Cl-c) of MATSKJI et al. ( 1982). Solid and dashed curves as in Fig. 1. The T(r) value for the predicted first nearest neighbor Si-0 interaction is 22.0.

steps until significant ion diffusion has occurred. However, the energy of the system started from the random configuration is greater than the energy of the system started from the crystalline starting configurations by 70 kJ/mol suggesting the liquid structure derived from the random positions is at a metastable equilibrium. The mean squared displacements (msd) of the ions in these simulations are vastly different because of the difference in the number of timesteps used for each run. The msd for each ion is on the order of 0.5A for the system derived from the crystalline configurations. On the other hand, the msd are on the order of tens of angstroms for the longest run ( lO,OOO-timesteps) started from the random configuration. This contrast points out that large msd of the ions (i.e., significant ion diffusion) is not sufficient to demonstrate that the calculated structure represents the lowest energy state possible. The calculated structure for the random start after 5000 timesteps is shown in Fig. 8. The O-O and Si-Si first nearest neighbor peaks are very broad and there is essentially no intermediate range structure. Further MD runs involving an increased number of timesteps

FIG. 8. The structure of silica liquid at 2OOOKcalculated from a random starting configuration using a full Ewald potential and the repulsive parameters of SOULES(1979). Solid and dashed curves as in Fig. 1.

1216

R. L. Erikson and C. J. Hostetler

and different temperature profiles did not change the calculated structure of the randomly started liquid. It is not understood at present how some studies fr.g.. MITRA, 1982) obtained a vitreous SiOz structure from a random ion configuration that is in excellent agreement with experiment. MITRA (I 982) also noted that the model structure of SiO;! melt at 2500K was highly distorted relative to SiO2 glass at 298K. This difference has not been confirmed ex~~men~lly (WASEDA, 1980) and is not evident in the minimum energy structures shown in this study. These discrepancies may be caused by the potential model used and emphasize the importance of using internal energy reversals to identify the lowest energy structural state. integration of silica pu~~rn~r~h duta to evaluate ionic model The ordered arrangement of ions in the crystalline polymorphs in conjunction with their physical properties provides additional information to evaluate the pamme~~~tion of intemction models used previously to calculate liquid structure. An accurate potential model should be consistent with the structural and thermodynamic properties of both Si02 liquid and crystals. Energy minimization (EM) lattice calculations of several crystalline SiO2 polymorphs were used to determine the suitability of the full Ewaid ionic model to represent the interactions in crystais as well as tiquids. Our calculations differ from those of other studies that have used the energy minimization technique in that the ion positions are fixed by the symmetry of the polymorph and the cohesive energy is minimized by varying the volume of the unit cell. Of the potential models investigated, three models used a Born-MayerHutins Form of the repulsive terms (WOODCOCKCT al., 1976; SOULES,1979; MATSUI et ul., 1982) and two models used the Pauling form (HOSTETLER,1982; MITRA, 1982). All models used formal ionic charges for Si and 0 except the model of MITRA ( 1982) where the fractionaf charges 4-2.272and - 1.136 were used, respectively. The results comparing the minimum energy structures of the polymorphs to experimental bulk properties for each set of repulsive parameters are listed in Table 4. The differences between the lattice energies represent the relative stabilities of the po~ymo~hs at zero pressure and temperature. Of the parameter sets studied, only the sets of MITR4 ( 1982) and MATSUI PI al. ( 1982) predict that alpha-quartz is the stable phase. These two parameter sets predict lattice transition energies of 12 and 20 M/mole for alpha-quartz/coesite and 2 i 3 and 32 1 W/mole for alpha-qua~fstishovite, respectively. The closest experimental equivalent is the difference between the enthatpies of formation (5 kJ/ mole for alpha-quartz/coesite, 40 M/mole for alphaquartz/stishovite) at 298K (ROBIE et al., 1978). The comparison between the calculated and experimental lattice parameters show that the models predict the ao of the ~lymo~hs within several percent (i.e., 0.1 to 0.2 angstrom). The unit cell edges predicted by

A Le+.an and Pravitc (1981) e Sinflair and lUmpmod (1978) f Lieberman et al. (1976)

the WOODCOCKet al. (1976) parameters arr‘ unifor, mally too small while the SOULES( 1979) and MATSCX d al. f 1982) pr~ictions are too large. The HOS’TF~Z~K (1982) and MITRA (1982) parameter sets are the most accurate for alpha-quartz and coesite, but do not predict the stishovite unit cell size with the same accuracy. For all of the parameter sets examined, the predicted bulk moduli showed the correct variation ~i.~,.,stishovite > coesite > alpha-quartz) but were larger than the ex~~mentaily measured values. The parameter set of MITRA ( 1982) most successfully reproduced the magnitudes of the polymorph bulk moduli. Of the properties examined, the bulk modulus appears most sensitive to variations in repulsive parameterization. These results are consistent with other studies that have used the EM technique for modeling the structure of silicates. The ionic model is adequate in reproducing structural characteristics to several percent. The predicted u. values (and unit cell volumes) can be either smaller or larger than experimental values independent of whether the Pauling or Born-Hug&s-Mayer form of the repulsive terms are used. The major characteristic responsible for the poet agreement between experimental and predicted elastic and thermodynamic data is the unrealistically large Si0 bond strength in the ionic model. All of the ionic models studied that use formal ionic charges are consistent with an Si-0 bond strength of approximately 120 X lo-” ergs/bond. Evidence from the~~herni~~~ calculations, viscosity data, and quantum mcchanicai

Ionic models of SiO? liquid

calculations of SiO and SiOz molecules all suggest that the Si-0 bond strength is between 7 X lo-” and 13 X IO-l2 ergs/bond (e.g., PAULING and PAULING, 1975; INGERY et al., 1976; BREWER and ROSENBLATT, 196 1; URBAIN et al., 1982; PACANSKY and HERMANN, 1978). This explains why the fractional charge model of MITRA (1982) with an Si-0 bond strength of 30 X lO-‘2 ergs/bond was the most successful in repro-

ducing the properties of Si02. Data from studies of other silicates also suggest that fractional charge ionic models more accurately reproduce structural as well as elastic properties (PRICE and PARKER, 1984; PRICE et al., 1985). CONCLUSIONS The ionic model can reproduce many of the structural features of amorphous silicates when the Ewald transformation is used to calculate the coulomb interactions. Approximations introduced in previous studies do alter the structural characteristics even though the internal energy differences are small. Eliminating the inverse lattice summation and using cut-offs in the potential functions result in the prediction of irregular tetrahedra and incorrect silicon-oxygen coordination behavior. These discrepancies are not present when the approximations are not used. The final state of the liquid can also be affected by the starting configuration. Structures predicted using random starting configurations were shown to be in a metastable equilibrium relative to structures derived from several crystalline starting configurations. Although this aspect of the calculations may be dependent on the potential model used, our results indicate that a two-phase procedure can be used to determine if a calculated structure is biased by the starting configuration. While not particularly important to the study of the Si02 composition, this behavior may be very important for more complicated melt compositions that do not have a cubic, crystalline analogue to use for a starting configuration. The major weakness of the ionic model is the large Si-0 bond strength which results in poor prediction of the elastic constants of the silica polymorphs. The use of a fractional charge model is an improvement over the fully ionic formalism but requires a larger number of adjustable parameters and may not be relevant to more than one composition. Improving the agreement with the elastic data will also be important in modeling the behavior of silicate melts at high pressures. We suggest that a generalized bonding model with adjustable Si-0 bond strength could be used to predict more accurately the structural and thermodynamic relationships of the polymorphs and liquid in the SiOz system. Acknowledgements-The

authors appreciate the support of D. S. Goldman, G. L. McVay, and L. R. Pederson, and assistance from T. Campbell. Critical comments on the manuscript from three anonymous reviewers are also appreciated. This work was supported by the U.S. Department of Energy.

1217

Basic Energy Sciences Office, Materials Division, under Con1830.

Itract DE-AC06-76RL0

Editorial handling: J. M. Ferry

REFERENCES ALDERB. J. and WAINWRIGHTT. E. (1959) Studies in molecular dynamics. I. General method. J. Chem. Phys. 31, 459-466.

ANASTASIOUN. and FINCHAMD. (1982) Programs for the dynamic simulation of liquids and solids, II. Rigid ions using the Ewald sum. Comp. Phys. Comm. 25, 159-I 76. ANGELLC. A., BOEHML., CHEESEMAN P. A. and TAMADD~N S. (198 I) Far IR spectra and electrical conductivity of Li and Na glasses by laboratory and computer simulation experiments. Solid State Ionics 5, 659-662. ANGELLC. A., CHEESEMANP. A. and TAMADDONS. (I 982) Pressure enhancement of ion mobilities in liquid silicates from computer simulation studies to 800 kilobars. Science 218,885-887. BACON,G. E. (1972) Coherent neutron scattering amplitudes. Acta Cryst. A28, 357-358.

BARTHT. F. W. (1932) The cristobalite structures: I. High cristobalite. Amer. J. Sci. 23, 350-356. BOCKRISJ. O’M., TOMLINSONJ. W. and WHITEJ. L. (I 956) The structure of the liquid silicates: Partial molar volumes and expansivities. Trans. Faraday Sot. 52, 299-3 IO. BREWERL. and ROSENBLAT?G. M. (1961) Dissociation energies of gaseous metal dioxides. Chem. Rev. 61, 257263. BRUSHS. G., SAHLINH. L. and TELLERE. (1966) A Monte Carlo study of a one-component plasma. J. Chem. Phys. 40,2102-2118.

CATLOWC. R. A., THOMASJ. M., PARKER S. C. and JEFFERSON D. A. (1982) Simulating silicate structures and the structural chemistry of the pyroxenoids. Nature 295, 658662.

DEMPSEYM. J. and KAWAMURAK. (1984) Molecular dynamics simulation of the structure of aluminosilicate melts. In Progress in Experimental Petrology (ed. C. M. B. HENDERSON),pp. 49-56. The National Environment Research Council. DEMPSEYM. J., KAWAMURAK., and HENDERSONC. M. B. (1984) Molecular dynamics modelling of akermanite and diopside-composition melts. In Progress in Experimental Petrology (ed. C. M. B. HENDERSON),pp. 57-59. The Natural Environment Research Council. DOANN. V. (1984) Etude en dynamique mol&culaire des defauts crCts par irradiation dans la silice vitreuse. Phi/. Msg. A49,683-687.

EWALDP. P. (192 1) Die Berechnung optischer und elektrostatischer Gitterpotentiale. Annal. Physik 64, 253-287. GAROFALINIS. H. (1982) Molecular dynamics simulation of the frequency spectrum of amorphous silica. J. Chem. Phys. 76,3189-3192.

GAROFALINIS. H. (1983) A molecular dynamics simulation of the vitreous silica surface. J. Chem. Phys. 78, 20692072.

HOSTETLERC. J. (1982)Equilibrium properties of some silicate materials: A theoretical study. Ph.D. dissertation, Univ. Arizona. KELSOJ. F., PANTANOC. G. and GAROFALINIS. H. (1983) A comparison of ion-scattering spectra and molecular dynamics simulations at the surface of silicate glasses. Su$ Sci. 134, L543-L549. KINGERY W. D., BOWENH. K., and UHLMANND. R. (1976) Introduction to Ceramics. John Wiley and Sons, 437~. LAGERG. A., JORGENSONJ. D. and ROTELLAF. J. (1982) Crystal structure and thermal expansion of a-quartz SiOz at low temperatures. J. Apply. Phys. 53, 675 l-6756. LFVIENL. and PREWITTC. T. (1981)/ High-oressure crvstal ~~ -.

121x

R. L. Erikson antri C. J. Hostetler

structure and compressibility of coesite. Amer. Mineral. 66, 324-333.

LEVIENL., PREWI~TC. T. and WEIDNERD. J. (1980) Structure and elastic properties of quartz at pressure. Amer. Mineral. 65, 920-930.

LIEBERMANR. C., RINGWOODA. E.. and MAJOR, A. (1976) Elasticity of polycrystalline stishovite. Earth Planet. .%I. Left. 32, 127-140.

MATSUIM. and MATSUMOTO T. (1985) Crystal structure and elastic constants of fi-Mg2Si04 under high pressure simulated from a potential model. Acfu Cryst. B41, 377-382. MATSUIY., KAWAMLJRA K. and SYONOY. (1982) Molecular dynamics calculations applied to silicate systems: Molten and vitreous MgSiOs and MgzSiOn under low and high pressures. Adv. Earth Planet. Sci. 12, 5 I l-524. MISAWAM., PRICE D. L. and SUZUKI K. (1980) The shortrange structure of alkali disilicate glasses by pulsed neutron total scattering. J. Non-Crystalline Solids 37, 85-97. MITRAS. K. (1982) Molecular dynamics simulation of silicon dioxide glass. Phil. Mag. B45, 529-548. MITRA S. K. and PARKERJ. M. (1984) Molecular dynamics simulation of a soda-silica glass containing fluorine. Phys Chem. Glasses 29,95-99.

MITRAS. K., AMINI M., FINCHAMD. and HOCKNEYR. W. (198 1) Molecular dynamics simulation of silicon dioxide glass. Phil. Mag. B43, 365-372. MOZZI R. L. and WARRENB. E. (1969) The structure of vitreous silica. J. Appl. Cryst. 2, 164-172. PACANSKYJ. and HERMANNK. (1978) Ab initio SCF calculations on molecular silicon dioxide. J. Chem. Phys. 69, 963-967.

PARKERS. C., CATLOWC. R. A. and CORMACKA. N. (1984) Structure prediction of silicate minerals using energy-minimization techniques. Acta Cryst. B40,200-208. PAULINGL. (1928) The influence of relative ionic sizes on the properties of ionic compounds. J. .4mer. Chem. Sot. 50, 1036-1045. PAULINGL. and PAULINC;P. (1975) Chemxtr>~. W. H. Freeman and Company. PEACORD. R. (1973) High-temperature single-crystal study of the cristobalite inversion. Zeit. Krist. 138, 274-298. POST J. E. and BURNHAMC. W. (1986) Ionic modeling of mineral structures and energies in the electron gas approximation: Ti02 polymorphs, quartz, forsterite, diopside. Amer. Mineral. 71, 142-150. PRICEG. D. and PARKERS. C. (1984) Computer simulations of the structural and physical properties of the olivine and spine1 polymorphs of Mg$iO.,. Phys. Chem. Miner& 10, 209-216. PRICE G. D., PARKER S. C. and YEOMANS J. (1985) The

energetics of polytypic structures: a computer simulation of magnesium silicate spinelloids. Acta Cryst. B41, 231239. RIEBLINGE. F. (1966) Structure of sodium alummosilicate melts containing at least 50 mole % Si02 at 15OOC.I. C’hem Phys. 44,2851-2865.

ROBIE R. A., HEMINGWAYB. S. and FISH~K J. R. (1978) Thermodynamic properties of minerals and related substances at 298.15K and 1 bar ( IO5Pascals) pressure and aI higher temperatures. U.S. Geoi. Surv. Bull. 1451’.456~. SINCLAIRW. and RINGW~OD A. E. (1978) Singie crystal analysis of the structure of stishovite. Nature 272, ‘714-7 15. SMYTHJ. R., ARTIOLIG., SMITHJ. V. and KVICK A. (1985) Crystal structure of coesite at 1SK from neutron diffraction data. EOS 66, 1116. SOULEST. F. (1979) A molecular dynamic calculation of the structure of sodium silicate glasses. J. Chem. Pt7ys 71,45704578.

SOULEST. F. and BUSBEYR. F. (198 1) Sodium diffusion in alkali silicate glass by molecular dynamics. J <‘hem !+I?.\ X,969-975.

SOULEST. F. and BUSBEYR. F. (1983) The rheological properties and fracture of a molecular dynamic simulation of sodium silicate glass. J. Chem. Phys. 78,6307-63 16. SOULEST. F. and VARSHNEYAA. K. (1981) Molecular dynamic calculations of a sodium borosilicate glass structure. J. .4mer. Cer. Sot. 64, 145-150. TOSIM. P. ( 1964) Cohesion of ionic solids in the Born model Solid State Physics 16, 1- 120. TOSI M. P. and FIJMI F. G. (1964) Ionic sizes and Born repulsive parameters in the NaCl-type alkali halides-II. The generalized Hug&s-Mayer form. .f Phys. <‘hem Solids 25,45-52.

URBAING., BOTTINGAY. and RICHETP. (1982) Viscosity ot liquid silica, silicates, and alumina-silicates. Grwhim. (‘or mochim.

Acta 46, 1061-1072.

WASEDAY (1980) The Structure of’h’on-~r!~.~tallrrr~,Uateriai; McGraw-Hill, 326~. WOODCOCKL. V., ANCELLC. PI. and CHEESEMANP. t 1976! Molecular dynamics studies of the vitreous state: Simple ionic systems and silica. J. Chem. Phys. 65, 1565- 1577. WRIGHTA. C. and SINCLAIRR. N. (1978) Neutron diffraction by vitreous silica. In The Phljsics of SiOz and it.5 Inter&w (ed. S. T. PANTELIDES),pp. 133-l 38. Pergamon Press. WRIGHTA. C. and SINCLAIRR. N. ( 1978) Neutron diffraction by vitreous silica. In The Physics of SiO, And Iis Ime$acr:$ (ed. S. T. PANTELIDES),pp. 133-138. Pergamon Press. WYCKOFFR. W. G. (1925) The crystal structure of the high temperature form of cristobalite (SiO3. Arm’,- 1 SC; 9. 448-459.