Rep:Mod:ym2814mgo

Introduction

The system: MgO

MgO lattice will be studied in this experiment and the key focus is its thermal expansion. Magnesium oxide is one of the candidate constituents of the Earth's lower mantle and the estimation of its thermal properties will lead to a better understanding of the mantle structure.

Methodology

The lattice structure of MgO is a face-centered cubic with lattice parameter a, b and c all equal to each other. A 1x1x1 grid size represents a primitive cell that contains only one MgO and grid size will be varied in this experiment. In order to obtain an accurate approximation for free energy, a reasonable grid size needs to be found out. The larger the grid size, the more accurate the calculations. Phonon dispersion curves and density of states are therefore calculated. The free energy of MgO is computed within the quasi-harmonic approximation (QHA), in which the anharmonicity resulted from the intrinsic phonon interaction can be neglected.

Thermal expansion is a fundamental property of materials and is extremely important when choosing materials for mechanical and structural applications. The same approximation is then used to optimize the structure of MgO lattice under different temperatures and calculate the coefficient of thermal expansion.

The expansion coefficient is defined as:

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

Where $\alpha$ is the expansion coefficient, V₀ is the initial cell volume and T is the temperature. The value of expansion coefficient can be determined from the temperature dependence of cell volume at constant pressure, which will be performed later in the investigation.

In the molecular dynamics simulation, the energy is minimized as a function of the position of the atomic position. The interionic potential is taken to be the sum of pairwise additive Coulomb, van der Waals, and repulsive interactions. The Helmholtz free energy is calculated by the equation showing below:

Equation 2: Helmholtz free energy

Aims

This exercise is to investigate the thermal expansion of MgO lattice by using both quasi-harmonic approximation and molecular dynamics simulations.

Tools in use

The softwares being used for all simulations are RedHat Linux, DLVisualize and GULP.

Results and Discussion

Calculating the internal energy of an MgO crystal

In order to compute the energy per unit cell of MgO, an initial calculation was run at 0 K and 0 Pa. The primitive unit cell of MgO is a rhombohedron, which has cell vectors a = b = c = 2.9783 Angstroms with internal angles 60 degrees. The conventional cell of MgO is a cube with side length of 4.212 Angstroms (lattice constant). According to the log file, the total lattice energy per unit cell is -41.07531759 electron volts, which is the energy required to form a crystal from infinitely-separated ions.

Calculating the phonon modes of MgO

It has been revealed that the vibration of atoms in a crystal is correlated, and a wave of allowed frequencies and amplitude is formed from a collective vibration. "Phonon" is the quantum of such lattice vibrations and "phonon dispersion" gives the wave vector dependence of phonon frequencies . Every possible vibration can be labeled with a k-vector. In order to understand the variation of frequencies with wavevector k, vibrational frequency is plotted against k and phonon dispersion curves can be therefore obtained.

Figure 1: Phonon dispersion curves of MgO

At the bottom of the panel, the special points along the conventional path in k-space, $W-L-\mathrm{\Gamma }-X-W-K$, are displayed, where $\mathrm{\Gamma }$ is the origin.There are 6 branches in phonon dispersion curves of MgO lattice, which corresponds to 2 atoms per cell and 3 dimensions.

In order to compute the free energy of the crystal, all of the vibrational bands need to be summed over, which means summing over all possible k-points. The number of points on a grid should be large enough to obtain more accurate results, however the more the points, the longer the computing time. The density of states (DOS), an average over all k-points yielding the number of vibrational modes at each frequency, is a useful object to summarize phonon dispersion curves and will be used to find the minimum grid size for a reasonable approximation for later calculations.

The density of states for 1x1x1 grid was first computed for a single k-point, which is showing below.

Figure 2: DOS for a 1x1x1 grid

There are 4 distinct peaks with two double in intensity compared to the other two. According to the log file, the vibrational frequencies are 288.49, 288.49, 351.76, 351.76, 676.23 and 818.12 cm^-1. These indicate that two sets of degenerate bands should be seen in the dispersion curves at 288.49 cm^-1 and 351.76cm^-1 respectively. From the dispersion curves showing in figure 1, it can be seen that point L shows the same characteristics.

