Temperature dependence of thermodynamic quantities can be calculated from simulations performed with temperature[111, 112]. For example, the free energy of a solid can be determined from the local atomic configuration, hence a minimisation of the free energy with respect to atomic coordinates give both the equilibrium structure and free energy of the solid which could contain defects. This requires thermal statistics to be gathered in either a molecular dynamics or Monte Carlo simulation to calculate the canonical average of the thermal property.

It is not necessary to do this here since the phonon modes for a given defect concentration have been calculated directly for a given defect concentration (which is inversely proportional to the size of the supercell containing the defect). Thermodynamic properties using the vibrational spectra for the crystal can be found simply from the partition function calculated from Bose-Einstein statistics as demonstrated in Chapter 4 where a P-T phase diagram was found. The partition function allows the vibrational free energy and vibrational entropy to be calculated.

By taking the differences between the perfect crystal and the crystals containing defects (scaled to the same number of atoms) the effects that the defects have on the thermal properties can be seen. The plots of these excess free energies, total vibrational energy and vibrational entropies varying with temperature are shown in Figure 7.11.

**Figure 7.11:** Differences in (a) vibrational free energy, (b) vibrational
entropy and (c) total vibrational energy are shown as temperature
varies. The zero of each scale is taken to be the value of that quantity
of the perfect diamond structure at each temperature. (Full line -
hexagonal interstitial, dashed line - tetrahedral interstitial, dotted
line - vacancy).

These can be compared to Figure 7.12 which shows the total vibrational energy, free energy and entropy for the perfect crystal.

**Figure 7.12:** Plots of total vibrational energy, vibrational free energy and
vibrational entropy per atom against temperature. Only the plots for the
perfect crystal are shown since the differences in the curves for the
defects are too small to be seen on this scale (see previous
figure).

These quantities for the defected crystals are not shown in Figure 7.12 since they are indistinguishable on that scale. It is found that the tetrahedral interstitial has a free energy higher than that of the hexagonal interstitial at any temperature. Note that this excess free energy is times smaller than that of the defect formation energy and hence the tetrahedral interstitial is stable with respect to the hexagonal site. This leads to the conclusion that the defect formation energies for interstitials and vacancies within a crystal in the harmonic approximation are almost independent of temperature. The free energy of a silicon crystal at ambient pressure was found from the Stillinger-Weber potential[77] by summing the logarithms of the frequencies obtained from the dynamical matrix method used here. Broughton and Li[113] find the sum of the logarithms to be 4.7533 per particle, although a full free energy analysis with temperature was not given. A similar calculation using the dynamical matrix for the perfect crystal, described above, gives the value of 4.6442. Calculations on the hexagonal and tetrahedral interstitials and the vacancy give (in the same reduced units) 4.6449, 4.6452 and 4.6444 respectively. It can be seen that the introduction of defects to the crystal only changes the free energy slightly. It is the small difference in this free energy in addition to the difference in defect formation energies that determines properties such as diffusion in defects through the crystal.

The free energies of the defects can also be found, including the anharmonic effects, by running the molecular dynamics simulations at finite temperature over an extended period of time to collect thermal averages. These simulations comprised of the crystal with defect running for 12000 molecular dynamics time steps of 1.0fs which is longer than 10 periods of the lowest frequency vibrations. This allowed the evaluation of the anharmonic terms to a reasonable tolerance. It was found that the defect formation energy is unchanged within a temperature range of 0-1000K. This confirms the validity of obtaining the above results using the harmonic approximation.

Finally, the configurational entropy of the defects can be calculated to give a comparison to the vibrational entropy. The Boltzmann definition of entropy, , where is the number of configurations, can be used to calculate the configurational entropy per atom, after simplification by Stirling's formula, by

where *m* is the number of possible different defect sites and *c* is
the defect concentration. The configurational entropies for the three
defects considered are shown in Table 7.2.

**Table 7.2:** Configurational entropies for the large
supercell and the smaller supercell. The units are
in eV/K per defect.

It is found that under the harmonic approximation, the contributions to the entropy of a crystal are given by the configurational entropy and the vibrational entropy in roughly equal proportions at non-zero temperatures.

In the converged calculations the differences in the interstitial
formation energy is 0.51eV (Table 7.1). From Figure
7.11, it can be seen that the vibrational free energy is
not enough to transform a tetrahedral interstitial to the hexagonal
site within the temperature range shown. The relative stability of the
anharmonic case is tested by running a molecular dynamics simulation at
increasing temperature, allowing the bonding to change throughout the
simulation. A transition
from one defect site to another is not found until 1360K. Higher
temperatures were required to allow the interstitial atom to move
freely throughout the crystal. This also indicates that the harmonic
approximation taken here is quite accurate up to relatively high
temperatures. Recent work using the Stillinger-Weber potential on self
interstitial diffusivity in silicon[113, 105] has found a
simple migration path for several self interstitials over a temperature
range of 733-1473K. In the work by Maroudas and Brown[105] the
relaxation was done in a somewhat different manner from the molecular
dynamics simulations here. Energy minimisation was carried out by
quenching the system by a steepest descents algorithm after relaxing
the system by a lengthy Monte Carlo method at 500K with system size of
the host crystal ranging from 60 atoms to 512. The work of Maroudas and
Brown and the calculations given here suggest that the tetrahedral
interstitial configuration is the most stable (by 1.74eV and 0.51eV,
respectively) with respect to the hexagonal interstitial, but the
Stillinger Weber formation energies are both significantly larger by
about 2eV. However, they find that a self-interstitial (which they name
an extended interstitial), which has lower symmetry than either of the
hexagonal or tetrahedral sites, is an even more stable configuration.
This is similar to the *ab initio* calculations, above, where a low
symmetry configuration was also found to be stable. Their calculations
suggest a migration mechanism of the extended interstitial where the
extra atom moves from this low energy state to an identical one by
passing through the state of a tetrahedral interstitial, thus requiring
less energy than the tetrahedral-hexagonal transition that was
calculated above. The migration paths have not been examined here, but
they agree that diffusion of the tetrahedral interstitial via a path
though the hexagonal state does not occur until a much higher
temperature than the Maroudas-Brown diffusion mechanism[105].

Thu Oct 31 19:32:00 GMT 1996