Mod:MgOjh6215

Thermal Expansion of MgO

Written by James Hale.

Introduction

The Structure of MgO

Figure 2 - 2x1x2 Primitive cell of MgO

Figure 1 - Conventional Cell of MgO

The MgO lattice is a face centred cubic (fcc) ionic lattice, which is represented in three different ways during this experiment. The conventional cell (figure 1), the primitive cell (figure 2) and the super cell. The conventional cell contains 4 MgO units and has the lattice structure: a=b=c, α=β=γ=90°. It has the same density but 4 times the volume of the primitive cell, which only contains 1 MgO unit and has the lattice structure: a=b=c, α=β=γ=60°. Even though the primitive cell contains fewer MgO units, the conventional cell is more widely used due to the simplicity of its symmetry. The supercell is used for the molecular dynamics calculations and is the conventional cell doubled in all 3 dimensions to make a 2x2x2 cell. This has 8 times the number of atoms and the lattice parameters defining the cell are twice as large (2a_Conventional = a_Supercell).

It is important to note that the computer will output the results for a primitive cell on completing calculations using the quasi-harmonic approximation. Therefore, the results will have to be converted to those of a conventional for direct comparison with other data.

Quasi-harmonic Approximation

The quasi-harmonic approximation is an improvement on the harmonic approximation for modelling lattice dynamics. For the harmonic approximation the atoms in the lattice only interact with their nearest neighbours. The energy of a 1 dimensional monatomic chain can be modelled by the following Taylor expansion of the interactions between the nearest neighbours:

$E = N\phi + \sum_{s\geq 1} \frac{1}{s!} \frac{\partial^s \phi}{\partial u^s} \sum_n (u_n - u_{n+1})^s ~~ Equation~1$

The equation states that the energy of a 1 dimensional monatomic chain, $E$, is equal to: the rest energy of every atom; where $N$ is the number of atoms and $\phi$ is the energy between 2 neighbours as a function of the distance between them, and the energy of the vibrations within the monatomic chain. The vibrational energy is equal to the Taylor series for the displacement of each atom from its rest position, where: $s$ is the order of the Taylor series, $u_n$ is the displacement of the nth atom and $u_{n+1}$ is the displacement of the nth +1 atom. This represents the correction of the energy between 2 neighbours from the rest energy where the distance between two neighbours is equal to the lattice constant, $a$ . For this equation to be formed the displacement of each atom, $u$ , must be small in comparison to the lattice constant, $a$ otherwise the series may not converge. Usually $u$ is much smaller than $a$, and so a quadratic approximation can be applied ($s=2$). This approximation can be made not just because the higher order terms are small and can be ignored but because harmonic equations have exact solutions, whereas anharmonic equations give numerical solutions. The additional level of accuracy is not required in this investigation and so using the anharmonic method would lead to an unnecessarily high computational cost. Therefore, the harmonic energy of a 1 dimensional monatomic chain is given by equation 2:

$E^{harmonic} = \frac{1}{2} \frac{\partial^2 \phi}{\partial u^2} \sum_n (u_n - u_{n+1})^2 ~~ Equation~2$

To solve the problem of having open ends of the monatomic chain, the Born-von Karman periodic boundary condition is imposed. Harmonic motion along the continuous chain is sinusoidal and the motion of the entire system can be modelled by travelling waves. The time dependent motion of the nth atom is dictated by the linear superposition of the traveling waves and is shown by equation 3:

$u_n(t) = \sum_k u_k exp(i[kx-\omega_kt]) ~~ Equation~3$

Where $u_n$ is the displacement of the nth atom, $t$ is time, $k$ is the wave vector and is equivalent to $\frac{2\pi}{Wavelength}$, $\omega_k$ is the angular frequency and $u_k$ is the amplitude. $x$ is the displacement of the wave can only be equal to $na$ where $n$ is an integer. The wave equation is for each wave can be solved if the wave vector corresponding to that wave is considered individually. This can be simplified to equation 4:

