Rep:Mod:ia2514 MgO

Abstract

In this experiment, the thermal expansion of a magnesium oxide crystal was studied within the quasi-harmonic approximation and employing molecular dynamics using the molecular modelling code GULP and DLVisualize and the results obtained from the two methods were compared. The internal energy, phonon dispersion and the density of states for MgO were computed. The Helmholtz free energy was calculated within the quasi-harmonic approximation for various k-space grid sizes, leading to the selection of a 32x32x32 grid size for the consequent structure optimisation of MgO at various temperatures in the 0-1000K range. Molecular dynamics simulations were run for the same temperatures as in the quasi-harmonic approximation and the data obtained from the two methods was used to plot the cell volume as a function of temperature and to compute the thermal expansion coefficient, the quasi-harmonic approximation results being in better agreement with literature values. The quasi-harmonic approximation was therefore found to be more accurate than the molecular dynamics method in the temperature range considered for the simulations. The benefits and disadvantages of the two methods were discussed.

Introduction

Magnesium oxide has a rocksalt structure with a FCC lattice and an experimentally measured lattice constant of 4.217 $\AA$. By assuming periodic boundary conditions (and thus no crystal defects), the properties of the crystal can be studied computationally by considering only one (or a finite number) of unit cells. One of the behaviours that can be simulated in this way is the thermal expansion of the crystal, allowing for the study of the variation of the cell size with temperature and the calculation of the thermal expansion coefficient:

$\alpha =\frac{1}{V}\frac{\partial V}{\partial T}$ (1)

The coefficient is dependent on temperature and for small temperature changes, it can be approximated by:

$\alpha =\frac{1}{V_0}\frac{\mathrm{\Delta }V}{\mathrm{\Delta }T}$ (2)

One of the approximations used to model particle vibrations is the harmonic approximation, which implies that near the equilibrium bond distance particles undergo harmonic motion, a result based on the Taylor expansion of the potential energy. While it is a good starting point for calculations of lattice dynamics, an exactly harmonic potential does not account for the change in equilibrium bond distance with increasing temperature and therefore cannot explain thermal expansion. The origin of the thermal expansion can be understood by considering an anharmonic potential where the increase in temperature leads to an increase in the equilibrium bond length since the potential curve is no longer a symmetric parabola. The quasi-harmonic approximation extends the harmonic approximation by assuming that phonon frequencies depend on volume and by treating the equilibrium volume as a free parameter which is obtained at each temperature by minimising the free energy. The change in volume with temperature therefore renders the quasi-harmonic approximation suitable for the study of thermal expansion.

Another method which can be used to predict the thermal expansion of a solid is the molecular dynamics method, which is employed for studying the motion of particles by treating them classically. This method is implemented by updating the positions and velocities of each atom at each time step by computing the forces that act on it due to interactions with the other atoms.

Regardless of the method used, performing energy and vibration calculations for a crystal requires a compromise between accuracy and time efficiency. Although the calculations would theoretically require considering an infinite lattice and summing over all the possible k-points in the reciprocal lattice, approximations can be made in order to reduce CPU time. By assuming periodic boundary conditions, only a finite number of unit cells can be studied and an optimal k-space grid size can be used for sampling.

Methodology

All of the calculations on the MgO crystals were performed using the molecular modelling code GULP with DLVisualize as the graphical user interface in the RedHat Linux environment. The ionic interactions were modelled using a force field with a Buckingham repulsive potential between the Mg ion with a +2e charge and the O ion with a -2e charge. The obtained data was further processed using Excel and Python.

The internal energy of MgO per primitive unit cell was computed in a single point calculation. The phonon dispersion curves were computed at 50 points in the k-space along the W-L-G-W-X-K path and the calculation also allowed for the visualisation of the vibrations as animations. A calculation for the phonon Density of States was run for a 1x1x1 k-space grid and the resulting DOS graph was compared with the phonon dispersion curves in order to identify the k-point at which the calculation was computed. The density of states was then computed for 1x1x1, 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 32x32x32 and 64x64x64 grids in order to decide on an optimal grid size for the coming calculations.

The Helmholtz free energy of the MgO crystal was calculated at 300K for 1x1x1, 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 32x32x32 and 64x64x64 k-space grids. The structure of MgO was optimised by minimising free energy for temperatures in the range 0-1000K in steps of 100K and by sampling over a 32x32x32 grid. The calculations were run both by starting from the same initial geometry for all temperatures or by starting from the geometry as obtained from the preceding optimisation (with increasing temperature) and, as expected in the case of an optimisation, the results were the same regardless of the starting geometry. The free energy was extracted from the log files and plotted as a function of temperature.

