In this section, the zone centre phonon frequencies of diamond, BC8 and ST12 structure in silicon are calculated from ab initio molecular dynamics simulations at non-zero temperature. In principle, any phonon mode may be calculated in an MD simulation but very large supercells are required to give a normal mode of any general k-vector, but the molecular dynamics method has the advantage that it includes anharmonicity in the calculation. It should be noted however that there are several other methods in which to calculate phonon modes for arbitrary wavevectors such as the frozen phonon approximation[61, 62, 63] or linear response schemes[64].
The unit cells used in these calculations are the cubic 8 atom diamond
cell which will give the phonon spectral density at the 
and X
points of the primitive cell, the 16 atom cubic representation of
BC8 also giving phonons at the and X points and 12 atom
tetragonal cell for ST12 giving only the point modes.
A molecular dynamics timestep consists of relaxing the electrons to the
ground state for a fixed ionic position and then allowing the ions to
move under the influence of the Hellmann-Feynman forces, integrating the
ionic equations of motion with a timestep of 0.5fs. To ensure accuracy
of the forces the energy was converged to better than 
eV per
ionic degree of freedom, ensuring that the system is kept close to the
Born-Oppenheimer surface at all times.
The initial configuration was that of the fully relaxed structures. Each atom is given a random velocity with the ensemble average corresponding to a temperature of 500K. The average temperature was maintained close to this by using a Nosé thermostat, that is, the instantaneous temperature is allowed to vary but the time averaged value is maintained near to 500K. This thermostat introduces a coupling of the system to an external heat reservoir which has the effect of keeping the average temperature of the ions close to a fixed value. This also corrects for the errors in integrating the ionic equations of motion using a finite timestep. The unit cell dimensions are kept constant at the equilibrium lattice parameters therefore thermal expansion and fluctuations are not incorporated into the simulation.
The phonon modes can be extracted from the velocity autocorrelation function
where 
is the ionic velocity at time t and the triangular
brackets denote the ensemble average. The denominator is simply a
normalising factor. A(t) for silicon in each of the three structures
is shown in Figures 4.19, 4.20 and
4.21. It can be seen that A(t) oscillates about zero
indicating that each structure is mechanically stable.
Figure 4.19: Velocity autocorrelation function for silicon in the diamond
structure.
Figure 4.20: Velocity autocorrelation function for silicon in the BC8
structure.
Figure 4.21: Velocity autocorrelation function for silicon in the ST12
structure.
The spectral density can the be obtained simply by Fourier transforming the velocity autocorrelation function. This is shown in Figures 4.22, 4.23 and 4.23. In each case the first 0.3ps of the simulation were not included in the Fourier transform to allow for equilibration of the system.
Figure 4.22: Phonon density of states obtained by molecular dynamics for
silicon in the diamond structure. The phonon modes are those at the
and X points.
Figure 4.23: Phonon density of states obtained by molecular dynamics for
silicon in the BC8 structure. The phonon modes are those at the
and X points.
Figure 4.24: Phonon density of states obtained by molecular dynamics for
silicon in the ST12 structure. The phonon modes are those at the
point.
These figures should be compared to the phonon calculations in Chapter
4, although they show the density of states allowed in the limit of high
temperature. A direct comparison can be made only if the density of
states found from the lattice dynamics are weighted by the Boltzmann
factor, 
. Similar to the empirical calculations
it would be possible to calculate the vibrational free energy and
entropy from the density of states. A similar calculation at various
pressures would then lead to a prediction of the full
pressure-temperature phase diagram of these structures. Unfortunately,
this is not a good approximation to the vibrational properties given by
the full vibrational density of states since the vibrational
wavevectors in the calculations are extremely limited. In particular,
the zone centre modes are probably the worst sampling point to represent
the entire vibrational spectrum.
The calculation to obtain the vibrational free energy and entropy from the vibrational density of states will be fully described in Chapter 4, but for completeness, these properties are calculated here for the each of the MD calculations performed here at ambient pressure. These results are given in Figure 4.25.
Figure 4.25: Thermodynamic properties calculated from the MD simulations.
The solid line refers to the BC8 structure of silicon and the dashed
line to ST12. The zero of each quantity is that of diamond and the
quantities are given per atom in each case.
Unlike in the empirical calculations of Chapter 4 where the full phonon density of states is evaluated, the results from this limited spectral density predict that the BC8 structure has a lower free energy. At ambient pressure BC8 is also energetically more stable than ST12 silicon. This indicates that ST12 will never be stable at any temperature. When all phonons are taken into account, the opposite situation is found where ST12 silicon is found to have a greater free energy than the BC8 structure. This difference is also found to increase with temperature. This indicates, that with suitable heat treatment, ST12 will become the stable structure.
This difference in results show that one has to be careful when sampling the Brillouin zone, in that too few k-points or a badly chosen set can introduce large errors and give contradictory results.
