Rep:Mod:Bnn13

Introduction

Magnesium oxide is one of the most popular natural compounds and is mainly found in minerals. The compound is made up from the strong ionic bonding between Mg²⁺ and O^2- ions, therefore it is stable and is usually a typical example used for lattice structure studies. Primitive lattice cell of MgO displays a triclinic structure. Each primitive lattice cell contains two atoms. On the other hand, the conventional cell shows a perfect cubic structure, containing eight atoms (four magnesium atoms and four oxygen atoms). The initial volume of a cell is 18.6804 ų with total lattice energy of -41.0753 eV. The conventional lattice cell is widely chosen for study due to its better symmetry. The structures of conventional, primitive lattice cell and supercell are demonstrated in Figure 1a-c.

Figure 1a. The conventional strucutre of MgO	Figure 1b. The primitive cell of MgO	Figure 1c. The supercell of MgO

The Fourier transformation of real space lattice is represented by reciprocal space. The reciprocal space shows significant advantages in studying atomic arrangement by relating to X-ray diffraction. The primitive cell of reciprocal space, Brillouin zone, provides fundamental information about the structure. Therefore, the phonon dispersion curve and density of states of MgO were measured in the first Brillouin zone.

In this study, free energy and thermal expansion of MgO were investigated using two methods, quasi-harmonic approximation (QHA) and molecular dynamic simulation (MD). Quasi-harmonic approximation is built up on the assumption that the vibrational atoms behave as harmonic motion neglecting phonon-phonon interaction. The Helmholtz free energy of the system was calculated using Equation 1.

$F=Eo+\frac{1}{2}\sum _{j,k}\hslash {\omega }_{j,k}+k_{B}T\sum _{j,k}ln[1-exp(\frac{-\hslash {\omega }_{j,k}}{k_{B}T})]$ (Equation 1)

Represented in the above equation, the first term of the equation, $Eo$, is lattice potential energy at absolute zero and the sum $\frac{1}{2}\sum _{j,k}\hslash {\omega }_{j,k}+k_{B}T\sum _{j,k}ln[1-exp(\frac{-\hslash {\omega }_{j,k}}{k_{B}T})]$ is responsible for the vibrational contribution to total energy. The zero-point energy term $\frac{1}{2}\sum _{j,k}\hslash {\omega }_{j,k}$ is the summation of the permissible quantum states resulting in the normal mode energy, while $k_{B}T\sum _{j,k}ln[1-exp(\frac{-\hslash {\omega }_{j,k}}{k_{B}T})]$ is the contribution of thermal energy. When temperature is infinitesimally small, the thermal energy term approaches zero and becomes negligible. This results in the temperature-dependence of the accuracy of the calculations using QHA that will be discussed in the subsequent sections.

Conventionally, the expression for free energy of the system adopts a time-dependent quadratic function in QHA because the higher order polynomial terms (obtained from Taylor’s expansion) are usually neglected.

The physical origin of thermal expansion is from the harmonicity of interatomic potential. When temperature increases, the crystal lattice expands lead to the increase in vibrational amplitude. Corresponding to the anharmonic potential, the rising of average interatomic separation results in thermal expansion. On the other hand, the average interatomic separation remains constant in harmonic potential, therefore, there is no thermal expansion.

Linear and volumetric thermal expansion coefficient (${\alpha }_L$ and ${\alpha }_V$, unit $K^{-1}$) of MgO were identified using Equation 2 and 3 respectively, in which $\frac{\partial L}{\partial T}$ and $\frac{\partial V}{\partial T}$ are the first derivative of the temperature-dependent function of length and volume. $L_o$ and $V_o$ corresponds to the initial length and volume of material at temperature of absolute 0 K.

${\alpha }_L=\frac{1}{L_o}\frac{\partial L}{\partial T}$ (Equation 2)

${\alpha }_V=\frac{1}{V_o}\frac{\partial V}{\partial T}$ (Equation 3)