Equation (3) was used in order to compute the lattice constant for the conventional unit cell (for which a_c=b_c=c_c) from the primitive unit cell parameters (a_p=b_p=c_p) and the results were plotted against temperature.

$a_p=\frac{1}{2}\sqrt{a_c^2+b_c^2}\Rightarrow a_c=\sqrt{2}{a}_p$ (3)

The cell volume as a function of temperature was fitted to two different polynomial functions in two separate ranges and the result was used to calculate the thermal expansion coefficient using both equations (1) and (2) and the accuracy of the approximation given in equation (2) was discussed.

Molecular dynamics simulations were run for a supercell containing 32 MgO units at temperatures in the range 0-1000K in steps of 100K with a time step of 1.0 femtoseconds, 500 equilibration and production steps and 5 sampling and trajectory steps. The conventional cell volume at different temperatures was calculated from the average of the last 10 supercell volumes divided by 32. The data was compared to the one obtained from the quasi-harmonic approximation. The data points were fitted to a line and the thermal expansion coefficient calculated using both equations (1) and (2) and the results were compared amongst themselves to investigate the accuracy of method (2). The values obtained from both the quasi-harmonic approximation and molecular dynamics were compared to literature values.

Results and discussion

MgO internal energy, phonon dispersion and DOS

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

The internal energy of the MgO crystal as extracted from the single point calculation was, as expected, -41.07531759 eV per unit cell (-3963.1403 kJ/mole unit cells).

The MgO phonon dispersion (the frequency of the phonons as sampled over 50 point along the W-L-G-W-X-K path) is shown in figure 1. Since it would be difficult to plot the energies over the three dimensional k-space grid, the phonons are computed over a finite number of k-points on a given path that connects symmetry points within the lattice (in this case, W-L-G-W-X-K).

By comparing the phonon DOS graph obtained by sampling a 1x1x1 k-space grid (figure 2) with the MgO phonon dispersion (figure 1) it can be concluded that the DOS was computed for the L point since the frequencies of the peaks in the DOS match those in the dispersion curves for L and the double intensity of the DOS peaks between 200 and 400 cm^-1 is in agreement with the double degeneracy of the bands in the phonon dispersion graph at the same frequencies. This observation is supported by the coordinates of the k-point as given in the DOS log file: (0.5; 0.5; 0.5), coordinates which correspond to the same frequencies in both the dispersion curves calculation and DOS log files.

The MgO phonon DOS was also computed for 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 32x32x32 and 64x64x64 k-space grids and the results are shown in figure 3. It can be seen that as the grid size increases, the number of peaks increases (due to more and more vibrations being computed) and the graph eventually becomes smooth.

Figure 3. MgO Phonon DOS for increasing grid sizes from 1x1x1 to 64x64x64	Figure 4. MgO Phonon DOS for a 32x32x32 grid size (optimum size)

32x32x32 was chosen as the minimum grid size for which the DOS is a good representation of the infinite lattice (figure 4). For smaller grids the DOS plot is noisy and individual vibrations at the few k-points sampled can be distinguished, rather than observing a continuous curve, whereas for the larger grid size sampled there are no significant differences and therefore 32x32x32 was selected as the optimum size to yield an accurate representation of the real DOS, while being efficient in terms of CPU time.

This grid size is expected to give equally accurate results in the case of a similar oxide to MgO (both in terms of size and structure), such as CaO. For a Zeolite, such as Faujasite, which has a significantly larger unit cell than MgO, the reciprocal lattice is smaller than that of MgO (since its size is inversely proportional to the real space lattice) and a smaller k-space grid size might be enough to yield results of similar accuracy. The opposite is expected for metals, such as lithium, given that their lattice size is smaller than that of MgO and their reciprocal lattice larger; they would therefore require a sampling of a larger grid than 32x32x32 for a satisfactory outcome.

The density of states at a given energy E represents the number of energy levels (both occupied and unoccupied) between E and E+dE per energy interval dE. The phonon DOS can be thus easily computed from the phonon dispersion curves by considering the number of points on the graph at each energy value (within the dE interval) in the dispersion curves plot. In other words, the flatter the dispersion curve is (the closer the slope is to 0), the higher the DOS. This can be seen in the case of the MgO DOS (figure 4), where the highest DOS corresponds to a frequency of around 400 cm^-1, and not only do several bands cross this value in the dispersion curves plot, but there is also a relatively flat band around the $\mathrm{\Gamma }$ point at approximately 400 cm^-1 (figure 1).

