Modeling MgO

Introduction

The thermal expansion coefficient, α, is an important property to know and so it would be very useful if computers could calculate it using models. Magnesium Oxide (MgO) was the system that the calculations were run on and two different approaches where taken to calculate the coefficient. They both calculate the cell volume of the system and by plotting the volume vs temperature, the coefficient could be calculated by using this equationː

$α = \frac{1}{V_{0}} {(\frac{\partial V}{\partial T})}_{P}$

where T is the temperature, V is the volume and V₀ is the inital volume.

Quasi-Harmonic Approximation

A phonon is a quantum of energy associated with the vibrations of a crystal. Models, such as the quasi-harmonic approximation, can be based on phonons to describe the thermal expansion of a system. The harmonic phonon model is based on approximating the Leonard-Jones potential into a parabola with the minimum at the equilibrium bond length. However, it cannot explain thermal expansion because in that approximation, the equilibrium bond length does not change with a change in temperature. The quasi-harmonic model resolves this by shifting the parabola which allows the bond lengths to change according to the temperature.

Molecular Dynamics Model

This model does not use the harmonic approximation; instead it simulates the motions of the atoms in the system and averages over all the vibrational modes of the crystal. It does this by computing the forces on the atoms, accelerating the atoms using F=ma, moving the atoms by updating the velocities and positions and repeating until the properties settle. This can only be done in a supercell which is 2x2x2 because it needs to account for the movements of the atoms not being periodic.

Results

Phonon Modes of MgO

Phonon dispersion curves show the frequencies along a path in k-space with gamma as the origin. The phonon density of states (DOS) is the average of all the k-points giving the number of vibrational frequencies which can be viewed on the phonon dispersion curve (figure 1). A shrinking factor of 1 (1x1x1 grid) will yield a DOS curve (figure 2) and by using figure 1, it can be seen that L is the k-point for the 1x1x1 grid because L has four bands, each frequency correlating to the frequencies of the four peaks in figure 2. A 3D system will produce six vibrational bands, however the L point has only four bands because two bands are degenerate. This can be seen on the DOS curve because two of the peaks near 400 cm^-1 are roughly double the size of the other two peaks.

Figure 2: 1x1x1 grid.	Figure 3: 2x2x2 grid.
Figure 4: 3x3x3 grid.	Figure 5: 4x4x4 grid.
Figure 6: 8x8x8 grid.	Figure 7: 16x16x16 grid.
Figure 8: 32x32x32 grid.	Figure 9: 64x64x64 grid.

As the shrinking factor gets larger, the number of peaks in the DOS curves increases because there are more k-points to compute. Eventually, there are so many peaks that they all converge to create a curve as seen on the 16x16x16 grid, figure 7. The peak at approx. 400 cm^-1 is greater than the other peaks because there are a great number of degenerate bands at that frequency. The larger the shrinking factor, the longer the time it takes to compute the calculations and so the minimum shrinking factor needed for a reasonable approximation of the DOS would be 32 because it gives the same maximum as that of a 64x64x64 grid (396 cm^-1) but is faster because it only takes 14 seconds to run, whereas a shrinking factor of 64 takes 681 seconds to complete the calculations.

Now that an ideal grid size has been found for MgO, assumptions can be made to other systems about their grid sizes. For CaO, calcium has the same number of valence electrons and charge as magnesium and so it can be assumed that the same grid size can be used for that system. A zeolite would have a bigger primitive cell than MgO in real space and so will have a smaller lattice parameter in reciprocal space due to the inverse relationship of the lattice parameter in real space and the lattice parameter in reciprocal space; therefore a smaller grid size is needed with fewer k-points to obtain an accurate DOS curve. A metal system, such as lithium, would have a smaller primitive cell in real space, however it would not have a larger grid size because the free electrons between the ions cause less vibrations due to repulsive forces and so the atoms have vibrational energy levels that fluctuate less. This means that the DOS curve does not need to have as many k-points to be as accurate; the grid size will be smaller than that of MgO.

Quasi-Harmonic Approximation

The Helmholtz free energy for MgO can be computed using the quasi-harmonic approximation to an accurate degree from the grid size that was optimized previously. The free energy values converge to a value that the 32x32x32 and 64x64x64 grid size both share, the values first increase as the grid size increases and then decrease to a plateau. Assuming that the 64x64x64 grid is the real value, it can be seen from table 1 that the 1x1x1 grid is not accurate to 1 meV but 2x2x2 is accurate to 0.5 meV, 4x4x4 is accurate to 0.1 meV and 8x8x8 is accurate to 0.01 meV.