The second method of study was molecular dynamic simulation. In this method, the classical Newtonian mechanic is applied as the fundament. All the atoms interactions are allowed to take place in a fixed short time scale, which is one femtosecond in this study. The time average values of forces, charge and position of atoms were obtained and contributed to the calculations of energy and volumetric thermal expansion. In order to perform a valid MD simulation, a large system of MgO consisting eight formula units (2x2x2 MgO cells in three dimentions) corresponding to 32 MgO units were required. This is because the investigation of a single MgO primitive cell will be similar to the sistuation that all the atoms are in phase. Therefore, the simulation of a larger cell size is required as it permits representative information about the system. However, the limitation incurred by this is a long operating time. In contrast to QHA, MD does not consider the contribution of lattice energy at absolute zero. This will lead to the difference between two methods in calculating thermal properties of MgO compound. Both calculations were operated in RedHat Linux background using DLVisualize (DLV) and General Utility Lattice Program (GULP) software. DLV provides an excellent environment for quantum mechanic related modellings while GULP is suitable for classical simulation. In this study, QHA involves quantum related calculations and MD involves classical related calculation, hence both DLV and GULP are needed.

Results and Discussion

Phonon dispersion curves

Figure 2. Phonon dispersion curve of MgO in the first Brillouin zone

The surface phonon dispersion curves for MgO was calculated in the Brillouin zone with the N-point of 50 at a pressure of 0 Pa. The phonon dispersion curve (Figure 2.) shows six branches in the first Brillouin zone (-π/a ≤ k ≤ π/a), corresponding to six energy states. Because MgO contains two atoms in its primitive lattice lattice cell, its dispersion curve displays both acoustic branches and optical branches. The acoustic branches are originated from the in-phase vibration of the neighboring Mg and O atoms. In contrast, out-of-phase vibrations of the atoms results in optical branches, thus induces the polarisation within the lattice bonds. Indeed, this oscillation can be generated by infrared excitation. As illustrated in the phonon dispersion curve, the acoustic branches converge to the origin $\mathrm{\Gamma }$, while the optical ones do not. The symmetry points are W(1/2, 1/4, 3/4), L(1/2, 1/2, 1/2), $\mathrm{\Gamma }$ (0, 0, 0), X(1/2, 0, 1/2), W(1/2, 1/4, 3/4) and K(3/8 3/8 3/4)

The phonon density of states

While the phonon dispersion curve is the plot of frequency against the numbers of energy states summing over all k-points, the phonon density of states are in fact the inversion plot of dispersion curve as they display the numbers of accessible vibrational states per unit energy respected to each k-point. Several density of states are shown as the following.

Figure 3a. DOS with shrinking factor 1x1x1	Figure 3b. DOS with shrinking factor 2x2x2	Figure 3c. DOS with shrinking factor 4x4x4	Figure 3d. DOS with shrinking factor 8x8x8
Figure 3e. DOS with shrinking factor 16x16x16	Figure 3f. DOS with shrinking factor 32x32x32	Figure 3g. DOS with shrinking factor 64x64x64

As mentioned above, the dispersion curve shows distinctively six different branches, therefore, six peaks of frequency are expected for the DOS of each k-point. However, with the shrinking factor of 1x1x1 and k-point value of one, the DOS graph showed only four different peaks of frequency (Figure 3a). Therefore, there were two doubly degenerate energy levels. In fact, the frequencies at this state were 288.49, 288.49, 351.76, 351.76, 676.23 and 818.82 cm^-1. By examining the dispersion curve, the respective k-point value was L (0.5, 0.5, 0.5).

For the first few DOS curves (Figure 3a-d.) of low shrinking factor, the DOS curve appeared as the combination of a few discrete peaks at different frequencies. With the increase of the grid number, the number of k-point also increases. Because DOS was calculated as an accumulation frequency peaks of all the k-points, the rise of number of k-point resulted in the increase of vibration mode and hence led to a smoother and more accurate DOS curve (Figure 3e-g.). A peak maximun shows at a frequency of 400 cm^-1 in every DOS plots indicating a majority number of accessible vibration states corresponding to that frequency.