The Helmholtz free energy of MgO within the quasi-harmonic approximation

Figure 5. MgO Helmholtz free energy per unit cell as a function of grid size

The variation of the Helmholtz free energy of MgO at 300K with increasing k-space grid size is shown in figure 5. For the small grid sizes, the free energy increases abruptly, then starts decreasing slowly and reaches a plateau at -40.9265 eV. The obtained values along with the differences in energy from this value at which the free energy tends with increasing grid size are given in table 1. It can therefore be seen that from the sampled grid sizes, the smallest one that can be used for calculations accurate to 1 meV or 0.5 meV is the 2x2x2 grid and for calculations accurate to 0.1 meV the 3x3x3 grid is the smallest one that should be used. The 32x32x32 grid size, just as before, is the optimal one for accuracy and time efficiency.

Table 1. MgO Helmholtz free energy for various grid sizes

Grid size	Helmholtz free energy per unit cell /eV	Deviation from infinite grid size energy /meV
1x1x1	-40.9303	3.8180
2x2x2	-40.9266	0.1260
3x3x3	-40.9264	0.0510
4x4x4	-40.9265	0.0330
8x8x8	-40.9265	0.0050
16x16x16	-40.9265	0.0010
32x32x32	-40.9265	0.0000
64x64x64	-40.9265	0.0000

Free energy, lattice constant and thermal expansion coefficient as functions of temperature within the quasi-harmonic approximation

Figure 6. MgO Free energy as a function of temperature	Figure 7. MgO lattice constant as a function of temperature (as calculated within the quasi-harmonic approximation)

The variation of the Helmholtz free energy of MgO as a function of temperature is shown in figure 6. As expected, given that A=U-TS, the free energy decreases with increasing temperature. The entropy increases with temperature, and although the internal energy also varies with temperature, the -TS term dominates resulting in an overall decrease in free energy.

As expected, the quasi-harmonic approximation does not fail to predict the thermal expansion of MgO since, as seen in figure 7, the lattice constant is increasing with temperature. While the quasi-harmonic approximation is likely to accurately predict the thermal expansion at low temperatures, at higher temperatures (closer to the melting point), phonon-phonon anharmonic interactions become significant and the phonon modes as computed within the quasi-harmonic approximation fail to accurately represent the motions of the ions. However, in the present case, the considered temperature range (0-1000K) lies significantly below the melting point of MgO (3100K), and therefore the anharmonic interactions can be considered negligible.

Figure 8. MgO primitive cell volume as a function of temperature (as calculated within the quasi-harmonic approximation)	Figure 9. MgO thermal expansion coefficient as a function of temperature (as calculated within the quasi-harmonic approximation)

The plot of MgO cell volume as a function of temperature was fitted to a third order polynomial (-6.27x10^-10xT³+1.07x10^-6xT²-8.80x10^-5xT+1.88x10) in the 0-400K range and to a line (5.47x10^-4xT+1.87x10) in the 400-1000K range (figure 8). The two obtained functions were used to calculate the partial derivative of volume with respect to temperature and thus compute the thermal expansion coefficient using equation (1). The thermal expansion coefficient was also calculated by assuming linearity within a 100K temperature range and by employing equation (2). The results yielded by the two methods are presented in figure 9. The largest deviation in the thermal expansion coefficient computed using equation (2) from the results obtained through equation (1) occurs at 0K and has a value of 0.48x10^-5. Above 300K, the difference between the two values at any temperature does not exceed 13% of the value of the thermal expansion coefficient calculated using the derivative of the fitted data.

The comparison of two of the thermal expansion coefficient values obtained within the quasi-harmonic approximation with experimental values (table 2) reveals that the computed values are of the right order of magnitude and that the values calculated by assuming linearity over 100K temperature intervals are closer to literature values, yet still approximately 30% lower than expected. This observation does however suggest that using equation (2) in the temperature range 400-1000K gives more accurate results than those obtained by attempting to fit the data to a straight line and compute the thermal expansion coefficient using its slope. The latter approach leads to a decrease in the thermal expansion coefficient with increasing temperature, a trend which goes against the one observed in the experimental data.

Table 2. MgO thermal expansion coefficient data comparison

Temperature /K	Thermal expansion coefficient /10-5K-1
	Calculated by equation (1)	Calculated by equation (2)	Literature value
900	2.85	3.05	4.38
1000	2.84	3.15	4.47

Quasi-harmonic approximation versus molecular dynamics

