Whether the interatomic forces are calculated by an empirical potential or by ab initio techniques, the methodology for moving the ions in a molecular dynamics simulation is the same. Each method produces the force exerted on particle i by the rest of the system. The ions are large enough that quantum mechanical techniques are unnecessary - only the Newtonian equations of motion need to be integrated.
Throughout this work, periodic boundary conditions are used on a cell containing N ions. Some initial starting configuration is assumed and is calculated for each i using one of the methods discussed above. The initial positions and velocities, , are given and the equations of motion
are then integrated numerically (for example using a Runge-Kutter algorithm) over a small timestep . The size of for atomic systems is of the order of sec=1fs. This allows the positions to be updated and the process repeated. One of the more intuitive methods of finding the a minimum energy structure is to extract kinetic energy (thus eventually finding a 0K configuration) by quenching the system. That is, removing kinetic energy by stopping the ions when the force exerted on them is in the opposite direction to their velocities, . In practice, this can be slow and therefore more efficient methods are employed, as discussed below.
This method will continue until the forces are sufficiently small, where the ions are in a local minimum of energy. However, if this is the result that is required there are more efficient methods to achieve this such as steepest descents or conjugate gradients which is able to find local minima of functions (such as the energy as a function of atomic coordinates) much faster than a molecular dynamics technique.
Introducing temperature into the system is also desirable. An intuitive method for stabilising the temperature is to scale the velocities of the particles at every timestep so that the total kinetic energy of the system is given by 3/2 per degree of freedom. Although this ensures that the temperature is constant at every timestep by adding or removing energy from the system, it is not physically correct. Temperature is a statistical quantity and therefore an average over many timesteps is required to define the temperature. A method which allows the instantaneous temperature, T, to fluctuate but where <T> remains constant was suggested by Nosé. This couples the temperature to an external heat reservoir, allowing the temperature to become a dynamical variable of the system. It is analogous to attaching a `spring' to the temperature so that T can fluctuate around a mean value. This implies that once the system has reached a thermal equilibrium, the total energy remains constant although ionic, kinetic and potential energies vary.
To use the Nosé thermostat, the Newtonian equations of motion have to be integrated at each timestep requiring knowledge of the instantaneous force on each particle. As described above, calculation of the forces by ab initio methods is extremely compute intensive, whereas the use of empirical potentials requires only the evaluation of simple functions. Therefore ab initio molecular dynamics is very limited so that even for very small unit cells it is only practical to simulate a real time of 1-10ps. For this reason, empirical potentials are still very useful and are the most practical method for simulating longer timescales. A comparable simulation using an empirical potential would be able to run for the order of several tens of nanoseconds.
A Nosé thermostat is used in the finite temperature molecular dynamics simulations in later chapters and all zero temperature configurations have been found by using conjugate gradients to move the ions.
It is also possible to perform constant pressure molecular dynamics and therefore allow the size of the unit cell to change (and also fluctuate with temperature). There are fundamental difficulties in doing this by ab initio methods, which are summarised below, and therefore this method has only been implemented in the empirical calculations in Chapters 4 and 6. Allowing for fluctuations in the size of the unit cell requires the particles coordinates to be expressed as fractional coordinates, , and the parallelipiped defining the shape of the cell can be expressed in matrix form, , (the volume then being ) giving the relation
The dynamics of the boxmatrix can then be included in the Lagrangian of the total particle-boxmatrix system
where the first two terms are the Lagrangian of the particles, the third is the term for constant pressure, P, and the final term defines the kinetic energy for the box and introduces the box mass, M, which is a fictitious mass that is required in order to integrate the equations of motion for the box over a given timestep.
This allows a full relaxation of the cell shape and lattice parameters to be performed simultaneously with the particle relaxation for a given P. Unfortunately this is not easily achieved with the ab initio techniques. This is because the basis set in which the wavefunctions are expressed are plane waves whose energy depends on the size of the unit cell. Therefore changing the size of the cell will change the basis set of the wavefunctions and therefore introduce the further complication of Pulay stresses. Due to the difficulties in calculating the Pulay stresses caused by the changing basis set, a full unit cell relaxation using ab initio methods is obtained by performing many simulations at different matrices and plotting a hydrostatic curve through the phase space of given by the lowest total energy found by varying constrained to a given volume det( ). The effective pressure can then be found for a given W by the tangent to the energy volume curve.