The minimum grid size of 16x16x16 was required for a good approximation of the density of states. The DOS of smaller grid size displayed a prominent difference compared to the DOS of grid size 16x16x16, whereas the larger ones show neglectable changes. Hence, the shrinking factor of 16x16x16 was expected to show a great balance between the accuracy and time efficiency of subsequent calculations.

This optimal grid size of 16x16x16 does not only work well for MgO but it is also suitable to approximate the DOS of other compounds which have similar lattice structure (face centred cubic) such as NaCl and CaO. Especially for the case of CaO, its lattice size is 4.815 Å at temperature of 298 K, which is similar to the value of MgO, 4.212 Å, at the same temperature. In addition, the ionic bond strength of CaO is similar to that of MgO as its lattice structure is made of ions Ca²⁺ and O^2-. Therefore, the grid size of 16x16x16 is expected to work well for DOS calculation of CaO.

However, this grid size might not be appropriate for the DOS calculation of other metal such as sodium, lithium, etc. Unlike MgO, these metals display different crystal structures (body centered cubic). Moreover, metals, like lithium, have their ions freely flow in a pool of delocalised electrons, which efficiently damp the vibration of metallic ions (Li⁺) resulting in a flatter dispersion curve. Therefore, a smaller grid size might be required to obtain an optimum calculation.

For more complex systems such as Zeolite family, the compounds contain the complexes of a few ions (for instance Faujasite contains Na⁺, Mg²⁺, Ca²⁺ ions). The complexes often have larger unit cell (a = 24.7 Å for Faujasite ), thus smaller reciprocal lattice space ($a^{*}=\frac{2\pi }{a}$, $a^{*}$ and $a$ are the reciprocal and real lattice parameters respectively) . As a result, a smaller grid size might be required to achieve a good approximation of DOS.

Quasi-harmonic approximation of free energy and thermal expansion

Helmholt free energy and thermal expansion of MgO were firstly idetified using QHA. The calculations were run with different grid sizes at temperature of 300 ^oC and pressure of 0 GPa in order to determine the optimun shrinking factor for subsequent calculations. A plot of free energy against grid size is shown in Figure 4.

The energy was fluctuated for the first three shrink factors. With the subsequent values of grid sizes, the energy constantly decreases (more negative) and converged to a constant value of -40.926483 eV (Figure 4). The grid size size of 3x3x3 allowed the calculation to be accurate to 1 meV (-40.926432 eV). An accuracy of at least 0.5 meV per cell was engendered by the minimum grid size of 4x4x4 (-40.926450 eV). Furthermore, a grid size of 5x5x5 was required for an accuracy of 0.1 meV (-40.926463 eV).

Similar to the calculation of DOS, a same grid size as for MgO is required for the energy calculation of CaO and a smaller one is needed in the cases of both lithium and a Zeolite because of the reasons discussed in the previous section.

In addition, a graph of entropy against grid size was generated (Figure 5) excluding the first there values of grid size due to their fluctuation. While the free energy decreases and converts to a constant, the entropy rapidly increases with the grid size at first and slowly converged to the value of 23.3634 J/(mol.K). This observation was in a good agreement with the fact that with a increasing number of k-point, there are more available energy states and, therefore, a more disorder system.

Figure 4. Helmholtz free energy curve respected to grid size	Figure 5. Plot of entropy against grid size

A minimum shrinking factor of 18x18x18 resulted in an energy of -40.926483 eV, close to the converged value. Therefore, it was chosen and expected to engender sufficiently good calculations for thermal expansion. A series of calculations of energy per fomula unit and lattice constant were carried out with temperatures varying from 0 K to 1000 K. The plots of free energy and lattice constant respected to temperatures are displayed in Figure 6 and 7 respectively.

Figure 6. Energy plot respected to temperature	Figure 7. Lattice constant plot respected to temperature