The following graphs show how DOS varies with grid size.

Figure 3: DOS for a 2x2x2 grid

Figure 4: DOS for a 4x4x4 grid

Figure 5: DOS for a 8x8x8 grid

Figure 6: DOS for a 16x16x16 grid

Figure 7: DOS for a 32x32x32 grid

Figure 8: DOS for a 64x64x64 grid

As the shrinking factor increases, the DOS becomes smoother due to more possible vibrations are sampled and more features appear. The 32x32x32 grid is chosen as the minimum grid size for a reasonable approximation because there is no much change between the 32x32x32 grid and 64x64x64 grid and the former one takes a shorter time to compute.

Calculating the Free Energy in the Harmonic Approxiamation

The density of states is calculated for testing molecular set and, based on it, the characterization and prediction of thermal properties of materials are evaluated within the quasi-harmonic approximation (QHA). In order to justify the choice of 32x32x32 grid size from a numerical point of view, the free energy of MgO was computed within quasi-harmonic approximation and the same shrinking factors are used as last part.

Shrinking Factors	Helmholtz Free Energy /eV	Total Lattice Energy /eV
1x1x1	-40.930301	-41.07531759
2x2x2	-40.926609
4x4x4	-40.926450
8x8x8	-40.926478
16x16x16	-40.926482
32x32x32	-40.926483
64x64x64	-40.926483

As grid size increases, the free energy first decreases dramatically from -40.930301 eV to -40.926609 eV, which is followed by some fluctuations before convergence. Smaller changes in free energy values can be seen for larger grid sizes and there is no difference in free energy between shrinking factor of 32 and 64, which is the evidence for complete convergence. Therefore 32x32x32 grid size is accurate enough for approximation.

According to the energy differences, 2x2x2, 4x4x4 and 8x8x8 grid size can be used for calculations accurate to 1 meV, 0.5 meV and 0.1 meV respectively.

The Thermal Expansion of MgO

In order to calculate the coefficient of thermal expansion, the structure of MgO will be optimised with respect to free energy from 0 K to 1000 K in steps of 100 K. The structure with minimum energy will be calculated at each temperature.

Figure 9: Thermal expansion simulation at different temperatures

The plots of free energy, lattice constant and cell volume against temperature are shown below:

Figure 10: Free energy against temperature

Figure 11: Lattice constant against temperature

Figure 12: Cell volume against temperature

It can be seen that free energy decreases with temperature and lattice constant increases with temperature. As for the cell volume, it increases linearly from 300 K to 1000 K with gradient equals to 0.0005. The initial cell volume is 18.836496 ų. Then the expansion coefficient can be calculated by equation 1.

$\alpha =\frac{1}{1.836496}\times 0.0005=2.65\times 10^{-5}{K}^{-1}$

In literature, the expansion coefficient at 300 K is 3.06 x 10^-5 K^-1, which is very closed to the computed value.

Molecular Dynamics

A larger cell that contains 32 MgO units was used for molecular dynamics simulation. The system is allowed to evolve in time according to Newton's second law F=ma. The time step used was 1 femtosecond and the temperatures vary from 100 K to 1000 K. The cell volume was calculated as the average of last 5 timesteps.

Figure 13: Cell volume against temperture for both QHA and molecular dynamic simulation (MD)

The QHA and MD share the same value of expansion coefficient due to the same gradient (0.0005). It seems that QHA is more accurate when the temperature is below 300 K because zero-point energy is taking into account. However, MD will be a more reliable method for higher temperatures.

Conclusion

The grid size will affect the accuracy of the calculation. The accuracy would be perfect for an infinite grid but the calculation would take an infinite amount of CPU time to run. By performing the phonon dispersion curves and density of states, 32x32x32 grid size was chosen for quasi-harmonic approximation and a supercell containing 32 MgO units was chosen for molecular dynamics.

The thermal expansion of MgO lattice was simulated by both QHA and MD and comparison has been made between two methods. Same expansion coefficient has been obtained from the two methods, which is 2.65x10^-5 K^-1. However, QHA will be more reliable at relatively low temperatures while MD is more accurate for higher temperatures.