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 be obtained for systems of only a few hundred or thousand atoms as long as structural units bigger than the system itself are not involved in the processes investigated. To avoid either creating surface or setting up a wall (Figure 1), periodic boundary conditions are usually imposed in numerical simulations. The cubic box is replicated throughout space to form an infinite lattice. As an atom moves in the original box, its periodic image in each of the neighboring boxes moves in exactly the same way. As an atom leaves the central box, one of its images will enter through the opposite face.
To simulate accurately a noncrystalline configuration in the cell, a larger number of atoms are nonetheless preferable. When a small cell of around several hundred atoms is studied, periodic boundary conditions, for instance, prevent structural units larger than the cell from forming. As already alluded to, any structural unit or spatial fluctuation larger than the wavelength which is greater than half of the cell size cannot be calculated appropriately because the repeatedly arranged fragments would be involved in the results. On the other hand, the maximum number of atoms which can be currently processed is of the order of 105–106. To check whether the number of atoms or the size of simulation cell is sufficient or not, the best way is to double the number of atoms and to maintain the same density through an adjustment of the cell size and then to make sure that there is little change in the calculated structural data.
Figure 1 Schematic representation of periodic boundary conditions.
2.2 The Importance of Interatomic Potentials
Strictly speaking, in numerical simulations the Hamiltonian should be calculated from the quantum wave‐function equation, Hψ = Eψ, where ψ and E are the wave‐function and energy of the system, respectively. As expounded in Chapter 2.9, such calculations require so much computing work that they are currently restricted to smaller systems typically made up of a few tens of atoms. In the present chapter, simulations made within a classical framework will thus be considered instead. They rely on the fact that the differences between vibrational and electronic energies and frequencies are so large that atomic vibrations may be considered to take place within a fixed electronic configuration. This is the celebrated Born–Oppenheimer approximation whereby the Hamiltonian of a system is expressed as the sum of kinetic and potential energies, which are functions of the selected set of coordinates q(i) and momenta p(i):
(3)
In Eq. (3) the kinetic energy, K(p), is simply expressed as a function of the mass and of the velocity calculated for each atom. The potential energy, Up(q), is much less readily determined because it strongly depends on the specific interactions between the various kinds of atoms present in the system. Quite generally, however, interaction potentials are markedly asymmetrical in terms of interatomic distances because they rise extremely steeply when atoms get mutually very close whereas they reach a constant value – the dissociation energy – when atoms become so distant that they can be considered as no longer bonded. As introduced in 1929 to represent isolated diatomic molecules, the Morse potential accounts well for this asymmetry:
(4)
where E0, k, and r0 are the dissociation energy, a measure of the bond strength, and the equilibrium interatomic distance of the molecule, respectively, three parameters that are determined from vibrational spectroscopy data.
In a condensed phase, potentials are much more complicated since a given atom interacts with a great many others over distances that can be large. To keep the number of parameters as small as possible in the expression of potential energies, one thus groups into the same term all interactions between given pairs of like or unlike atoms regardless of their mutual distances. Although the Morse potential remains a good starting point for systems where bonding is covalent, other kinds of analytical expressions are generally used for potential energies in the MC and MD simulations dealt with in this chapter. As borne out by the variations with composition of macroscopic properties, atomic interactions have the simplifying feature that they are primarily pairwise in oxide or salt systems. This feature is embodied in the most popular potentials used for these systems, namely, the Buckingham,
(5)
and the Born–Mayer–Huggins potentials,
(6)
where rij is the distance between two atoms, and A, B, r0, C, D are parameters inherent to each interaction (Figure 2).
In both potentials, the first, second, and third terms represent repulsive interaction, dipole–dipole dispersion, and dipole–quadrupole dispersion, respectively. In more precise formulations, ternary and higher‐order effects must be accounted for so that the potential energy is made up of terms depending on the coordinates of individual atoms, pairs, triplets, etc.
Figure 2 Examples of potential energy models: Morse and Buckingham potentials used in B2O3 simulations for B─O and for O─O and B─B, respectively [6].