Molecular-dynamics simulations of (NaO2)x(SiO2)1-x glasses: Relation between distribution and diffusive behavior of Na atoms

Molecular-dynamics simulations of „NaO2…x„SiO2…12xglasses: Relation between distribution and diffusive behavior of Na atoms Jaime Oviedo and Javier Ferna ´ndez Sanz* Departamento de Quı ´mica Fı ´sica, Facultad de Quı ´mica, E-41012 Sevilla, Spain ~Received 7 April 1998! Molecular-dynamics simulations of sodium silicate glasses in a range of alkali concentration going from 1.8% to 33.33% molar of Na2O are reported. Our simulations show that there is a tight relationship between Na atoms and nonbridging oxygens which are mainly located in the first Na coordination shell. In the whole range of composition, Na and nonbridging oxygen atoms appear to segregate giving rise to a heterogeneous distribution. For the higher alkali concentrations, formation of microchannels is observed. The activation energies for alkali diffusion have been computed and found in agreement with the experiment. The mechanism for the diffusion has also been investigated and found to occur through the microchannels, if present, or across the network for low compositions, but the Na motions always appear to be assisted by nonbridging oxygen atoms. @S0163-1829~98!09637-4# I. INTRODUCTION Silicate glasses constitute one of the most widely used materials. However, in spite of its importance, the structure of modified oxide glasses is not well known yet. Considerable efforts have been devoted to elucidate this by using new techniques, such as x-ray absorption fine-structure spectroscopy ~XAFS!,1–3 magic angle spinning NMR ~MAS-NMR!,4–7 isotopic substitution neutron scattering,8as well as the traditional spectroscopies of IR,9Raman,10 and x-ray photoemission spectroscopy ~XPS!.11,12 However, the experimental data are not conclusive, and from these studies both nonrandom13,14 and random distributions of modifier cations15 have been suggested. It is worth noting that neither inhomogeneous nor homogeneous distributions have been directly observed, the atomic arrangement depending on the interpretation of the experimental data. With respect to transport properties, it has been shown that the activation energy for diffusion in these glasses decreases drastically from '1.2 eV at low alkali concentrations to about 0.6 eV at around 15 mol % of alkali oxide. Above 15 mol % of alkali oxide, it decreases very little when the alkali oxide concentration is increased up to 33 mol %.16 Also, the diffusion mechanism associated to the alkali ion transport has been the subject of several conjectures. Recently Greaves13 suggested a modified random network ~MRN!model for silicate glasses to account for the structural and transport properties. This model results from XAFS data in which the number of nonbridging oxygens ~NBO!around alkali ions was found to be larger than that of bridging oxygens ~BO!.1–3 In this model, there is a process of microsegregation in the glass and the modifying cations form clusters in the SiO2rich matrix, to which they are predominantly linked via the NBO’s. At higher concentrations, the clusters percolate through the glass bulk giving rise to pathways for ionic conduction. These suggestions are, as yet, unconfirmed. Molecular dynamics ~MD!simulations are a valuable tool for insight into glass structure and to help interpret experimental results. Previous MD studies have reported both homogeneous17 and heterogeneous18–20 distribution of alkali cations. There have also been studies of transport properties and, thus, Huang and Cormack reported on the diffusion coefficients of Na for a range of concentrations. However, these simulations were constrained since in the calculations the silicate network was frozen and, in fact, it was not until recently that Smith, Greaves, and Gillan21 reported quantitative results for the sodium diffusion coefficient. It has been suggested that the alkali ion motion is characterized by jumps from site to site in the silicate matrix.21,22 However, these studies have been limited to compositions with a high alkali concentration. The purpose of this paper is to report the results of an investigation of the structure and dynamics of sodium silicate glasses (Na2O)x(SiO2)12xin a wide range of concentrations and temperatures through unconstrained MD simulations. Our first aim is to analyze the effect of Na doping on silicate structures as concentrations change. At the same time, a quantitative determination of diffusion constants and activation energies as well as the diffusion mechanisms over this range of concentrations will be reported. These simulations will allow us to test previous conjectures about alkali transport mechanism. II. COMPUTATIONAL DETAILS As is well known, potential functions play a key role in MD simulations. In the present work the interatomic potentials employed for SiO2were those derived by Vashishta et al.23 This potential includes a three-body term to take into account the Si tendency towards tetrahedral coordination and has been shown to give good results for vitreous and molten states. The reported effective charges were 20.8 and 1.6 for O and Si, respectively, although, in order to maintain electroneutrality, it is necessary to vary these charges with Na concentration. For this purpose we have chosen to assign charge 11 to Na cation and to vary only the O charge to keep our system neutral. The interaction for the Na-Si, Na-O, and Na-Na pairs has been described by means of a simple pair potential consisting PHYSICAL REVIEW B 1 OCTOBER 1998-IIVOLUME 58, NUMBER 14 PRB 58 0163-1829/98/58~14!/9047~7!/$15.00 9047 © 1998 The American Physical Society of an electrostatic contribution and a short range repulsive term. The short-range term was represented through a Born- Mayer analytical potential: Vij5Aije2rij / r ij1qiqj rij , whose parameters Aand r for the Na-Na, Na-O, and Na-Si pairs were adjusted in order to reproduce the first coordination shell distance @2.3–2.4 Å ~Refs. 1, 3, and 24–26!# observed in the radial distribution function ~RDF!obtained from MD simulations for a system with a Na2O concentration of 10.74% molar ~system C in Table I!. Although this choice could be considered arbitrary, we hope it will be valid to perform a systematic study of vitreous states in a Na2O concentration range. On the other hand, after a preliminary analysis it was found that the short range term in the Na-Na and Na-Si potentials were negligible and, therefore, a good description of the system could be achieved considering only the parameters for Na-O interaction pair, whose final values were: A(Na-O)53060310212 erg and r (Na-O)50.265 Å. The simulations were performed for five Na2O concentrations ranging between 1.8% M ~system A!and 33.33% M ~system E!. The experimental density at room temperature for every Na2O concentration was used and for the sake of simplicity they were assumed to be constant in the range of temperatures considered. The main physical aspects of the system are summarized in Table I. MD simulations were undertaken in the microcanonical ensemble, with constant NVE, where N is the number of atoms, V is the volume, and E is the total energy of the system. Equations of motion were numerically integrated using the leapfrog algorithm with a time step of 0.5 fs. Periodic boundary conditions were handled according to the minimal image and the electrostatic contributions were computed through the Ewald sum technique. During our calculations energy stays constant better than one part in 105. For a given alkali concentration, the initial configuration was obtained by randomly adding the appropriate amount of Na2O to vitreous SiO2. This configuration was thermalized at 3000 K in a 30 000 steps simulation in which the velocities were rescaled. After this thermalization, another run of 5000 steps was performed in order to check that the system reached equilibrium, and finally, a production run of 15 000 steps ~7.5 ps!was carried out for statistics. From this initial state at 3000 K we generated others at 2500, 2000, 1500, 1000, 600, and 300 K. For every temperature the same simulation runs schedule was used. From these simulations we obtained the static properties; however, we found that the statistics was relatively poor for the calculation of the diffusion coefficients, mainly for the systems with low amounts of Na2O. For this reason, we increased the production runs up to 45 ps for the temperatures higher than 1500 K ~there is no appreciable diffusion below this temperature!. Also, additional states at 3500 and 4000 K were created. III. RESULTS AND DISCUSSION Structure We will start this section analyzing structural results. In general, the radial distribution functions, RDF’s, corresponding to the pairs Si-O, Si-Si, and O-O are almost the same as those found for pure SiO2. On going from system A to E, there is a small shift of the Si-O RDF first maximum from 1.61 to 1.58 Å. This shift can be easily understood taking into account that increasing the Na concentration leads to more negative charge on the oxygen atoms, and therefore to an increment in the Coulombic term of the interaction potential. FIG. 1. Radial distribution functions for Na-O, Na-Si and Na-Na pairs. TABLE I. Physical description of (NaO2)x(SiO2)12xsamples. System ~concentration % Na2O molar! A~1.8%!B~4.8%!C~10.74%!D~19.4%!E~33.33%! Atoms 4 Na2O–216 SiO2 11 Na2O– 216 SiO2 26 Na2O– 216 SiO2 52 Na2O– 216 SiO2 108 Na2O– 216 SiO2 ~g/cm3!2.21 2.24 2.31 2.40 2.56 Box side ~Å!21.499 21.634 21.886 22.379 23.350 qo20.810 20.850 20.913 21.015 21.200 9048 PRB 58 JAIME OVIEDO AND JAVIER FERNA ´NDEZ SANZ Concerning the structure of Na sites, the RDF’s for Na-O, Na-Si and Na-Na obtained at 300 K are reported in Fig. 1, while the first maxima of the distributions are summarized in Table II. As can be seen, the first coordination shell for Na is well defined and, in general, the maximum of the Na-O RDF slightly shifts towards lower distances when doping. This is not the case for system A, whose first maximum in the Na-O distribution appears at shorter distance as a well resolved peak. This feature will be explained later. With respect to the values obtained, they fall in the range 2.3–2.4 Å showing the coherency of the fitting procedure used in the determination of the Born-Mayer potential parameters. On the other hand, when temperature rises, a thermal broadening of the RDF’s is observed and, as shown in Table II, the maxima position shifts slightly toward lower values. However, the shape of the bands persists, indicating that the cation environment is maintained. Similar behavior is found in the Na-Si RDF’s, although the shifts induced by both temperature and alkali concentration are somewhat larger. With respect to Na-Na pair, the RDF is reported at the bottom of Fig. 1 where there is a prominent band indicating a clear coordination shell in the 3.5–4.1 Å distance range. This result could appear amazing, since assuming a homogeneous distribution, the first Na-Na maximum would be at larger distance ~for example, 5.86 Å for system C and even larger for systems A and B!. This means that Na cations tend to be together and there is a spatial heterogeneous distribution. This result supports the modified random network model proposed by Greaves13 based on XAFS experiments carried out with a disilicate composition sample. Also, our value for the Na-Na distance corresponding to the same composition ~system E, 3.59 Å!agrees reasonably with the 3.2 Å distance obtained from the XAFS data. From the RDF’s, and taking as cutoff the distance corresponding to the first minimum, we have estimated the coordination numbers for Si and O atoms. These values for the A-E compositions are reported in Table III for temperatures 300 and 3000 K. This analysis shows that, whatever the composition is, most of the silicon atoms are fourcoordinated, the number of three-coordinated silicons increasing moderately when the temperature rises. With respect to the oxygen coordination number, it appears clearly that when the Na concentration increases, the number of bridging oxygen ~BO!atoms ~i.e., those bound to two silicon atoms!decreases monotonously, while, concomitantly, the number of nonbridging oxygen ~NBO!atoms rises. Roughly, it is found that there is one NBO by each Na atom. On the other hand, the temperature dependence is low and the general trend is an increment of under-coordinated species when temperature increases. The coordination environment of Na cation can be observed in Fig. 2, where the coordination number Nagainst the distance is reported for the Na-O and Na-Na pairs. The Na-O curves show an inflection of about N55,6 for all compositions except that of system A, for which N52,3. Such a fivefold, sixfold coordination agrees well with the XAFS data. Moreover, there is not a regular pattern for the coordination polyhedra in disagreement with the formation of trigonal bipyramidal units around the Na cation suggested by Smith, Greaves, and Gillan.21 On the other hand, the Na-Na curves show a regular increment of the coordination number as Na concentration rises. In order to make clear the relationship between NBO’s and Na cations, an analysis of the type of oxygen atoms bound to the Na cations was performed. In Table IV, the oxygen atoms have been grouped depending on whether they are of NBO or BO type, and within each class whether they are bound or not to a Na cation. Inspection of Table IV shows that, as expected, increasing the alkali concentration gives rise to an increment of both BO and NBO percentages TABLE II. First maximum ~Å!in the RDF at 300 K and 3000 K ~between parentheses!. System ABCDE Na-O 2.31 ~2.23!2.40 ~2.39!2.38 ~2.39!2.38 ~2.36!2.37 ~2.32! Na-Si 3.76 ~3.69!3.56 ~3.57!3.54 ~3.45!3.41 ~3.33!3.01 ~3.08! Na-Na 3.66 ~4?!3.55 ~4.16!4.16 ~4.04!3.98 ~3.95!3.59 ~3.53! TABLE III. Coordination ~%!numbers for Si and O atoms at 300 K ~3000 K!. Coord. A B C D E 3 3.21 5.07 3.66 0.20 0.45 Sin~10.09!~11.42!~8.94!~4.29!~3.22! 4 96.79 94.93 96.34 99.8 99.5 ~89.86!~88.52!~91.00!~95.69!~96.77! 1 3.44 8.58 14.85 21.69 34.07 On~7.34!~11.21!~18.01!~22.35!~35.09! 2 96.54 90.28 82.73 77.48 62.59 ~92.06!~87.72!~80.52!~76.35!~60.57! 3 0.02 1.14 2.20 0.62 0.19 ~0.57!~0.94!~0.70!~0.52!~0.73! PRB 58 9049MOLECULAR-DYNAMICS SIMULATIONS OF... related to Na cations. However, whatever the composition is, most of the NBO’s are bound to Na centers ~last column!. This tight relationship between NBO’s and Na cations appears to be unaffected by the lowering of the cutoff distance as shown in the table, where results for cutoffs of 3.3, 2.9, and 2.4 Å are reported for C composition. Contrarily, the number of BO’s bound to Na centers drops significantly when the cutoff is lowered. In summary, from this analysis and the Na-Na RDF’s we can conclude that both Na and NBO’s tend to be grouped into the bulk of the SiO2matrix, which is typical of a microsegregation phenomenon. This behavior can easily be seen in Fig. 3 where snapshots of the MD simulations for systems A, C, and E are presented. The trend for the Na and NBO atoms to be grouped appears clearly in the three images. Also, it is found that above the percolation limit composition ~system C!, the microsegregation process gives rise to formation of microchannels. In order to get a deeper insight of the structure, the oxygen contributions to the RDF arising from NBO’s and BO’s has been separated and are reported in Figs. 4 and 5. The first peak of the NBO-NBO RDF ~Fig. 4, top!appears at 2.60 Å and corresponds to two NBO’s bound to the same Si ~this peak is not present in system A because such a situation is unusual at very low Na concentration!. More interesting is the second peak, approximately at 3.3 Å which has been attributed to the presence of trigonal bipyramids.21 The analysis of the data reveals that it corresponds to O-Na-O units belonging to highly distorted polyhedra. On the other hand, the RDF for BO-BO pairs ~Fig. 4, bottom!are found to be similar to the total O-O RDF. Nevertheless, when Na concentration is high ~systems D and E!, the peak at 4.9 Å disappears and the curves show again a band at 3.3 Å indicating the formation of BO-Na-BO units. With respect to NBO-Na and BO-Na RDF’s shown in Fig. 5, it turns out that the Na cations bind NBO’s giving rise to a well defined band, whatever the concentration is ~2.30, 2.33, 2.34, 2.30, and 2.28 Å for systems A–E, respectively!. For BO’s the first peak becomes clear only at high concentrations ~2.30 and 2.28 Å for systems D and E!, the formation of this peak being coherent with the presence of the 3.3 Å peak at BO-BO RDF for systems D and E. On the other hand, and in light of these partial contributions of the oxygen atoms to the RDF, a new analysis of the coordination numbers could be performed. However, it is difficult to assign a cutoff distance to every peak and cutoffs are different for NBO and BO contributions, although we can conclude that Na cations are bound to 2–3 NBO’s and complete their coordination up to 5–6 with BO’s. The combination of NBO-Na and BO-Na RDF’s allows us to explain the good resolution observed for the first peak in the O-Na RDF for system A in Fig. 1, as FIG. 2. Plot of the Na coordination number N(r) as a function of the distance. FIG. 3. Snapshots for some compositions. Top: system A. Middle: system C. Bottom: system E. TABLE IV. Proportion of bridging ~BO!and nonbridging oxygen atoms in % at 300 K ~3000 K for C8!. Oxygen type BO BO ~Na!NBO NBO ~Na! A 95.61 0.95 0.25 3.19 B 81.81 9.62 0.01 8.56 C(r53.3 Å) 50.42 34.51 0.01 15.06 D 13.13 64.97 0.00 21.90 E 0.80 61.97 0.00 37.22 C(r52.9 Å) 65.13 19.81 0.12 14.94 C(r52.4 Å) 83.41 1.53 5.17 9.89 C8(r53.3 Å) 48.96 34.47 0.83 15.74 9050 PRB 58 JAIME OVIEDO AND JAVIER FERNA ´NDEZ SANZ well as its low coordination number, since for this system only the NBO-Na pair makes a significant contribution to the total RDF. A further structural property refers to the angular distribution functions ~ADF!. The ADF for O-Si-O always shows a sharp peak at 109° reflecting that the tetrahedral coordination of Si atoms is preserved when composition is varied. For the Si-O-Si ADF a broad distribution is found, and on going from system A to E, the maximum shifts from 144° for system A, the same value as that of pure SiO2,27 to 142°, 140°, 137°, and 136° for systems B, C, D, and E, respectively. This lowering indicates that the presence of the O-Na-O triads slightly distorts the relative position of the SiO4tetrahedra. The ADF’s associated with the Na atoms appear to be more interesting, and are reported in Fig. 6 for systems A, C, and E at 300 K. The distributions for O-Na-O show a first well resolved band around 56°, followed by a broad band with maxima around 90° and 115°. The presence of the peak at 56° confirms the absence of trigonal bipyramids. A careful analysis of Na environments shows that the first peak is FIG. 4. Radial distribution functions for NBO-NBO and BO-BO pairs. FIG. 5. Radial distribution functions for Na-NBO and Na-BO pairs. FIG. 6. Angular distribution functions for Si-O-Na and O-Na-O triads at different compositions. PRB 58 9051MOLECULAR-DYNAMICS SIMULATIONS OF... mainly associated to situations in which Na is close to a BO and a NBO belonging to the same SiO4 tetrahedron ~NBO- Na-NBO triads are much less common and its associated angle is 70°!. On the other hand, the broad bands correspond to oxygen pairs belonging to different Si atoms. It should be noticed that the first maximum, ;90° is in agreement with a Na-O distance of 2.3 Å and the maximum observed at 3.3 Å in the O-O RDF. In its turn, the Si-O-Na ADF shows a broad band in the 80°–120° range that becomes narrower as Na concentration rises. We have also analyzed ADF’s with more than one Na atom, although only those of Na-O-Na and Na-Na-Na will be briefly considered. For the Na-O-Na a broad band with a maximum at ;95° is found. Taking 2.3 Å for the Na-O distance, this angle gives a Na-Na interdistance of 3.4 Å in agreement with the maximum observed in the Na-Na RDF. The Na-Na-Na ADF features a poorly resolved distribution which broadens when temperature rises. In general, repetitive Na-Na-Na patterns are not found in the present simulations. Na ion diffusion The diffusion coefficients Dfor systems A–E have been determined in the @2000–4000 K#temperature range from the plots vs time of the medium square displacements ~MSD!. The variation of logarithm of Dagainst 1/Tis plotted in Fig. 7. Inspection of this figure shows a quasilinear behavior in agreement with the Arrhenius law observed21 for diffusion coefficients for alkalis in disilicates. According to this, the diffusion coefficient is given by: D5D0eEa/kBT where D0is the preexponential factor, kBis the Boltzmann constant and Eais the activation energy. The slope of a least-squares linear fit to our data gives the following activation energies for the compositions here considered: 1.70, 0.99, 0.85, 0.84, and 0.66 eV for systems A, B, C, D, and E, respectively. When these values are compared with the experiment an excellent agreement is found. First, a lowering of the energy activation is found and this lowering appears to occur quickly. For Na concentration larger than that of B, there is a relative stabilization, the activation energy diminishing more slowly until the disilicate composition, system E, for which the experimental value is of 0.68 eV.16 Our values are also in good agreement with those obtained from previous MD simulations for the disilicate composition ~0.56 eV!.21 On the other hand, at the level of impurities, the activation energy has been estimated to be about 1.2–1.3 eV.16 Likely, our value for system A could be overestimated, since as shown in Fig. 7, only three points have been used in the fitting. This is because at lower temperatures, no significant diffusion is observed. In any case, it clearly appears that the slope for these three points is larger than those of higher concentrations. In summary, our calculated values compare well with both experimental and theoretical measurements, and, as far as we can ascertain, this is the first time that these diffusion coefficients have been calculated in this range of concentrations. The mechanism for Na diffusion will be considered now. In Fig. 8, typical trajectories for some compositions and temperatures are depicted. In the top of this figure, the trajectories for low concentration ~systems A and B!at 2500 K describe a quasicontinuos motion of the Na cation through the network. This type of motion is observed in the whole temperature range and would agree with the ‘‘network hopping’’ mechanism proposed by Greaves and Ngai28 for low alkali concentration, in which it was suggested that the hopping would occur without the participation of NBO. However, the analysis of the coordination shell for Na along these trajectories shows that when Na moves the description of its environment is essentially the same, i.e., new NBO’s are being formed along the Na motion through breaking of Si-O bonds. On the other hand, for high alkali concentrations ~C and higher!, the trajectories depicted in the bottom part of Fig. 8 are indicative of sudden jumps between sites. Inspection of the Na environment along these trajectories reveals that the Na motion appears to occur through the microchannels present in the material. This mechanism is in agreement with the ‘‘intrachannel hopping’’ proposed by Greaves and Ngai28 and confirmed later by Smith, Greaves, and Gillan21 from MD simulations of disilicate. It is worth noting that in our simulations, these kind of trajectories are observed only at temperatures below 2500 K. When the temperature rises, patterns similar to those of systems A and B are found. This FIG. 7. Plot of the diffusion coefficient logarithm ln Dagainst T21.FIG. 8. Representation of Na trajectories for different compositions. 9052 PRB 58 JAIME OVIEDO AND JAVIER FERNA ´NDEZ SANZ behavior does not mean that there is a change in the mechanism, which still goes through the microchannels, but it reflects simply that the residence time in the sites where Na cation is stabilized is much shorter, according to the larger values of the diffusion coefficients. Finally, it should be noted that the diffusion behavior is perfectly coherent with the computed activation energies. For low alkali concentrations, Na motion necessarily involves breaking of Si-O bonds giving rise to NBO’s, and the energetic cost would account for the higher activation energy for diffusion. However, for alkali concentrations higher than the percolation limit, the Na diffusion occurs assisted by the NBO’s already present in the microchannels. The motion is easier, with lower energetic barriers, and becomes significantly faster when the kinetic energy of the cations reaches a given threshold value. IV. CONCLUSIONS In the present work MD simulations of sodium silicate glasses in a range of alkali concentrations and temperatures are reported. From the study of several structural properties and direct inspection of configurations from MD simulations the following points can be concluded: ~i!on increasing the Na2O concentration, an increment in the number of NBO is found. Roughly, there is one NBO for each Na atom; ~ii! NBO’s are found to be close to the Na atoms, belonging to their first coordination shell; ~iii!the Na atoms ~and therefore the NBO’s!tend to be together. These three findings are characteristic of a microsegregation process in which, starting from a homogeneous alkali distribution, the system evolves towards a heterogeneous one. Beyond a given Na2O concentration ~.10% in this simulation!, formation of microchannels is also observed. These results would support the modified random network model proposed by Greaves.13 Activation energies for alkali diffusion have been computed and found in quite good agreement with the experimental ones, both in their magnitude and in their dependence on the Na2O concentration. The alkali diffusion always is assisted by NBO’s and for high concentrations it takes place through the microchannels. In such a case it is relatively easy since NBO’s are always present in these regions. For low concentrations, the alkali diffusion involves breaking of Si-O bonds, in agreement with the larger activation energy. ACKNOWLEDGMENTS This work was supported by the DGICYT ~Spain, Project No. PB95-1247!, and by the European Commission ~Contract No. ERBCT1-CT94-0064!. J.O. thanks the Junta de Andalucı ´a for financial support. *Corresponding author. Electronic address: [email protected] 1G. N. Greaves, A. Fontaine, P. Lagarde, D. Raoux, and S. J. Gurman, Nature ~London!293, 611 ~1981!. 2G. N. Greaves, in Glass Science and Technology, edited by D. R. Uhlmann and N. J. Kreidl ~Academic, London, 1990!, Vol. 4A, p. 1. 3S. N. Houde-Walter, J. M. Inman, A. J. Dent, and G. N. Greaves, J. Phys. Chem. 97, 9330 ~1993!. 4R. Dupree, D. Holland, and M. G. Mortuza, J. Non-Cryst. Solids 116, 148 ~1990!. 5H. Maekawa, T. Maekawa, A. Kawamura, and T. Yokakawa, J. Non-Cryst. Solids 127,53~1991!. 6W. Hater, W. Mu ¨ller-Warmuth, M. Meier, and G. H. Frischat, J. Non-Cryst. Solids 113, 210 ~1989!. 7I. Farnan, P. J. Grandinetti, J. H. Baltisberger, J. F. Stebbins, U. Werner, M. A. Eastman, and A. Pines, Nature ~London!358,21 ~1992!. 8P. H. Gaskell, in The Physics of Non-Crystalline Solids, edited by L. D. Pye, W. C. LaCourse, and H. J. Stevens ~Taylor & Francis, London, 1992!,p.15. 9E. I. Kamitsos, A. P. Patsis, and G. D. Chryssikos, in The Physics of Non-Crystalline Solids ~Ref. 8!, p. 461. 10N. Umesaki, N. Iwamoto, M. Tatsumisago, and T. Minami, J. Non-Cryst. Solids 106,77~1988!. 11R. Bru ¨ckner, H. U. Chum, H. Goretzki, and M. Sammet, J. Non- Cryst. Solids 42,49~1980!. 12D. Sprenger, H. Bach, W. Meisel, and P. Gu ¨tlich, J. Non-Cryst. Solids 159, 187 ~1993!. 13G. N. Greaves, J. Non-Cryst. Solids 71, 203 ~1985!. 14J. F. Stebbins, J. Non-Cryst. Solids 106, 359 ~1988!. 15R. Gua ˚ker and S. Urnes, Phys. Chem. Glasses 14,21~1973!. 16G. H. Frischat, in Ionic Diffusion in Oxide Glasses ~Trans. Tech., Aedermannsdorf, 1975!. 17A. A. Tesar and A. K. Varshneya, J. Chem. Phys. 87, 2986 ~1987!. 18C. Huang and A. N. Cormack, J. Chem. Phys. 93, 8180 ~1990!; 95, 3634 ~1991!. 19C. Huang and A. N. Cormack, in The Physics of Non-Crystalline Solids ~Ref. 8!. 20B. Vessal, G. N. Greaves, P. T. Marten, A. V. Chadwick, R. Mole, and S. Houde-Walter, Nature ~London!356, 504 ~1992!. 21W. Smith, G. N. Greaves, and M. Gillan, J. Chem. Phys. 103, 3091 ~1995!. 22T. F. Soules and R. F. Busbey, J. Chem. Phys. 75, 969 ~1981!. 23P. Vashishta, R. K. Kalia, J. P. Rino, and I. Ebbsjo ¨, Phys. Rev. B 41, 12 197 ~1990!. 24B. E. Warren and J. Biscoe, J. Am. Ceram. Soc. 21, 259 ~1938!. 25G. N. Greaves and K. L. Ngai, in Proceedings of the International Conference on Defects in Insulating Materials, edited by O. Kanert and M. Spaeth ~World Scientific, Singapore, 1993!,p. 53. 26M. Misawa, D. C. Price, and K. Suzuki, J. Non-Cryst. Solids 37, 85 ~1980!. 27R. L. Mozzi and B. E. Warren, J. Appl. Crystallogr. 2, 164 ~1969!. 28G. N. Greaves and K. L. Ngai, Phys. Rev. B 52, 6358 ~1995!. PRB 58 9053MOLECULAR-DYNAMICS SIMULATIONS OF...