The plot of free energy against temperature displays a concave curve, demonstrating a decline of energy respected to temperature (Figure 6). This greatly coincides with the theory shown by the equation A = U –TS. With the positive value of entropy, increasing temperature resulted in a more negative value of Helmholtz free energy. On the other hand, the plot of lattice constant against temperature is convex, demonstrating the growth of lattice size with temperature.

The coefficient thermal expansion of MgO was calculated using the plot of lattice constant against temperature (Figure 7). The plot of lattice constant vs. temperature shows the best fit curve with the polynomial function of $y=10^{-8}\times x^2+9\times 10^{-6}\times x+2.9858$. Because linear thermal expansion is related to lattice constant as stated in Equation 2, the first derivative of the best fit line was identified, which was as $\frac{\partial L}{\partial T}=2\times 10^{-8}\times x+9\times 10^{-6}$. The value of $L_o$, which is indeed the y-intercept of the curve, is shown to be 2.9858 in the plot. Hence, the coefficient of thermal expansion is obtained as an temperature-dependent equation, ${\alpha }_L=\frac{1}{2.9858}(2\times 10^{-8}\times x+9\times 10^{-6})$, in which the thermal expansion coefficient ${\alpha }_L$ has the unit of $K^{-1}$ and the temperature $x$ has the unit of $K$. As a result, respected to a low temperature of 130 K, the calculated linear coefficient of thermal expansion was $3.8850\times 10^{-6}{K}^{-1}$, which greatly agreed with the reported value of $3.86\times 10^{-6}{K}^{-1}$ measured at a same temperature by Schulz et al.. At high temperature, 2000 K for instance, the linear coefficient value was noted as $1.6411\times 10^{-5}{K}^{-1}$. This value was in fact considerably higher than the recorded one at the same temperature, $12.6\times 10^{-6}{K}^{-1}$

At high temperature, especially near MgO melting point, which is at $3250\pm 20K$, the calculated value of linear coefficient of thermal expansion displayed a significant deviation from reported values. These deviations were due to the contribution of prominent phonon-phonon interaction resulted in anharmonic behaviour of MgO phonon modes. Moreover, the main assumption of quasi-harmonic is that all the atoms were assumed not interacted and their vibrations were treated as harmonic motion. Consequently, the bond length would be expected to increase infinitely with temperature without taking into account of bond dissociation. In reality, bond dissociation will occur when temperature exceeds a certain value. Therefore, quasi-harmonic approximation failed to estimate accurately the energy, lattice parameters and thus thermal expansion at high temperature.

Molecular dynamic simulation of free energy and thermal expansion and their comparison with QHA

The calculation of MgO thermal expansion was repeated with molecular dynamics method. The calculations were carried out using supercell structure containing thirty two MgO unit cells because the simulation of one single primitive cell would simply show an in-phase vibration. An isothermal-isobaric ensemble, in which the number of particle (N), pressure (P) and temperature (T) of the system were retained constant, was investigated in this measurement. Volume (V) of the system was allowed to vary in each calculation. A number of calculations were carried out at different temperatures. For a valid comparison between MD and QHA, the energy and volume per formula unit were obtained by dividing the reported values from MD method by thirty two.

Figure 8. The plot of energy against temperature obtained by MD	Figure 9. The combined graphs of energy vs. temperature by both QHA and MD

The free energy per unit cell obtained by MD was directly proportional to temperature as shown in the Figure 8. The best fit line was represented as the equation $y=0.0005\times x-41.092$. This affirmed the linear relationship between energy and temperature in molecular dynamic simulation taking the assumption of Newtonian mechanic behavior of atoms interaction. The combined graphs of energy against temperature achieved by both QHA and MD is presented in Figure 9. Compared to the values calculated by QHA, which constantly decreased with temperature, the values by MD linearly increased when temperature raised. The difference became more significant at higher temperature. At room temperature, the calculated energy values achieved from QHA and MD methods were - 40.9532 eV and - 40.943 eV, respectively. These calculated values demonstrated a good agreement between QHA and MD and coincide with the literature value of - 40.7233 eV

