Introduction

The aim of this experiment is to get familiarised with DLVVisualize and the GULP by modelling the sodium chloride like structure MgO and computing some of its properties. The accuracy of those properties are questioned and investigated by altering the grid shrinkage of the computed system, to find the ideal balance between computational cost and significant results. A comparison between the quasi-harmonic approximation and a molecular dynamics method is made when computing the thermal expansion coefficient for the crystal and their limitations are explored. While the computations are simple from the users point of view, obtaining and demonstrating the underlying theoretical knowledge is key in this report. Details of the methodologies and theories involved are further explained and discussed in each section of the experiment.

Quasi-harmonic approximation

The QHA is an adaption of the harmonic approximation which is based on the harmonic oscillator, but adds a temperature dependent entropic term to the concept. This allows the variation of the equilibrium distance of the crystal's vibrations under a change in temperature, which in the simple harmonic approximation stayed constant. Physically this allows the system to model and compute the thermal expansion coefficient, i.e. the change in volume of the crystal with change in temperature. The Helmholtz free energy under the QHA is given by $F (T, V) = E_{0} + \frac{1}{2} \sum_{k, j} ℏ ω_{k, j} + k_{B} T \sum_{k, j} \ln [1 - \exp (- \frac{ℏ ω_{k, j}}{k_{B} T})]$ , where the first two terms relate to the internal energy and the vibrational zero point energy (dependent on the frequency of the j^th phonon mode at k, $ω_{k, j}$ ), whereas the third term corresponds to the vibrational degrees of freedom and is entropically driven and dependent on temperature.

Molecular Dynamics

The Molecular Dynamics method is based on a classical approximation which assumes the atoms as elastic hard spheres. They are assigned an initial configuration in space and a random velocity vector which overall averages out to a corresponding temperature. An algorithm is then looped over the system and the interactions ("accelerations") update the velocities and consequently the spacial arrangements of the atoms. This is done for a particular number of time steps at a time step of 1 femto second in order to capture vibrations of all frequencies.

Initial Calculation of MgO

The internal energy of one unit cell of MgO was calculated in this exercise. This was done by assuming a fully ionic lattice by assigning a +2 and -2 charge to Mg and O respectively which is justifiable as this behaviour was observed using Hartree-Fock and density-functional quantum-mechanical calculations. ^[1] Their interaction consists of a long range Coulombic potential and a short range interaction described by the Buckingham repulsive potential of $ϕ_{i j} (r) = A_{i j} e x p (\frac{- r_{i j}}{ρ_{i j}}) - \frac{C_{i j}}{r_{i j}^{6}}$ . This is called the pair potential approximation and assumes that the potential only depends on the internuclear distances (using the two-body potential term) rather than having an additional angular dependence (more body terms are neglected). This interaction energy is added to the the self-interaction energy of the atoms, which arrises from the atoms being brought together from infinite distance of the gaseous state. In the calculation for MgO had a value of -41.07531759 eV per primitive unit cell.

Calculating the Phonon Modes of MgO

Computing the Phonon dispersion curves

The phonon dispersion curve for MgO was calculated over the W-L-G-W-X-K path in the k-space, which correspond to the lattice vibrations of the crystals. As it is very difficult to plot all vibrations of all the k vectors in a three dimensional space, the previously mentioned path was chosen (defined by symmetry points) and 50 points along it were plotted on the 2D graph seen below. 6 branches are generated which was predicted by the fact that there are two atoms in the first Brillouin zone and they have each three degrees of freedom. ^[2] The three branches with zero frequency at G are acoustic branches, which originate from the coherent movement of the atoms within a primitive unit cell. The zero frequency is reached as the wavelength of the vibration goes to infinite (i.e. when all the atoms in the crystal move in the same direction). The other three branches are optical branches and originate from the atoms within a cell moving towards each other, at lower wavelengths and thus normally at higher frequencies than the acoustic branches.


Figure 1. Phonon dispersion for the MgO	Figure 2. MgO DOS for a 1x1x1 k-space grid

The Phonon Density of States

The Density of States was computed for the grid sizes 1x1x1, 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16 and 32x32x32. As the size of the grid increased, so did the available states for more frequencies as displayed by the appearance and eventual broadening of additional peaks. This is due to the computing of additional k-points which leads to a wider range of available frequencies. Judging by eye the change between the grid size of 16x16x16 and 32x32x32 only differ slightly by an additional smoothing of the curves. So it can be said that the a sufficiently reasonable approximation can be made using a grid size in-between, i.e. a shrinking factor of 24. At these higher grid sizes the individual vibrations of the k-points are no longer distinguishable and a continuous curve is observed instead.

Figure 3: 2x2x2 grid.	Figure 4: 3x3x3 grid.	Figure 5: 4x4x4 grid.
Figure 6: 8x8x8 grid.	Figure 7: 16x16x16 grid.	Figure 8: 32x32x32 grid.

