Rep:Mod:TAB9192

A Computational Study of the Thermodynamic Properties of MgO

Abstract

The volumetric coefficient of thermal expansion ($\alpha$) of MgO was computed through the use of both molecular dynamics (MD) simulations and phonon calculations under the quasi-harmonic approximation (QHA). Good agreement between both methods and existing literature was found for temperatures within a range of 100-1000 K. At high temperatures, the MD simulations showed anomalous behaviour in $\alpha$ at temperatures between $3000-3050\mathrm{K}$, matching the experimental melting range at $0\mathrm{P}\mathrm{a}$ closely. The volume at $300\mathrm{K}$ was determined as a function of pressure, and used to calculate the bulk modulus at $0Pa$ (${\mathrm{K}}_0$) of $222.9\pm 0.1\mathrm{G}\mathrm{P}\mathrm{a}$ with a pressure derivative $\mathrm{K}{\text{'}}_0$ of $3.960\pm 0.004$. Again, both MD simulations and QHA-based calculations were performed.

Introduction

Solid MgO adopts the well-known rock-salt structure, which consists of two superimposed but offset face-centred cubic lattices. A detailed theoretical understanding of its thermodynamic properties is paramount to a variety of applications. For example, MgO is a commonly used refractory material, in part due to its chemical stability at high temperatures. In addition, the high-pressure melting of MgO in the earth's lower mantle plays a major role in the seismology of so-called ultra-low velocity zones. However, experimental investigations at such conditions rely on either shockwaves or laser-heated diamond anvil cells, both of which require highly sophisticated equipment and are thus very expensive. A convenient alternative is the use of computational methods. The present work therefore characterises the isothermal compression and thermal expansion of MgO through the use of both MD and QHA-based methods, comparing their accuracy and range of validity.

Methodology

The software used to perform all calculations and simulations was the General Utility Lattice Program (GULP).. DLVisualize was used to provide the user interface.

Lattice Energy of MgO

Both the MD simulations and the phonon calculations under the QHA require a pair potential that determines the inter-particle interactions. Throughout this study, the Coulomb-Buckingham potential was used, taking the following general form:

$U(r)=A\exp (-Br)-\frac{C}{r^6}+\frac{q_1{q}_2}{4\pi {ϵ}_{0}r}$

Where the constants $A$, $B$ and $C$ are chosen to provide the best match to experimental data, and are shown in table 1.

Table 1. Coulomb-Buckingham parameters used.

Pair	A	B	C	Cutoff / $\AA$
O-O	$0.130\times 10^4$	$0.300$	$0.00$	$12$
Mg-O	$0.228\times 10^4$	$0.149$	$27.9$	$12$

The cutoff indicates the value of $r$ above which the short-range Pauli-repulsion and dispersive potential interactions are discarded (the first two terms of the Coulomb-Buckingham potential). Note that the Coulombic interactions were summed using the Ewald method, where interactions are calculated in Fourier-space, due to its significantly faster convergence for periodic systems .

Phonon Calculations under the QHA

Phonons are the collective vibrational modes of a solid. For small distortions from some equilibrium interatomic spacing $r_0$, the potential can be approximated by a parabola, as can be easily verified by Taylor-expansion around $r_0$:

$V(r)\approx V(r_0)+[\frac{\mathrm{d}V}{\mathrm{d}r}{]}_{r=r_0}(r-r_0)+\frac{1}{2}[\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{r}^2}{]}_{r=r_0}(r-r_0{)}^2$

Since $r_0$ is a potential energy minimum, $[\frac{\mathrm{d}V}{\mathrm{d}r}{]}_{r=r_0}=0$, and so, defining $V(r_0)=0$ yields $V(r)\approx \frac{1}{2}[\frac{{\mathrm{d}}^{2}V}{\mathrm{d}{r}^2}{]}_{r=r_0}(r-r_0{)}^2=\frac{1}{2}k(r-r_0{)}^2$.

This is the origin of the 'harmonic' part of the QHA. However, this is clearly not a widely applicable model, as the expactation value of the bond length is equal to $r_0$ for all excited states, and so no thermal expansion could take place.

Under the QHA, the Helmoltz free energy per unit cell consists of contributions from the lattice energy (described in the previous section), the zero-point energy of all phonons in the Brillouin zone, and the energy of occupied excited states (including entropic contributions). The latter two of these are averaged over all posible phononic modes, for example by integration over all possible frequencies, weighted by the (normalized) phonon DOS: $F(V,T)=U+\int \mathrm{d}\omega g(\omega )[\frac{\hslash \omega }{2}+k_{\mathrm{B}}T\log (1-e^{-\frac{\hslash \omega }{k_{\mathrm{B}}T}})]$

The origin of the 'quasi'-part of the QHA is that the volume of the crystal is allowed to vary with a change in conditions in such a way as to minimise the Helmholtz free energy. In general, the force constants of the phonons decrease with increasing volume, and so the vibrational frequencies $\omega$ are lower (a relationship known as Badger's rule ). On one hand, this causes the energy of occupied states (the final term) to become more negative and lowers the zero-point energies, but results in a less negative lattice energy for volumes larger than the equilibrium volume at $0\mathrm{K}$. There are therefore two opposing contributions to the Helmholtz free energy with a change in volume, which will therefore have a minimum at some volume, yielding the equilibrium structure at the given temperature and pressure.

At each volume, this expression is the analytical result derived from the partition function of a harmonic oscillator, averaged over all possible modes. Clearly then, the phonon DOS, $g(\omega )$, contains all of the information required to calculate the Helmholtz free energy at a given temperature. Computationally, for any real system this integral must be approximated by a summation over a finite number of frequencies, as the phonon DOS cannot generally be determined analytically. The phonon DOS is itself summed over a selection of points in the Brillouin zone, which are determined by the shrinking factors. For a 3D lattice, this is a set of three numbers defining the 'resolution' of a mesh of k-points in the Brilluoin zone at which the phonon DOS is calculated. However, a larger number of sampling points comes with greater computational expense, and so, prior to any studies of the thermodynamic properties of a real system, it is important to determine the convergence of the phonon DOS with the shrinking factor, such that a reasonable compromise between accuracy and speed can be reached. This was assessed both visually from plots of the phonon DOS at different shrinking factors, or by monitoring associated properties such as the Helmoltz free energy.

Molecular Dynamics Simulations

Molecular dynamics (MD) simulations are a purely classical approach to modelling the behaviour of systems. The basic principle is that the time evolution of a system obeying Newtonian mechanics is simulated through the use of some iterative algorithm, which integrates the equations of motion at each timestep in order to update the velocities and positions. In this case, the so-called leapfrog algorithm was used, which takes its name from the fact that it offsets the moments in time at which positions and velocities are calculated.

The Coefficient of Thermal Expansion

The thermal expansion of a material can be characterised by its coefficient of thermal expansion: $\alpha =\frac{1}{V}(\frac{\partial V}{\partial T}{)}_P$

The division by the volume of the material is required for $\alpha$ to be an intensive quantity. This generally varies with temperatue, and can thus be computed by differentiation of a polynomial fit to the volume-temperature data (note hat this renders extrapolation meaningless).

The Bulk Modulus

The isothermal bulk modulus (the inverse of the compressibility) of a material can be obtained by fitting the Birch-Murnghan equation-of-state (BMEOS) to volume-pressure data at constant temperature :

$P(V)=\frac{3}{2}{K}_0\left[\right(\frac{V_0}{V}{)}^{\frac{7}{3}}-\left(\frac{V_0}{V}{)}^{\frac{5}{3}}\right][1-\frac{3}{4}(4-K{\text{'}}_0)[(\frac{V_0}{V}{)}^{\frac{2}{3}}-1]]$

Where $V_0$ is the volume at $P=0\mathrm{P}\mathrm{a}$, $B_0=-V\left(\frac{\partial P}{\partial V}\right)$ is the isothermal bulk modulus in the absence of an external pressure and $B{\text{'}}_0=(\frac{\partial B}{\partial P}{)}_{P=0}$ its pressure derivative.

Results and Discussion

The Total Lattice Energy

The total lattice energy per primitive unit cell was determined to be -41.07531759 eV.

The Phonon Dispersion

The phonon dispersion was calculated at 50 points along the W-L-G-X-W-K path in the Brillouin zone, the graph of which is shown in figure 1. Figure 2 contains an analogous (albeit the path in the first Brilluoin zone somewhat differs) graph derived from neutron scattering data, obtained by Sangster et al. . As can be expected from the general rule that the total number of branches is given by $d\times m$, where $d$ is the dimensionality of the system and $m$ the number of atoms in the unit cell, there are six branches, since the primitive unit cell contains one Mg and one O ion. In addition, $d$ of these branches should be acoustic, in other words, have a minimum at $k=\mathrm{\Gamma }=0$. Clearly, this is also the case.

Figure 1: The phonon dispersion calculated at 50 points along the W-L-G-X-W-K path. Three acoustic and three optical branches are clearly visible.

Figure 2: The phonon dispersion as obtained experimentally by Sangster et al. Clearly, the qualitative agreement is excellent.

The Phonon DOS

Figures 3-8 illustrate the convergence of the phonon DOS with increasing mesh resolution (shrinking factor). As can be seen, the change in visual appearance (mainly an increase in the smoothness) of the DOS-curve as a result of a change in the mesh resolution decreases with increasing resolution, to the point where there is no significant visual difference between subsequent DOS-curves at resolutions of 30x30x30 and 40x40x40, and so the convergence has satisfied at least the visual criterion.

Figure 3: Phonon DOS with a mesh resolution of 1x1x1. By comparison of the frequencies present, it can be seen this corresponds to the symmetry point L in the phonon dispersion graph.	Figure 4: Phonon DOS with a mesh resolution of 2x2x2.
Figure 5: Phonon DOS with a mesh resolution of 10x10x10.	Figure 6: Phonon DOS with a mesh resolution of 20x20x20.
Figure 7: Phonon DOS with a mesh resolution of 30x30x30.	Figure 8: Phonon DOS with a mesh resolution of 40x40x40.

Comparison of the phonon DOS for evaluated at only a single point in the Brillouin zone (figure 3) with the phonon dispersion curve (figure 1) shows that this point must correspond to the symmetry point L, as this matches the observed frequencies.

In this context, the phonon DOS is defined as the number of states per unit increase in energy at a given energy i.e. $DOS(E_0)=\frac{\mathrm{d}n}{\mathrm{d}E}$. Using the chain rule, this can be rewritten as $DOS(E_0)=\underset{E(\overrightarrow{k^{\prime }})=E_0}{\int }\mathrm{d}\overrightarrow{k^{\prime }}\frac{\mathrm{d}n}{\mathrm{d}|\overrightarrow{k^{\prime }}|}\frac{\mathrm{d}S(\overrightarrow{k^{\prime }})}{|{\mathrm{\nabla }}_{\overrightarrow{k^{\prime }}}E(\overrightarrow{k^{\prime }})|}$, where the integral along the surface of $E(\overrightarrow{k})=E_0$ arises from the need to consider the contributions of all k-values for which $E(\overrightarrow{k})=E_0$. This expression can of course generally not be computed analytically, but the important thing to note is that the DOS varies inversely with ${\mathrm{\nabla }}_{\overrightarrow{k}}E(\overrightarrow{k})$, i.e. the slope of the phonon dispersion curve. In this specific case, the frequency $\omega$ is used in place of the energy, however, the same relationship holds true, as the zero-point energy of a particular phonon mode is given by $E=\frac{\hslash \omega }{2}⟹E\propto \omega$

Free Energy in the Quasi-Harmonic Approximation

Following the above discussion of the visual convergence of the phonon DOS curve, it should be noted that for the purpose of calculating thermodynamic quantities, it is far more important to know the level of accuracy that can be expected at a given mesh resolution. Thus, the Helmholtz free energy per unit cell of the MgO lattice at $300\mathrm{K}$ was calculated at various mesh resolutions. The results are summarised in table 2 and illustrated in figures 9 and 10.

Table 2. Helmholtz free energy dependence on mesh resolution.

Mesh resolution	F/eV	$\mathrm{\Delta }F$/eV
1x1x1	$-40.930301$	N/A
2x2x2	$-40.926609$	$0.003692$
3x3x3	$-40.926432$	$-0.000177$
4x4x4	$-40.926450$	$-0.000018$
5x5x5	$-40.926463$	$-0.000013$
6x6x6	$-40.926471$	$-0.000008$
7x7x7	$-40.926475$	$-0.000004$
8x8x8	$-40.926478$	$-0.000003$
9x9x9	$-40.926479$	$-0.000001$
10x10x10	$-40.926480$	$-0.000001$
15x15x15	$-40.926482$	$-0.000002$
20x20x20	$-40.926483$	$-0.000001$
25x25x25	$-40.926483$	$0.000000$
30x30x30	$-40.926483$	$0.000000$
40x40x40	$-40.926483$	$0.000000$

Figure 9: Variation of the Helmholtz free energy with mesh resolution (where a single number Y indicates a YxYxY grid).

Figure 10: Variation of the fractional change of the Helmholtz free energy, $|\frac{\mathrm{\Delta }F}{F}|$ , with mesh resolution.

The most distinctive feature of figures 9 and 10 is the rapid convergence of the fractional change in Helmholtz free energy to zero, and that of the Helmholtz free energy itself to a constant value of -40.926483 eV . To investigate the accuracy at each mesh resoltion it is necessary to compare not this change in Helmholtz free energy, which depends also on the choice of shrinking factors at which the free energy was evaluated, but instead the difference between the final converged free energy (taken to be -40.926483 eV) and the free energy at a particular mesh resolution. A plot of this is shown in figure 11.

Figure 11: Variation of the difference between the converged Helmholtz free energy and that for a given mesh resolution with mesh resolution.

Clearly then, an accuracy of 1 meV and 0.5 meV can be achieved with a 2x2x2 grid, and one of 0.1 meV with a 3x3x3 grid. Whilst the discrepancy between the rapid convergence of the Helmholtz free energy, and the much more obvious changes in the visual appearance of the phonon DOS curve might seem counterintuitive, it is important to note that the energetic contribution (both zero-point and excited-state) from phonons is, particularly at 300 K, very small compared to the lattice energy. Therefore, the Helmholtz free energy can be determined with reasonable accuracy even when the phonon DOS curve does not strongly resemble its converged form.

However, as the calculations were not particularly computationally expensive, a mesh resolution of 40x40x40 was used throughout this investigation, to ensure that a lack of sample points in the Brillouin zone was not a significant contributor to the uncertainty in the final results. This grid size would certainly be appropriate for similar compounds such as CaO. Since the dimension of the first Brillouin zone (BZ) is proportional to $\frac{1}{a}$, where $a$ is the characteristic length of the Wigner-Seitz cell, a material with a considerably larger Wigner-Seitz cell than MgO, such as a Zeolite (e.g. Faujasite), would have a smaller first BZ and thus require fewer sample points over which to evaluate the phonon DOS at similar levels of accuracy, i.e. a smaller mesh resolution. For a metal such as Lithium, another factor comes into play. Namely, the delocalised electrons are able to screen the repulsions between approaching charged ions, resulting in decreased phonon dispersion ('width' of the branches). As a result, the energies of the phonon modes at various k-points are more similar than in an insulator, and so an average over a fewer number of sample points describes the phonon DOS curve with equal accuracy. However, the Wigner-Seitz cell of Lithium is smaller than that of MgO, and so, conversely to the case of a Zeolite, this would tend to increase the mesh resolution required for the same level of accuracy. The question of the balance between these two competing factors is outside the scope of this investigation.

The Thermal Expansion of MgO

Figures 13 and 14 show the behaviour of the lattice parameter (the edge-length of the primitive unit cell of MgO) and the Helmholtz free energy with increasing temperature, as calculated under the QHA.

Figure 13: The Helmholtz free energy as a function of temperature

Figure 14: Dependence of the lattice parameter on temperature

As may be intuitively expected, the lattice parameter, and thus the volume of the unit cell, generally increases with increasing temperature. However, at temperatures below ~100 K, there is no significant thermal expansion. To explain this, it is necessary to consider the behaviour of all three terms in the expression of the Helmholtz free energy:

$F(V,T)=U+\int \mathrm{d}\omega g(\omega )[\frac{\hslash \omega }{2}+k_{\mathrm{B}}T\log (1-e^{-\frac{\hslash \omega }{k_{\mathrm{B}}T}})]$

At low temperatures, the last term can be discarded, as the argument of the logarithm approaches unity and is multiplied by the temperature. The observed lack of thermal expansion - and variation in the Helmholtz free energy - at low temperatures is therefore a result of the temperature-independence of the remaining terms. This equation also provides the physical basis for the observed positive thermal expansion. Generally, an increase in volume decreases the values of $\omega$, decreasing the zero-point energy but increasing the lattice energy. However, it also increases the ability of the system to access higher-lying excited phonon states, the entropic contribution of which is characterised by the last term. At larger temperatures, the contribution of this to the free energy increases, favouring a larger volume.

To compute the coefficient of thermal expansion, a fourth-order polynomial was fitted to the volume-temperature data, and its numerical derivative evaluated, from which the coefficent of thermal expansion $\alpha$ was then calculated. This analysis was performed for both the data obtained through MD simulations and calculations under the QHA, and the resulting curves of the temperature-dependence of $\alpha$ compared to another and to experimental data, as shown in figures 15 and 16.

Figure 15: The volume of the primitive unit cell of MgO as a function of temperature

Figure 16: Dependence of the coefficient of thermal expansion on temperature. Note the excellent qualiative agreement between the results obtained using the QHA and those obtained by Reeber et al.

As can be seen, the MD simulations do not reproduce the low-temperature tailing observed in the other datasets. This can be attributed to the fact that MD simulations operate exclusively in the realm of Newtonian mechanics, and thus it is no surprise that they fail to account for such fundamentally quantum mechanical effects as the zero-point energy. The agreement between the two methods is better between 400-1000 K, as the zero-point energies become less significant. However, as will be shown later, the discrepancy increases again at temperatures approaching the melting point.

An intruiging feature of figure 16 is that both of the theoretical methods used overestimate the thermal expansion coefficient relative to the experimental data. The most significant approximations used in this investigation are the modelling of MgO as a perfect crystal, and the use of the Coulomb-Buckingham potential to simulate the interatomic interactions. Of these two, the latter is most likely to account for the observed discrepancy (particular at the modest temperatures under consideration), since the defect energy of MgO is very large, and the Coulomb-Buckingham potential ignores the partially covalent character of the Mg-O bond, thus underestimating the bond strength.

The Melting of MgO

The QHA is inherintly unable to predict the melting of a crystal lattice, since at each volume and temperature, the vibrational potentials are modelled to be harmonic. There is considered to be no thermal motion, and so, whilst the lattice spacing may change, the structure will always be that of a perfect crystal. Thus, the behaviour of the volume of the crystal around the expected melting point of ~ 3000 K. The resulting data from the MD simulations is shown in figure 17.

Figure 17: The volume of the primitive unit cell of MgO as a function of temperature around the melting range. Of special interest is the onset of large fluctuations ~3020 K

The most notable feature of figure 17 is the onset of large deviations from the otherwise close to linear variation of the primitive unit cell volume with temperature. The similarity of these fluctuations to those observed by Vocadlo et al. , and their occurance at approximately 3025 K, which is very close to the experimental zero-pressure melting range of $3040\pm 100\mathrm{K}$ , indicate that they arise from the melting of the crystal lattice.

Interestingly, attempts to perform calculations using the QHA at temperatures above 2700 K did not produce usable data, as the energy did not converge. However, it can clearly be seen in figure 17 that the unit cell volume calculated under the QHA is significantly larger than that derived from the MD simulations at equivalent temperatures, much more so than the systematic discrepancy present in the range of 100-1000 K.

For all MD simulations, a supercell consisting of 2x2x2 conventional MgO unit cells - and thus 32 MgO formula units - was used with periodic boundary conditions. Particularly at larger temperatures, where thermal motion becomes more sinificant, this may result in significant deviations from the behaviour of a perfect infinite crystal. This presents the need for a compromise between computational expense and accuracy, analogous to the shrinking factors used in the phonon calculations under the QHA. It would therefore have been of interest to compare the obtained results with those from MD simulations of larger supercells. However, time constraints did not allow for this to take place.

The Bulk Modulus of MgO

The volume-pressure data shown in figure 18 was fitted to the BMEOS, yielding the two curves shown. The resulting values of the Bulk modulus and its pressure derivative at 0 Pa are compared to another and the values obtained experimentally by Rekhi et al. in table 3:

Figure 18: The volume of the primitive unit cell of MgO as a function of pressure, as calculated from MD simulations and under the QHA. The data was fitted to the BMEOS.

Table 3. Bulk modulus as calculated by MD simulations and QHA, compared to experimental values.

	This investigation			Rekhi et al.
MD or QHA	$K_0$ / GPa	$K{\text{'}}_0$	$V_0$ / ${\AA }^3$	$K_0$ / GPa	$K{\text{'}}_0$	$V_0$ / ${\AA }^3$
QH	$222.9\pm 0.1$	$3.960\pm 0.004$	$18.8884\pm 0.0003$	$203\pm 16$	$2.65\pm 0.45$	$18.82\pm 0.10$
MD	$223.3\pm 4.6$	$4.06\pm 0.22$	$18.84\pm 0.01$

Where the convidence intervals in the calculated parameters correspond to $\pm$ one standart deviation, and were obtained from the covariance matrices of the numerical fits.

Conclusions

In summary, the behaviour of MgO under both changes in pressure and temperature was investigated. First, the convergence of the phonon DOS - calculated under the QHA - with an increasing resolution of sample points in the first BZ, was assessed both visually and by monitoring of the Helmholtz free energy. A mesh resolution of 40x40x40 was then used to calculated the dependence of the primitive cell volume on the temperature, the results compared to those derived from MD simulations, and both used to compute the coefficient of thermal expansion. It was found that MD simulations were able to reproduce the melting of the MgO lattce, whilst calculations under the QHA were not. The bulk modulus of MgO in the absence of pressure was determined through the use of both MD and QHA-based methods, and the results found to be consistent with experimentally measured values.