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].
The first term u1 is discarded in standard simulations because it accounts for external fields (i.e. wall, electrical field, gravity, magnetic field, and centrifugal force). The second term u2 is the most important since it represents the relevant pair potentials. When determined empirically, it actually includes three‐body and many‐body effects, which is why models relying on simple pair potential model reproduce reasonably well liquid or glass structures, and why it is better in this case to denote it by the term “effective” pair potential.
As illustrated by Eq. (7), empirical potentials have a great flexibility since specific terms may be added if needed. When electrostatic interactions are important, a Coulomb charge–charge interaction may, for instance, be included in the form:
(8)
where zi, zj, and ε0 are the charge on atom i and j and the permittivity of free