The average values for the unit cell volume were calculated from the molecular dynamics calculations performed for temperatures in the 0-1000K range and their variation with temperature is shown in figure 10. As expected, the molecular dynamics method also predicts the thermal expansion of MgO since the cell volume increases with temperature. The MD data points follow a more linear trend than the QH data and the two approach convergence at high temperatures, the difference between the values calculated with the two methods decreasing with temperature.

Figure 10. MgO primitive cell volume as a function of temperature (quasi-harmonic approximation and molecular dynamics)

The thermal expansion coefficient was also calculated from the molecular dynamics data, using both equations (1) and (2) and the values were plotted in figure 11. The MD values follow a more linear trend than the QH values and their dependence on temperature is less significant, a trend which is in disagreement with the increase in the experimental values of the thermal expansion coefficient with increasing temperatures. The MD values obtained at 900K and 1000K are further away from the literature values than the QH values (table 2).

Figure 11. MgO thermal expansion coefficient as a function of temperature (quasi-harmonic approximation and molecular dynamics)

While both the quasi-harmonic approximation and the molecular dynamics method predict the thermal expansion of MgO and yield numerical results of the right orders of magnitude (even for a k-space grid size of 32x32x32 and a supercell containing 32 MgO units, respectively), over certain temperature ranges they fail to produce accurate values or the expected trends. At very low temperatures, molecular dynamics calculations are expected to fail since quantum effects are significant and MD treats particles classically. In contrast, at high temperatures, where the vibrations are no longer harmonic, the quasi-harmonic approximation is unsuitable since it does not account for the anharmonic phonon-phonon interactions. MD is therefore better suited for studying behaviour at phase transitions and therefore at higher temperatures where atoms behave more according to classical mechanics. On a large temperature interval, ranging from 0K to melting point, plots of the cell volume against temperature as obtained from the quasi-harmonic approximation and molecular dynamics are therefore predicted to start from different points at low temperatures, converge in the middle region and then diverge again towards high temperatures, where the quasi-harmonic approximation underestimates the actual cell volume which is accurately predicted by molecular dynamics. However, since the temperature range studied in the present simulations is considerably lower than the melting point of MgO, overall the quasi-harmonic approximation produces more accurate results that the molecular dynamics method, which is expected to return better outcomes at temperatures higher than 1000K.

An additional reason for the poor performance of molecular dynamics in the studied temperature range is the limitation due to the small size of the system under investigation (a supercell containing 32 MgO units). As concluded from the calculations within the quasi-harmonic approximation, a 32x32x32 k-space grid size is optimum, which is the equivalent of sampling a 32x32x32 MgO cell in molecular dynamics. As shown in table 1, the 16x16x16 grid size is expected to yield accurate free energy calculations within 0.001 meV, suggesting that a 16x16x16 supercell could produce reliable outcomes in molecular dynamics. However, performing molecular dynamics simulations on this supercell size would still be too computationally demanding.

Additionally to the aforementioned intrinsic shortcomings of the two methods, another approximation which limits the accuracy of the results in comparison with experimental data is the assumption that the MgO lattice is infinite and without defects.

Conclusions

The thermal expansion of magnesium oxide was investigated by applying two different methods: the quasi-harmonic approximation and molecular dynamics. Calculations within the quasi-harmonic approximation were conducted by sampling over a reciprocal space grid of size 32x32x32, which was determined to be the optimum size for accurate, yet time efficient results when computing the MgO phonon DOS and the free energy for various grid sizes. A supercell containing 32 MgO units was used for the molecular dynamics simulation as a compromise between reproducing as many vibrations as possible within a finite cell and computational time. Within the 0-1000K temperature range, the two plots of variation of cell volume with temperature as computed from the two different methods start converging towards higher temperatures. The thermal expansion coefficient as calculated within the quasi-harmonic approximation was found to be in better agreement with experimental data than the one computed through molecular dynamics, yet still about 30% lower than the expected values. Overall, within the given constraints (k-space grid size, MD supercell size), the quasi-harmonic approximation was deemed more suitable for the study of the thermal expansion of MgO in the 0-1000K range, since at low temperatures, quantum effects become significant and the molecular dynamics method does not account for them and the considered temperatures are low enough for anharmonic interactions (which are ignored in the quasi-harmonic approximation) to be negligible. Nevertheless, at higher temperatures, closer to the melting point, molecular dynamics is expected to simulate the behaviour of magnesium oxide more accurately, a prediction which could be investigated as a future extension to the present experiment.