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.
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.