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é[27]. 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[28]

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.

Thu Oct 31 19:32:00 GMT 1996