Furthermore, the calculations of lattice volume by MD was carried out by molecular dynamic method. The best fit line of volume against temperature was represented by a linear equation of $y=0.0006\times x+18.625$ (Figure 10). The combined graph of volume against temperature resulted from both MD and QHA is shown in Figure 11. The lattice volume shows a linearly increasing trend with temperature under MD and a parabola trend under QHA. The volume per formula unit of MgO approximated by QHA demonstrates a great agreement to that by MD in a range of temperature of 500 - 1300 K. The deviation between two methods occurs in the lower range of temperature (below 500K) and in the higher range of temperature (above 1300 K). The difference become more noticeable as temperature further increases.

Figure 10. The plot of volume against temperature obtained by MD	Figure 11. The combined graphs of volume vs. temperature by both QHA and MD

By applying Equation 3, the volumetric thermal expansion of MgO was calculated as a constant of $3.22148\times 10^{-5}{K}^{-1}$. A similar calculation was carried out for the volumetric thermal expansion coefficient by QHA in a range of temperature of 100 - 2900 K. The volumetric thermal expansion coefficient was obtained as a function of ${\alpha }_V=\frac{1}{18.841}(4\times 10^{-7}\times x+2\times 10^{-4})$. At the temperature of 273 K, the volumetric coefficients of thermal expansion achieved from QHA and MD were $1.1104\times 10^{-5}{K}^{-1}$ and $3.22148\times 10^{-5}{K}^{-1}$ respectively. Compared to the literature value of $3.15\times 10^{-5}{K}^{-1}$, MD showed a better performance in calculating volumetric coefficient of MgO compared to QHA.

Figure 12. The combined volumetric coefficient of thermal expansion by both QHA and MD

While thermal expansion was temperature dependent and increased linearly with temperature in quasi-harmonic approximation, it was temperature independent in molecular dynamic (Figure 12). The volumetric coefficents of thermal expansion calculated by QHA and MD showed a fair agreement in a small range of temperature, 1100 - 1400 K. A significant deviation in the calculated values of ${\alpha }_V$ by these two methods occured in lower and higher temperature regimes.

The difference in term of energy, thermal expansion and lattice volume obtained by two methods were the consequence of assumptions taken by QHA and MD. It is noteworthy that QHA only work well when the vibration is near equilibrium position, which is at the bottom of potential energy well, because it takes the assumption that atoms vibration as simple harmonic motion. Because QHA takes into consideration of zero-point energy, therefore, it make a good approximation for the values at low temperature. At high temperature, bond dissociation, which is neglected in QHA approximation, occurs. Furthermore, the phonon-phonon interaction becomes significant and contributes to thermal energy, QHA is no longer accurate in the calculation of MgO thermal properties. On the other hand, MD adheres to Newton’s laws and ergodic hypothesis. Due to this classical approximation, the total free energy and lattice volume were linearly proportional to temperature. Though MD was a good replacement for QHA to measure thermal expansion of lattice structure at high temperature, it showed limitation for the calculations at low temperature because Newtonian mechanic does not take into consideration of zero-point energy. Beyond the melting point temperature of the compound, both methods were not appropriate to approximate the thermal properties of MgO because the compound is no longer in solid state.

Conclusion

The calculations of thermal expansion and lattice volume respected to temperature by two method quasi-harmonic approximation and molecular dynamic simulation agreed to each other and coincided with literature value in the temperature range from about 500 – 1300 K. At lower temperature QHA demonstrated more accurate approximations at it considers the contribution of zero energy and quantum mechanics behavior of the system while MD dose not. On the other hand, at high temperature regime, QHA no longer accurately approximates the thermal properties of the system the phonon-phonon interaction became more significant. Moreover, bond dissociation, which QHA does not account for, might become an issue when the vibrations occured at too high temperature. In this case, the molecular dynamic became a better approximation method as it takes into consideration the contribution of all the atom interactions. Nearer to the MgO melting point and higher, the compound is no longer in its perfect solid state, therefore, both methods are no longer suitable to measure the thermal properties of the material.