There are several ways of calculating the elastic properties of a material from the knowledge of the interaction potential. They can be found either from gathering statistics from long molecular dynamics calculations using the fluctuation formula method, or directly from inverting the matrix of spring constants which can easily be calculated from the interaction potential. The advantage of using the molecular dynamics method is that the elastic constants can be calculated at temperature and the other thermodynamic averages can be calculated simultaneously. Obviously, inverting a matrix at every timestep in the simulation of temperature is computationally prohibitive, but the calculations presented here are carried out at 0K, therefore making the latter method more effective since the dynamical matrix has already been calculated.
A full derivation of the calculation of elastic coefficients can be found in reference . The final method for a crystal containing many atoms in a single primitive unit cell will be given here. The elastic coefficients can be calculated directly from as follows. Allow , , , , and to run over the co-ordinate axes x, y and z. Then define
There will be three modes that have zero frequency due to the three translational degrees of freedom of the entire crystal. Therefore will be singular since it will have linearly dependant rows and columns. Now we introduce the matrix which is defined to be the inverse of , where we allow , and for convenience of calculation define
For convenience, the k=1 rows and columns are removed from the dynamical matrix for the inversion. The choice of k is arbitrary and therefore can be taken to be k=1 without loss of generality.
Now define the following brackets where is the volume of the supercell:
The elastic coefficients are then
provided the strain energy of the crystal is invariant against rigid body rotations, the force on each atom is zero and the stresses on the crystal vanish. The atomistic relaxation simulation which relaxed the atoms into the lowest local energy configuration and the periodic boundary conditions on the supercells in conjunction with the Parrinello-Rahman Lagrangian ensures that these conditions are satisfied.
The elastic coefficients can then be expressed in their more familiar form by pairing the indices by the equivalences , , , , and .
The bulk modulus, (for a cubic crystal) is then
Table shows the independent elastic coefficients for the various supercells in the two simulations. The entire elastic matrix is calculated in all cases. It is found that , and . All other coefficients are almost zero. The bulk modulus is higher than that of experiment ( =0.61 eV/Å ) because the parameters of the potential were fitted for C=0, although the elastic coefficients are sensitive to the value of C. It can be seen, however, that the change in bulk modulus seems to be almost independent of the type of defect under consideration. The elastic coefficients are also calculated for the smaller supercells, which allows the change in elastic coefficients with respect to oncentration to be calculated. The coefficients are seen to rise with concentration, (Figure 7.13),apart from C which decreases very slightly. This shows the characteristic defect stiffening in silicon[74, 75].
Figure 7.13: Change in elastic coefficients with respect to concentration are shown for the hexagonal interstitial. It is found that the change in elastic coefficients with respect to concentration, , is fairly constant over the range considered. Note the change in scale for C .