$\omega_k = \left(\frac{4J}{m}\right)^{1/2} |sin(ka/2)| ~~ Equation~4~~~~$

Where $J$ is the force constant $J = \frac{\partial^2 \phi}{\partial u^2}$ and m is the atomic mass. This equation models the relationship between phonon frequency, proportional to the energy, and the wave vector. The equation can be used to the plot dispersion curve, which is the frequencies of phonon vibration at each k-point (figure 3). However, reciprocal space, which is equivalent the Fourier transform of the original lattice, must be implemented to plot the dispersion curve. The relation between the lattice constant in real and reciprocal space is:

$a* = \frac{2\pi}{a} ~~ Equation ~5$

Using reciprocal space allows for the wave vectors in higher dimensional lattices to be calculated. The density of states (table 1) at each k-point on the dispersion curve can be found and summing over all of the k-points, all of the vibrational contributions to the free energy can be found. Combined with the quasi-harmonic theory, reciprocal space can therefore facilitate the calculations of force constants and energy between atoms in a three dimensional lattice, ultimately leading to the calculation of the Helmholtz free energy of the lattice in equation 6 below:

$A = E_0 + \frac{1}{2}\sum_{k,j}\hbar\omega_{k,j}+ k_bT\sum_{k,j}ln[1-exp(-h\omega_{k,j}/ k_bT)]~~ Equation~6$

Where $A$, is the Helmholtz free energy, $E_0$is the internal energy of the lattice, $\frac{1}{2}\sum_{k,j}\hbar\omega_{k,j}$ is the vibrational ground state energy of the lattice and the final term is the contribution from higher vibrational states. The quasi-harmonic approximation is able to calculate the free energy of the lattice accurately because it is able to consider the contributions from the internal energy, which is trying to make the crystal smaller and the entropy, which is trying to increase the distance between the atoms to reduce the frequency of phonon vibrations. The quasi-harmonic approximation makes the assumption that all of the anharmonicity is contained within the resulting thermal expansion and therefore can assign volume dependent harmonic approximations to each lattice volume. This works well at lower temperatures because the higher order terms from equation 1 remain small. However, at high temperatures the anharmonic effects become larger and the higher order terms begin to make the quasi-harmonic approximation diverge from reality.

For diatomic lattices, the dispersion plots contain two types of branches. The first is the acoustic branch, where the different neighbouring atoms vibrate in phase and the second is the optic branch, where the different neighbouring atoms vibrate out of phase. The optic vibrations are of higher energy and are close to the optical part of the electromagnetic spectrum.

Molecular Dynamics Method

Molecular dynamics is another method that can be used to study the thermodynamic properties of a lattice. The quasi-harmonic approximation is unable to account for anharmonic effect or to model phase transitions. The molecular dynamics method can be used to circumvent these flaws. The molecular dynamics simulation involves solving Newton's laws of motions for a set of atoms in a lattice, which interact using a preset anharmonic potential. The simulation uses small time steps, on a femtosecond time scale in this case, to model the velocity and position of each atom. The forces applied on each atom from its neighbours will cause the velocity of the atom to change, via Newton's second law, and its position and velocity are updated in the next time step. A classical trajectory for each atom will be established at the end of the simulation and these are evaluated to identify the properties of the lattice. A set of velocities is assigned randomly to each atom in the lattice at the start of the simulation, with the overall energy of the atoms made to equal the effect of the set temperature. The simulation must be run for long enough for the starting velocities to be equilibrated across the lattice and the time step must be small enough so that an error does not occur in the lattice. Increasing both the total simulation time and decreasing the time step both lead to a higher computational cost so a balance with the desired accuracy should be found. It must be noted that this is a purely classical theory and that quantum effects will lead to errors in the simulation.

Thermal Expansion

Thermal expansion is an entropy driven process and is the tendency of a substance to increase in volume with increasing temperature. It is modelled by equation 7 below:

$\alpha_V = \frac{1}{V_0} \left(\frac{\partial V}{\partial T}\right)_P ~~ Equation~ 7$

