Lucas, Materials Science and Engineering, University of Arizona, Tucson, AZ, USAJ. C. Mauro, Materials Science and Engineering, The Pennsylvania State University, PA 16802, USA
2.8 Atomistic Simulations of Glass Structure and Properties
Akira Takada
Research Center, Asahi Glass Co. Ltd., Hazawa‐cho, Yokohama, Japan
1 Introduction
Following Zachariasen's epoch‐making continuous random network (CRN) model [1], the first generation of atomistic models of glass relied on craft construction of ball‐and‐stick representations of atomic structures. Perhaps the most successful model was constructed for SiO2 glass by Bell and Dean [2] who patiently assembled manually rods and balls to compose 188 tetracoordinated units with a total of 614 atoms. The disordered linkage of the SiO4 groups, which were assumed to be rigid, did satisfy Zachariasen's structural rules.
From the angles and distances they measured, Bell and Dean also determined the pair distribution functions (PDF) for Si─Si, Si─O, and O─O, from which they derived a radial distribution function (RDF) that was in good agreement with the experimental data over their full range of definition, i.e. out to a distance of about five times the Si─O bond length. Yielding in particular an average O─Si─O bridging angle of 153°, this model enabled the three‐dimensional atomic configuration of SiO2 glass to be visualized, but it suffered from an obviously large arbitrariness in its handmade combination of building blocks. In addition, the model did not lend itself to estimations of physical properties. An exception was the density, which was calculated for an internal portion of the model system consisting of 72 SiO2 units, to yield a value of 1.99 g/cm3 that proved to be much lower than the actual 2.20 g/cm3.
Atomistic simulations have rendered such mechanical models obsolete as they readily provide not only atomic coordinates but also predict physical properties for glass and melts of any composition under a variety of temperature, pressure, or energy conditions. Besides, the accuracy and versatility of these calculations have been improving steadily, thanks to faster computing processors, more efficient algorithms, and bigger systems investigated. Originally, these simulations were mainly developed to give exact solutions to problems in statistical mechanics which would otherwise have been intractable for complex states of matter such as liquids, glasses, or aggregates.
The first simulations were made with the Monte‐Carlo (MC) method (e.g. [3]), which had been devised in the 1930s to solve general mathematical and statistical problems. It took advantage of the first electronic computers to sample the configurations of the system according to Boltzmann statistics, and weight them evenly when calculating the associated properties of interest. Because of this reliance on Boltzmann statistics, however, only equilibrium states can be investigated in MC simulations. At the cost of much computational complexity, the advantage of molecular dynamics (MD) simulations (e.g. [4]) is to characterize at every moment the state of the system by the positions and momenta of its constituting atoms and, thus, to account for the dynamics of the system whether in equilibrium or not. For SiO2, the first MD simulation performed by Woodcock et al. [5], for instance, dealt with the anomalous properties of the melt and assigned increases in diffusion coefficients to pressure‐induced structural changes. Following this pioneering study on a real material, a great many MD simulations have been carried out on glass/melt systems since that time.
Interatomic potentials are critical ingredients in both MC and MD simulations (e.g. [4]). As a preamble, it is thus appropriate to begin this chapter with a brief review of the manner in which they are determined within the general framework of numerical simulations. The principles of MC and MD simulations will then be presented. Finally, simulations made on amorphous oxide glasses will illustrate the interest and diversity of results that can be obtained with MD simulations and their complementarity with those of experimental studies. Amorphous SiO2, silicates, or phosphates will be mentioned but special interest will be paid to the structure of B2O3 glass that actually represents stringent tests of simulations as it requires accurate descriptions of both short‐ and medium‐range order to be understood. Other important applications of atomistic simulation to transport properties are discussed in Chapter 4.6.
2 Basics of Numerical Simulations
2.1 General Features
Atomistic simulations rely on statistical mechanical models. As such they may be performed within three main kinds of statistical ensembles. The canonical NPT ensemble (constant number of atoms, pressure, and temperature) is chiefly used in heating or quenching cycles. In contrast, the micro‐canonical NVE ensemble (constant number of atoms, volume, and energy) is primarily used when properties are calculated within the precise framework of statistical mechanics. As for the grand canonical μVT ensemble (constant chemical potential, volume, and temperature), it is typically used to investigate chemical equilibrium in systems that can exchange energy and matter with a reservoir.
For an isolated macroscopic system made up of a very large number N of atoms, each having three degrees of freedom, the microscopic state at a given instant is completely specified by the values of 3N coordinates r(i), collectively denoted by r N , and 3N momenta p(i), denoted similarly by p N . The values of the variables r N and p N define a point in a 6N‐dimensional space, called the phase space, symbolized by ΓN. If H(r N , p N ) is the Hamiltonian of the system, the path followed by this point in the phase space is determined by Hamilton's equations:
(1)
(2)
where i = 1,…,N. In principle, 6N coupled equations subjected to 6N initial conditions should be solved to specify the values of all r(i) and p(i) at a given time.
The first difficulty encountered is that the timescale of microscopic processes is ultimately controlled by the 10−14–10−15 second period of atomic vibrations. In atomistic simulation, a time integration with similar discrete steps thus is required to describe accurately the time evolution of a system. Even with today's most powerful computers, this constraint restricts simulations to low‐viscosity conditions at which relaxation times (cf. Chapter 3.7) are lower than about 10−6 seconds to be consistent with the calculational time steps.
It follows that the glass transition cannot be investigated as a number of calculation steps of the order 1018 would be required to deal with a relaxation time of ~103 seconds. Likewise, crystal nucleation and growth or liquid–liquid phase separation take place too slowly to be subjected to atomistic simulations. Besides, simulated melts are quenched at cooling rates at least six orders of magnitude faster than the highest rates of ~106 K/s achievable practically. The fictive temperatures of the simulated glasses that then be up to 1000 K higher than those of real glasses, which one should keep in mind when making any kind of comparison between both kinds of materials.
Besides, a second difficulty stems from the fact that N is of the order of the Avogadro number (6.02*1023) for any macroscopic system. Experience shows, however, that accurate results can
