From the phonon calculation, the temperature dependent vibrational free energy and total vibrational energy associated with all three structures can be calculated. Since the phonon density of states has been evaluated, the partition function of the vibrations can be simply obtained using Bose-Einstein statistics from which the thermal properties of the structures in the harmonic approximation can be found. Graphs of this temperature dependence are shown in Fig. 5.12 at ambient pressure. The graphs show the total vibrational free energy, entropy and total vibrational energy per atom. On this scale, the quantities are almost indistinguishable between the three structures, but it is this difference in the free energy between the BC8 and ST12 structures that will govern the temperature at which each structure will be the most stable structure. Also shown in Figure 5.12 is the differences in the thermal properties. The zero of each scale is the value that the diamond structure takes at each temperature. Notice that on the plot of differences in vibrational free energy it can be seen that vibrational free energy for ST12 is significantly lower than that for BC8 which contrasts with the structural cohesive energy which is lower in BC8. Also note that the difference in vibrational free energy for the BC8 and ST12 structures increases with temperature, therefore whenever this difference is larger than the difference in structural energy, the ST12 structure will be favoured.
Figure 5.12: The three plots on the left show the vibrational properties of the diamond, BC8 and ST12 structures as a function of temperature at ambient pressure. The plots on the right show the differences in these quantities with the zero of each taken to be that of the diamond structure. The solid line is the difference BC8-diamond and the dashed line ST12-diamond.
Another interesting point to note from the plots of differences in thermal properties is that at low temperatures (up to approximately 200K) the diamond structure has lower vibrational free energy and entropy than the BC8 structure. This can be attributed to the fact that there is a greater density of low frequency modes in the diamond structure than in BC8. At low temperatures these low frequencies will be favoured over the more energetic bond stretching modes, thus lowering the free energy. By scaling the phonon density of states plots by it can be seen that low temperatures favour higher occupation of the lower frequency modes.
The lattice dynamics techniques described in Section 5.2 have been applied to these relaxed structures which has enabled the calculations to span the entire pressure-temperature regime in the evaluation of free energy. The total Gibbs free energy in the harmonic approximation can now be calculated by
at all points in the pressure-temperature phase space for the ST12 and BC8 structures. For a given temperature and pressure, the structure which has the lowest total free energy, G, will be favoured. The structural (cohesive) energy and PV can be obtained from either the ab initio or empirical molecular dynamics and the energy associated with the vibrations, and entropy, S, are calculated from the statistics of the lattice vibrations. Note that the temperature dependence on volume is not included here since the harmonic approximation is used in the lattice dynamics and therefore does not incorporate thermal expansion. The error in this assumption is relatively small since the volume of a tetrahedral semiconductor has a much greater response to pressure than to temperature.
Anharmonic contributions to the total free energy have been estimated from molecular dynamics simulations at a range of temperatures using the empirical potential. These results were obtained from runs with the molecular dynamics code, MOLDY, using 768 and 1024 atoms in a supercell for the ST12 and BC8 structures respectively, running for 25000 timesteps of 1.0fs. This is approximately 100 times longer than the period of the lower frequency modes therefore producing good statistics on the vibrations. Similar runs for diamond yield similar results. The absolute contribution to the energy arising from anharmonic effects was found to be negligibly small, while calculating them is extremely computationally demanding. Consequently, the harmonic approximation for the vibrational properties are used throughout.
It is possible to derive the specific heat capacity from these calculations, and in each case this was found to be Jmol , indistinguishable from the harmonic case. There is, of course, no electronic contribution to the heat capacity in this model. One useful anharmonic quantity which can be derived from the molecular dynamics is the thermal expansivity of the two structures. The behaviour of the third derivative is not included in the fit, so it is not expected that the expansivity will be very well reproduced by the potential. However, it is expected that the trend observed from the diamond phase through BC8 to ST12 to be correct. The values are found to be and for BC8 and ST12 respectively. These are somewhat smaller than the value for diamond of (about 20% higher than experiment). As yet there has not been any experimental measurement of expansivity in the metastable phases.
The empirical model predicts that the difference in structural energy, , to be smaller than that of the ab initio calculations for silicon. This will affect the position of the stability field boundary in the phase diagram. For internal consistency the empirical potential was used throughout to calculate a full pressure-volume phase diagram for ST12 and BC8 silicon which is shown in Fig. 5.13. Note that this is a metastable phase diagram since the diamond structure is stable throughout the P-T region shown.
Figure 5.13: Relative stability fields of the ST12 and BC8 structures as determined from the empirical potential using the vibrational free energies as a function of temperature and the structures found under pressure from the molecular dynamics simulations.
The main contribution to the vibrational free energy comes from the low vibrational frequencies. These are mostly the bond bending modes which are determined by the C parameter in the final term of the potential. Reduction in the C parameter softens these modes and therefore increases the vibrational free energy. Hence for lower C, the minimum temperature at which the ST12 structure is favourable is lowered. The first two terms in the potential essentially determine the depth, curvature and position of the minimum of the structural energy of a tetrahedrally bonded crystal, therefore a simple rescaling could be used to give a germanium-like potential which in effect would weaken the bond stretching and lower the corresponding phonon frequencies. As noted above, this would retain the general shape of the phase diagram, but lower the temperature required to favour the ST12 structure. However, the bond bending term in germanium will be even more significantly reduced and therefore reduce the low bond bending modes even further, an effect which can be traced back to the smaller overlap of orbitals on distortion. This also lowers the temperature and the pressure at which the ST12 structure will become stable over the BC8 structure. It makes ST12 germanium the metastable depressurised structure at room temperature. It therefore seems likely to be the case that BC8 germanium will be favoured only at low temperatures and low pressures, perhaps forming on rapid depressurisation. Recently this was found to be the case where BC8 germanium was recovered from rapid depressurisation from the metallic phase. The effects of temperature were not reported but at room temperature BC8 germanium was found to be short lived, although increasing the pressure significantly increased its lifetime.