Where $\alpha_V$ is the coefficient of thermal expansion, $V_0$ is the initial volume and $\left(\frac{\partial V}{\partial T}\right)_P$ is the partial derivative of volume with respect to temperature at constant pressure. This equation can be evaluated graphically with the gradient of a temperature against volume graph equalling the expansion coefficient after it is correct by the initial volume term. Thermal expansion occurs because a larger volume reduces the frequency and thus the energy of vibrations. The vibrations increase in frequency at higher temperatures and the volume of the substance increases to account for their contribution to the free energy.

Thermal expansion is very important when considering the construction of buildings and machines. The thermal expansion may cause a structure to collapse or be needed to allow a mechanism to work more efficiently, like in combustion engines.

Software Used

Linux RedHat 5 is a free operating system. It was chosen to run DL Visualise because it is more compatible than Windows with the software. DL Visualise was the graphical user interface chosen for this investigation. It allows for the user to construct and visualise crystals with ease. DL visualise allows for the execution of GULP (General Utility Lattice Program) calculations, which is the program of choice for this investigation's computations. GULP is a purpose built program for analysing lattices and has a reduced computation time due to a more efficient use of lattice symmetry.

Experimental

All parts of the investigation were run on a linux operating system, using a DLV interface to observe the structure and properties of the MgO lattice. The lattice structures to be investigated were downloaded from 2 premade files: one conventional cell and one supercell. All calculations were run using GULP. To prevent excessive computational cost the shells were modelled as hard spheres of charge rather than shells of charge, which can be polarised and so Include Shells option was not selected.The pressure was always set to 0 Pa.

Calculating the Phonon Modes of MgO

The conventional cell, which can be simplified to the primitive cell, was used for this part of the investigation because calculations run with the conventional cell require much less computational power than with the supercell. The phonon dispersion curve was calculated first. The phonon energies were found at 50 points along the k-space between the k-points W-L-Γ-X-W-K. To achieve this the phonon dispersion was computed, with phonon eigenvectors calculated. The temperature was set to 300 K and the Npoints were set to 50. Some phonons were then visualized to check a successful calculation had been completed.

The phonon density of states was then found for k-point grid sizes ranging from 1x1x1 to 64x64x64 at 300 K. The grid sizes were adjusted using shrinking factors A,B and C. For example, when A=1, B=1 and C=1 a k-point grid is formed with dimensions 1x1x1.

Computing the Free Energy - The Quasi-harmonic Approximation

The log files produced during the density of states calculations were examined to find the free energy given by each grid size. The larger the grid size the more accurate the free energy calculation. The ideal grid size is a balance between computational cost and the appropriate level of accuracy. The converged value for the free energy and the ideal lattice size were found.

The Thermal Expansion of MgO

The thermal expansion of MgO was investigated using the quasi-harmonic approximation with the conventional cell and then the molecular dynamics method with the supercell.

For the quasi-harmonic approximation a 24x24x24 grid size was found to be the ideal grid size to provide the best balance between accuracy and computational cost. The density of dates calculation was still active in order to allow for the larger grid size. The conventional cell was used to maintain a low computational cost. The Gibbs energy of the MgO lattice was then optimised over the temperature range 0 K to 2900 K in 100 K steps. Even though the melting point of MgO is 3125 K, the quasi-harmonic approximation would give an error for temperature at 3000 K or above.

For the molecular dynamics method the supercell was used. Even though it requires greater computational cost, the molecular dynamics approximation would not work with a conventional cell because the every cell in the lattice would be moving in phase. This is not a valid representation of the MgO structure and therefore a larger starting cell must be used. The molecular dynamics calculations were run with pressure, the number of particles and temperature all fixed. The time step between each recorded position of the individual atoms was set to 1 fs and the number of equilibration steps, the number of steps allowed for the system to reach equilibrium, was set to 500. The production steps, the number of steps run after the system has reached equilibrium, was also set to 500. The number of sampling steps, the number of steps over which the means are taken, and trajectory write steps, the number of steps between which the geometry is saved in order to produce animations, were both set to 5. The temperature range for the calculations was 100 K - 3300 K with 100 K steps. A 50 K reading was also taken because the simulation did not work at 0 K. The vibrations were animated after each calculation to check that it had run correctly.