Table 1: Free Energies of MgO for different grid sizes
Grid size	Helmholtz Free Energy (eV)	ΔA (mEV)
1x1x1	-40.930301	3.82
2x2x2	-40.926609	0.126
3x3x3	-40.926432	0.0510
4x4x4	-40.926450	0.0330
8x8x8	-40.926478	0.0050
16x16x16	-40.926482	0.0001
32x32x32	-40.926483	0
64x64x64	-40.926483	0

The effect of temperature on the free energy and lattice constant can be seen in figure 10 and 11. The Helmholtz free energy, A, for a lattice in the quasi-harmonic approximation isː

$A = U + U_{0} - T S$

where U is the internal energy, U₀ is the zero-point vibrational energy, T is the temperature and S is the entropy. It can be seen from figure 10 that as the temperature increases, the free energy decreases; this is explained in the equation.

Temperature has an opposite effect on the lattice constant or cell volume as seen in figure 11.

The thermal expansion coefficient for MgO was calculated to be 2.1297×10^-5 K^-1 using the gradient of the cell volume vs temperature curve (figure 14) and assuming the intercept to be the initial volume. This value is similar to the experimental value, 3.706×10^-5 K^-1. ^[1] The reason for the difference in these values come from the many approximations that are taken. It is assumed that the crystal has no impurities or defects, this is usually not true for real crystals. The defects define the properties of the crystals and so could have the potential to change the way that temperature effects the cell volume.

Increasing the temperature to the melting point of MgO would cause the quasi-harmonic approximation to break down because it does not allow the bonds to break therefore at temperatures above the melting point, the phonon modes do not represent the actual motions of the ions.

For a diatomic molecule with a harmonic potential, increasing the temperature would not affect the intermolecular distance because the harmonic potential is symmetrical so the equilibrium bond lengths would stay the same as the energy increases as seen in figure 12. Lennard-Jones ^[2] explained that as temperature increases, the energy increases which causes the intermolecular distance of the ions to increase (figure 13) therefore ionic cohesion is what the thermal expansion originates from.

Molecular Dynamics

Using molecular dynamics to calculate the thermal expansion coefficient proved to be more accurate than using the quasi-harmonic approximation. The cell volume follows the same trend as the lattice constant when the temperature is increased because they are directly proportional to each other. Using the gradient of the curve in figure 14, the thermal expansion coefficient was found to be 3.2134×10^-5 K^-1. Molecular dynamics would also be more accurate at higher temperatures because it would allow the bonds to break unlike the quasi-harmonic approximation. However, molecular dynamics is not able to produce a value at 0 K, it breaks down and is not accurate at that temperature but crystals have cell volumes at 0, which is shown in the quasi-harmonic approximation.

Conclusion

It has been found that when using the quasi-harmonic approximation, a shrinking factor of 32 is best because it saves on time while being as accurate as a calculation based on a shrinking factor of 64. However, different systems will require different grid sizes; this depends on the lattice parameter in reciprocal space and if there are any free electrons to cause shielding effects like in metal systems. After finding the optimal grid size for MgO, it was used to calculate the Helmholtz free energy. The trend for increasing the temperature for the free energy is opposite to the trend for the lattice constant and cell volume. Using the gradient from the cell volume vs temperature curve and the equation for thermal expansion, α could be calculated. The value given by using molecular dynamics was closer to the literature value when in comparison to using the quasi-harmonic approximation to calculate the coefficient. However, at 0 K the quasi-harmonic approximation is more accurate whereas at temperatures above the melting point, the molecular dynamic model is more accurate. This is explained by the many approximations that the quasi-harmonic takes, the main being that it is simply a harmonic potential being shifted to account for the varying bond lengths with changes in temperature.

The next variable to be changed in the system could be the pressure. The pressure and temperature could be changed to calculate how much force the systems requires to be broken or changed. The pressure could also be increased to observe when defects begins to appear, but to calculate defects a more complex model would be required along with software that can stimulate the creation of defects.

References

↑ O. L. Anderson and K. Zou, J. Phys. Chem. Ref. Data, 1990, 19, 1, 69-81. DOI:10.1063/1.555873
↑ J. E. Lennard-Jones "Cohesion" , Proc. Phys. Soc, 1931, 43, 5, 461-482. DOI:10.1088/0959-5309/43/5/301

[1] O. L. Anderson and K. Zou, J. Phys. Chem. Ref. Data, 1990, 19, 1, 69-81. DOI:10.1063/1.555873

[2] J. E. Lennard-Jones "Cohesion" , Proc. Phys. Soc, 1931, 43, 5, 461-482. DOI:10.1088/0959-5309/43/5/301

[1]

[2]