The DOS of the grid with shrinking factor 1 corresponds to the k-point of L because there are four bands available with the two peaks at lower frequencies being degenerate (both have about double the size in comparison to the higher frequency bands). Additionally the frequency at which the DOS occur are the same as in the Phonon Dispersion diagram for L (just below and above 350 cm^-1 and two bands around 700 cm^-1. The connection between those two diagrams is that they represent the available energy states (in terms of frequency) in different energy intervals. For example when interchanging the axis of the Phonon Dispersion diagram and looking at a frequency range (i.e. around 400 cm^-1) shows that a lot of energy states are available for a wide range of k-points as indicated by the optical band in that region. This is also the reason there is a maximum in that frequency region when looking at the DOS diagram at higher shrinkage.

Further Speculations

The optimal grid size for MgO was found to be the 24x24x24. The same grid size would be applicable to CaO as well as its crystal structure, the relative charges of the ions and the ionic nature of the bonding are very similar. Zeolite however, has a much bigger primitive unit cell in real space with much larger cell parameters than MgO. Due to the inverse relationship between the cell parameters of real and reciprocal space, the lattice parameter in latter is much smaller compared to MgO. Thus a smaller grid size is sufficient to compute a DOS of similar accuracy. One would expect the opposite behaviour of a metal like lithium as it has a comparable small lattice parameter in real space as its primitive unit cell only consists of one atom. However, due to the nature of the delocalised electrons in the metal which "dampen" the vibrations of the crystal structure and thus minimise the repulsive forces. This leads to less fluctuation in the vibrational energy levels, so that a smaller grid size can be used in comparison to MgO as less k-points are necessary to get an adequate accuracy.

Computing the Free Energy - The Harmonic Approximation

Selecting a suitable k-space grid

Figure 9. Helmholtz free energy per unit cell vs grid size

As it is impossible to compute the sum over the vibrational modes of the infinite crystal a finite crystal is used instead. This section investigates the trade off between computation time and accuracy of the computed Helmholtz free energy. The free energies have been calculated for the same grid sizes as in the previous section and the results are displayed in the table below. In addition another column was added to show the deviation from the energy at infinite grid size which was extrapolated by graphing the data and found to be equivalent (to 7 decimal places) to the grid size of 32x32x32. From the table can be seen that the 2x2x2 grid is already accurate to 1 and 0.5 meV, while the 3x3x3 grid is accurate to 0.1 meV. As the grid with shrinkage of 32 had slightly larger computational time it was decided that the future calculations will take place using a grid of shrinkage 24.

Grid size	Free energy per primitive unit cell /eV	Deviation from the energy at an infinite grid size /meV
1x1x1	-40.9303	3.818
2x2x2	-40.9266	0.126
3x3x3	-40.9264	0.051
4x4x4	-40.9265	0.033
8x8x8	-40.9265	0.005
16x16x16	-40.9265	0.001
32x32x32	-40.9265	0.000

Thermal Expansion of MgO and Speculations

In this exercise the variation of the free energy upon a change in temperature of the system of size 24x24x24 was computed. A temperature range of 0 to 1000 K was chosen and the results for the change in free energy and lattice constant is displayed in the graphs below.


Figure 10. Energy vs Temperature	Figure 11. Lattice Constant vs Temperature

The Helmholtz free energy in the quasi-harmonic approximation is given by $F (T, V) = E_{0} + \frac{1}{2} \sum_{k, j} ℏ ω_{k, j} + k_{B} T \sum_{k, j} \ln [1 - \exp (- \frac{ℏ ω_{k, j}}{k_{B} T})]$ , where the first two terms relate to the internal energy and the vibrational zero point energy, whereas the third term corresponds to the vibrational degrees of freedom and is entropically driven and dependent on temperature. At higher temperature this term dominates and is also negative (as the natural logarithm of a number between zero and one will be negative). Thus it was expected that the free energy decreases with increasing temperature, which was observed. Physically this is explained by the availability of new vibrational ground states with increasing temperature and thus a higher multiplicity.

Contrary to the decrease of free energy, the lattice parameter for the cell increases with increasing temperature, i.e. the cell is expanding. By graphing the lattice volume which was obtained from the output files against the temperature, one can compute the thermal expansion coefficient which is given by $α_{V} = \frac{1}{V_{0}} {(\frac{\partial V}{\partial T})}_{P} \approx \frac{Δ V}{V_{0} Δ T}$ . The points between 0 and 300 K were neglected for the linear fit as they deviate due to their low contribution (at low temperatures) of the entropic term in the Helmholtz free energy equation mentioned above. It can be seen that they are dominated by the first two terms and rather behave linearly (at Temperatures below 100 K) and then have equal contributions of both parts of the equation (between 100 K and 300 K). A V₀ was extrapolated from the y-intercept of linear fit and the thermal expansion coefficient was calculated as α_V = 2.99 x 10^-5 K^-1 at 0 K or α_V = 2.96 x 10^-5 K^-1 at 300 K (obtained by dividing through the volume at 300 K rather than extrapolating V₀). The latter value is relatively close to literature ^[3] which obtained a thermal expansion coefficient of α_V = 3.12 x 10^-5 K^-1 at 300 K. Physically, increasing the temperature in a crystal increases the average displacement around their equilibrium position of vibrating atoms in a crystal, thus leading to an expansion of the crystal.

The main approximations in the calculation were firstly that the quasi-harmonic approximation holds true at increasing temperatures, which would imply that bonds don't break once the melting point is reached due to the parabolic shape of the potential energy curve. Around this temperature the phonon modes would inaccurately represent the actual vibrations as they would just continue to imitate the behaviour of a solid even though in reality the MgO would have turned liquid. In fact, if an exactly harmonic potential is used, the bond length would not increase at all with increasing temperature, due to the parabolic shape of the potential energy curve (the average value of the vibration remains always around r₀. However, using quasi-harmonic approximation adds the temperature dependent entropic term as discussed previously, which accounts for the thermal expansion. Secondly, it is assumed that the ions behave fully ionic as hard spheres, which disregards any covalent effects and any multi body interactions as described previously are also neglected, essentially making the atoms an ensemble of independent quasi-harmonic oscillators.

Molecular Dynamics and Speculations

Figure 13. Comparing the Volume vs Temperature of QHA and MD

Another method to calculate the thermal expansion coefficient using the classical interpretation of assigning velocity vectors to the atoms in the crystal structure and evaluating desired properties from the average behaviour of all atoms. In order to get an adequate sample set and remove any ideal behaviour caused by a very small unit cell number, a cell containing 32 MgO units was used. The obtained final volumes from the output files are plotted in the same volume vs temperature graph from before to compare the quasi-harmonic approximation to the molecular dynamics calculations. The MD method seems almost completely linear, has a steeper gradient and at higher temperatures almost converges with the QHA calculation. Thus the calculated thermal expansion coefficient is higher with α_V = 3.09 x 10^-5 K^-1 at 0 K which is less than 1% away from the previously mentioned literature value.

Both methods have their advantages and disadvantages, and are more significant at different temperature ranges. I.e. the classical approximation of the MD is justified at higher temperatures as vibrations are no longer harmonic and the QHA breaks down as discussed previously. At lower temperature however, quantum mechanic effects dominate and the QHA models them with good accuracy, as entropic contributions are minimal. So when especially looking at the temperature range of phase transitions MD is invaluable, but QHA is sufficient for the temperature range of this experiment (especially since it comes at a lower computational cost as the modelling of less cells are required).

Conclusion

The different phonon DOS were investigated for different grid sizes and it showed that the diagram became less noisy at larger grids. For the calculations of the free energy and lattice constants a grid size of 24x24x24 was chosen to be ideal as it showed a good balance between computational time and result accuracy. The consequent calculations as a function of temperature showed expected behaviour, i.e. decreasing free energies with increasing temperature due to the increasing negative entropic contribution and increasing lattice constants with increasing temperature due to the increase in average displacement from the equilibrium position from vibrations. Consequently the QHA was compared to the MD method for computing the thermal expansion coefficient. The QHA proved to show better results at low temperatures due to minuscule entropic contributions and dominating quantum mechanical effects, while the MD method was better applicable at higher temperatures where the classical approximation is more adequate and gives better results. Latter method in particular deviated from the literature value by only 1%, even thought the QHA was sufficient for this experiment. Further research could potentially include variations in pressure on the system, as well as the implementation of various defects and an investigation on how they influence those properties.

References

↑ Timothy S. Bush. "Self-consistent interatomic potentials for the simulation of binary and ternary oxides." Journal of Materials Chemistry 4.6 (1994): 831-837.
↑ Peckham, G. "The phonon dispersion relation for magnesium oxide." Proceedings of the Physical Society 90.3 (1967): 657.
↑ Anderson, Orson L., and Keshan Zou. "Thermodynamic functions and properties of MgO at high compression and high temperature." Journal of Physical and Chemical Reference Data 19.1 (1990): 69-83.

[1] Timothy S. Bush. "Self-consistent interatomic potentials for the simulation of binary and ternary oxides." Journal of Materials Chemistry 4.6 (1994): 831-837.

[2] Peckham, G. "The phonon dispersion relation for magnesium oxide." Proceedings of the Physical Society 90.3 (1967): 657.

[3] Anderson, Orson L., and Keshan Zou. "Thermodynamic functions and properties of MgO at high compression and high temperature." Journal of Physical and Chemical Reference Data 19.1 (1990): 69-83.

[1]

[2]

[3]