Results and Discussion

Calculating the Phonon Modes of MgO

Table 1 - Contains the density of states plots for varying grid sizes

Grid Size	Phonon Density of States Plot
1x1x1
3x3x3
24x24x24
32x32x32
64x64x64

Figure 3 - The phonon dispersion curve of MgO

The phonon dispersion curve is found by calculating all possible frequencies at different k-points between W-L-Γ-X-W-K. The density of states plot is the number of different frequencies and their probability of occuring at a single k-point in reciprocal space. The density of states plot for a 1x1x1 grid is given in table 1. By looking at the frequencies present (around 280, 350, 680 and 820 cm^-1) and the probability of each frequency occurring (2:2:1:1 respectively) it is possible to match the 1x1x1 density of states plot to point L, which has the k-point (0.5,0.5,0.5). This result was verified by checking the relevant log file. The 4 acoustic bands become doubly degenerate and account for the doubling of the probability of the 280 cm^-1 and 350 cm^-1 peaks with respect to the two optical peaks. There are twice as many acoustic modes as optical modes in this case because the acoustic modes can be both longitudinal and transverse phonons whereas the optical modes can only be transverse.

The density of states plots can be seen in table 1 below. Increasing the grid size increases the number of k points sampled, more phonons are displayed and distinct peaks become continuous curves. The number of phonons follows the 3N - 6 rule, where N is the number of atoms. The 6 is subtracted for the translational and rotational modes of the lattice. This is verified by the 1x1x1 grid size because there are 2 atoms present, 1 magnesium and 1 oxygen, and the number of phonons shown is 6. For the larger grid sizes it is difficult to see the individual phonons but the same trend continues. Once a large enough grid size is used the density of states plots no longer changes to a significant degree. This is illustrated in the difference between the 24x24x24, 32x32x32 and 64x64x64 grid sizes.

The minimum reasonable grid size to approximate the density of states was chosen to be 24x24x24. Choosing the ideal grid size is considering the computational cost with the increase in level of accuracy. The 24x24x24 grid size was chosen because there are no significant changes between the density of states at this grid size and those of 32x32x32 (the next largest grid size) and 64x64x64 (the largest grid size used). The optimal grid size of 24x24x24 is very similar to that of larger grids but much smoother and more detailed than the smaller grids of 3x3x3 and 1x1x1. The density of states is at a single k coordinate. The dispersion curve shows a variety of k coordinates and how they change periodically. Therefore, the density of states of a k-point is a single x axis value on the dispersion curve.

Computing the Free Energy - The Quasi-harmonic Approximation

Figure 4

Table 2

Grid Size	Free Energy / eV	Difference Between Measured and Converged Values
1x1x1	-40.930301	0.003818
2x2x2	-40.926609	0.000126
3x3x3	-40.926432	0.000051
4x4x4	-40.92645	0.000033
5x5x5	-40.926463	0.000020
6x6x6	-40.926471	0.000012
8x8x8	-40.926478	0.000005
10x10x10	-40.92648	0.000003
12x12x12	-40.926481	0.000002
16x16x16	-40.926482	0.000001
20x20x20	-40.926483	0.000000
24x24x24	-40.926483	Converged
32x32x32	-40.926483	Converged
48x48x48	-40.926483	Converged
64x64x64	-40.926483	Converged

Figure 5

The Thermal Expansion of MgO

Table 3

Temperature / K	300	800	1300	1800
Thermal Expansion Coefficient for the Quasi-harmonic Approximation x10-6 / K-1	26.7	32.7	38.5	47.9
Thermal Expansion Coefficient for the Molecular Dynamics Method x10-6 / K-1	30.1	32.6	23.8	39.7
Literature Values for the Thermal Expansion Coefficient x10-6 / K-1	31.2	42.6	47.1	51.3

Figure 6

Figure 7

Figure 8

Figure 9