Rep:Mod:SG94

Introduction

The software used is DLVisualise on Linux (RedHat). DLV is a graphical interface allowing the construction and studying of crystals. DLV features an interface to to molecular modelling code GULP (General Utility Lattice Program), which used in this experiment. This is a program which uses lattice dynamics to perform various types of simulations on 3D-periodic solids.

The energies and vibrations in a magnesium oxide (MgO) lattice are calculated to compute the free energy of the crystal and to observe the effects of temperature on the material (i.e. thermal expansion).

Two methods were used to observe the thermal expansion of the lattice: the quasi-harmonic (QH) model; and using molecular dynamics (MD). The QH approximation is a simple extension from the harmonic approximation which allows some anharmonic behaviour where the movement of atoms is in one direction. The MD techniques conducts the simulations by reproducing the actual; vibrations of the molecules and thus allows motion in three dimensions.

Internal Energy of MgO crystal

Fig1. Primitive cell of MgO

GULP was used to define the MgO crystal structure using a primitive unit cell (Figure 1). The cell vectors of this cell are (‎Å) and the LogFile can be found here:

0.00000	2.10597	2.10597
2.10597	0.00000	2.10597
2.10597	2.10597	0.00000

The Log File shows that the primitive cell has sides of length 2.9783 ‎Å and internal angles of 60 degrees. The total lattice energy is calculated by taking short-range repulsions and electrostatic potentials into account and summing over the lattice and represents the energy required to separate the ions of the lattice to infinite separation. The energy of the MgO lattice is -41.07531759 eV and corresponds to the energy of the lattice at 0 K (with no anharmonic contributions).

Phonon Modes of MgO

The lattice vibration modes (phonon modes) of MgO can be represented by a dispersion curve or density of states (DOS).

Phonon Dispersion Curves

Fig2. Phonon Dispersion of MgO along 50 points

In an infinite lattice, there are an infinite number of vibrations which can be simplified by considering the vibrations in the primitive unit cell of the crystal. Every possible vibration of the crystal has a k-label (symmetry points: $W,L,G,W,X,K$) which is related to its direction and wavelength. By summing over k, the total number of vibrations can be calculated. This is achieved by using GULP to calculate the phonon dispersion and by computing the phonons at 50 points along the path W-L-G-W-X-K in k-space.

The phonon dispersion curve (Figure 2) shows how the frequency of vibration changes with a change in k and the computed curve agrees with that found in literature .

As there is more than one atom in the primitive unit cell of MgO (N = 8*1/8 + 1 = 2), two types of phonons are displayed in the dispersion relation: acoustic and optical, corresponding to in-phase and out-of-phase movements respectively. As the calculations are taking three directions into account (x, y and z), there are a total of six phonon modes. This is explained by G. E. Peckham , where the three acoustic modes consist of one longitudinal and two transverse modes (which pass through the origin and represents the crystal behaving as a continuum). The three optical modes (3N-3) have non-zero frequencies and zero gradient at k=0. They are optical modes as they can interact with IR radiation.

By summing all of the possible k-points of the lattice, the free energy of the lattice can be computed.

Phonon Density of States (DOS)

The density of states summarises the dispersion curves by calculating the average number of vibrational modes at each frequency. Shrinking factors define the grid of k-points that will be used to compute the DOS curve and their effects on the curve were investigated.

A 1x1x1 grid computes the DOS from a single point and is shown in Table 1. It can be seen that there are four distinct peaks which corresponds to the k-point L. There are only four peaks present in the DOS rather than 6 as bands become degenerate at frequencies at 288.49 and 351.76 cm^-1. This explains why the intensities of these two peaks are double that of the other two at higher frequencies (676.23 and 818.82 cm^-1).

Fig3. Density of States at grid size 1x1x1

Fig4. Density of States at grid size 14x14x14

Fig5. Density of States at grid size 15x15x15

Fig6. Density of States at grid size 30x30x30

By increasing the grid size, more points in k-space are taken to compute a more representative DOS and the graph starts to form a curve at 16x16x16. Several grid sizes were investigated to find a compromise between accuracy of results and CPU time (below). The minimum grid size for a good approximation is the 15x15x15 grid. Below this size, there are cuts on the x-axis (not a curve) as the number of sampled k-points were insufficient; and above this grid size, the DOS becomes smoother and more representative with an increase in the number of sampled possible vibrations. It can be seen that there is very little difference in the DOS of the two largest grid sizes. This occurs when a sufficient number of points in k-space are taken at a given size; and a further increase in would make little/no difference. The CPU time for the 40x40x40 grid size was noticeably longer than the 30x30x30 with very little difference in the DOS and so it was decided that the latter was the optimum size used for future calculations.

The phonon dispersion curve states the possible frequencies of phonons at each k-label and by differentiating this curve, the DOS produced gives information about the entire set of phonon frequencies. The frequencies which more k-labels correspond to increases the density of states at that frequency.

Computing the Free Energy

Table 1. The Effects of Grid Size on the Helmholtz Free Energy of MgO

Grid size	Helmholtz Free Energy / eV	Accuracy
1x1x1	- 40.9303	0.1 eV
2x2x2	- 40.9266	0.1 eV
3x3x3	- 40.9264	1 meV
4x4x4	- 40.9265	0.5 meV
10x10x10	- 40.9265	0.1 meV
20x20x20	- 40.9265	0.1 meV
30x30x30	- 40.9265	n/a

The quasi-harmonic approximation was used to compute the free energy of MgO by summing all of the possible vibrational modes of the crystal at 300 K. The free energies computed at various sizes using GULP are displayed in Table 2. An increase in grid size lowers the free energy of the lattice by increasing the degree of accuracy. This is the case due to the increase in number of microstates sampled ($S=klnW$), increasing the entropic term and in turn, lowering the Gibbs free energy ($G=H-TS$).

Table 2 shows that 30x30x30 grid size is an appropriate choice for accurate results with an appropriate CPU time and found the free energy converge to - 40.9265 eV. These grid sizes are applicable to other similar crystal structures such as calcium oxide (CaO), where the lattice constant (4.800 Å) is similar to that of MgO (4.211 Å).

Running a calculation on structures with a significant difference in lattice constants would require a different sized grid for the same level of accuracy. For zeolite structures such as faujasite with a larger lattice constant (between 24.2-25.1 Å ), a smaller grid size (fewer sampled k-points) is required to achieve sufficiently accurate results.

Metal structures such as lithium have a smaller lattice constant of 3.49 Å , consisting of a smaller reciprocal space and so it is predicted that more k-points (smaller grid size) are required for the same level of accuracy. However, the lithium lattice consists of a regular lattice of positively charged ions in a delocalised sea of electrons which experiences less coulombic attraction compared to the ions in the MgO crystal. This higher symmetrical nature of the lithium lattice results in fewer symmetry points and so a fewer number of points need to be computed along the path in k-space.

Thermal Expansion of MgO

When a substance is heated, the energy is converted to kinetic energy in its molecules, resulting in the expansion of the material. The linear thermal expansion coefficient (CTE) indicates the degree of expansion in one dimension (i.e. the lattice constant) of a material upon heating and the equation used to calculate this value is given by:

${\alpha }_L=\frac{1}{L_0}\frac{dL}{dT}$

Where $L_0$ is the lattice constant at 0 K and $dL/dT$ is the rate of change of the lattice constant per unit change in temperature.

The volumetric CTE at constant pressure is given by:

${\alpha }_V=\frac{1}{V}{\left(\frac{\partial V}{\partial T}\right)}_p$

This equation is important for a gas as the volume will vary with a change in pressure but this is not important for MgO as it is a solid.

The crystal structure of MgO means that it has identical properties in all three directions, i.e. it is isotropic. This means that the volumetric CTE should be three times that of the linear CTE: ${\alpha }_V=3{\alpha }_L$

The CTEs of MgO were predicted using the the quasi-harmonic approximation and molecular dynamics.

Quasi-Harmonic Approximation

Graph 1. Effect of temperature on the Gibbs free energies of MgO

Graph 2. Effect of temperature on the lattice constants of MgO

This method optimises the MgO structure by minimising the Gibbs free energy through computing the internal energy and phonons at a sequence of geometries. By varying the temperature from 0 K to 1000 K in 100 K intervals, the effects of temperature on the Gibbs free energy and lattice constants are discovered (Graphs 1 and 2).

The negative correlation between Gibbs Energy and temperature can be explained through the definition of Gibbs where temperature is negatively proportional to the free energy; while the positive correlation between the lattice constant (bond length) and temperature proves the thermal expansion of MgO upon an increase in temperature.

$G=H-TS$

Where G is the Gibbs free energy, H is the enthalpy (J mol^-1), S is the entropy (J K^-1 mol^-1) and T is the temperature in Kelvin (K).

The volumetric expansion coefficient from 0 - 1000 K was calculated to be 2.2485x10^-5 K^-1 while the linear expansion coefficient for the same temperature range is 7.4400x10^-6 K^-1. The volumetric CTE is approximately three times that of the linear CTE and so the isotropic property is confirmed. From 0 - 300 K, the thermal CTE is 9.47x10^-6 K^-1 which is in general agreement with the literature value of 10.5x10^-6 K^-1.

This method explains the thermal expansion based on the harmonic approximation while including anharmonic contributions to an extent. This allows the harmonic approximation to hold for every value of the lattice constant (i.e. phonon frequencies are volume-dependent), and explains why the bond length increases with temperature (as volume increases). The quasi-harmonic model in reciprocal space reduces the size of the force constant matrix by approximating that the vibration of the atoms in the primitive unit cell has the same magnitude and direction and only differs in phase. This simplifies the model and reduces the effort required to calculate the free energy and other thermal properties of the MgO crystal. .

The relationship between the Gibbs free energy and temperature should be a linear one based on the derivation below:

$G=H-TS$

Which is the same as:

$G=U+pV-TS$

Where U is the internal energy; p is the pressure; and V, the volume.

As U is defined as:

$U=\delta U+\delta W=TdS-PdV$

The Gibbs energy becomes:

$G=Vdp-SdT$

So at constant pressure as stated in the QH model, the values should be a linear fit with gradient -S.

${\left(\frac{dG}{dT}\right)}_p=-S$

This trend however, is not reflected in Graph 2 and this can be justified through an over-simplification of the QH model, where there is an increase in anharmonic behaviour at increasingly high temperatures.

As the temperature approaches the melting point, the model breaks down as the harmonic approximation can no longer be applied at these lattice parameters. The QH model only takes linear movement into account and is therefore not applicable to molten MgO, where Mg²⁺ and O^2- behave as free ions.

In the harmonic potential, the displacement of the bond length from equilibrium is limited and so an increase in the lattice constant is observed up to the point where a further increase in temperature will not lead to bond cleavage (a larger lattice constant). The QH model on the other hand, includes an anharmonic factor, where the dissociation of MgO is possible and the lattice constant can approach infinity at increasing temperatures.

Molecular Dynamics

Fig7. 32 unit MgO structure used for MD simluations

The molecular dynamics method studies the physical movements of atoms and molecules where atomic and molecular interactions are allowed. The system is studied and evolves over a period of time according to Newton's second law: $F=ma$. The timestep at which the vibrations to be sampled was taken to be 1 femtosecond (fs) and the number of equilbration steps and production steps were set to 500. This situation obeys the ergodic hypothesis where when the number of sampled configurations are sufficiently large, the average calculated is representative for the time averaged system. Originally due to L. Boltzmann. Increasing the size of the crystal used in the MD calculations rather than a single primitive cell allows more flexibility when the system is vibrating. By using a structure containing 32 units (64 atoms) movements of atoms that are not in phase are allowed.

In both cases, the anharmonic effects become significant at high temperatures and Graph 3 shows that the cell volumes obtained from the QH model is in good agreement with the MD data. The values obtained from the QH model are generally higher than those from the MD simulations while the latter shows a stronger anharmonic effect at higher temperatures shown by the cross-over in cell volume around 930 K. The strong linear fit between the cell volume and the temperature using this method further suggests that the calculations are based on a less harmonic approximation. The QH model describes the thermodynamic properties of the crystal lattice using the Grüneisen parameter (γ), describing the effect of a change in temperature on the volume of the lattice, and is a measure of the anharmonicity of the system. This method measures the properties using a primitive cell, which is comparable of the larger system used in the MD simulations. A larger size allows movements that are not in phase and an increase in temperature and hence volume would result in more random movement, and so the values obtained from this value differ from those using the QH model.

Fig8. MgO structure post-simulations

Graph 3. Comparing the Quasi-Harmonic and Molecular Dynamics methods and their effects of temperature on cell volume

The thermal CTE calculated using the molecular dynamics method is 2.9993x10^-5 K^-1 between 100 - 1000 K, similar to the value computed from the QH model of 2.4984x10^-5 K^-1, and the value from 0 - 300 K (3.0280x10^-5 K^-1) agrees well with the literature value of 1.05x10^-5 K^-1. The deviation from the literature value can be justified using the size of the structure used for the simluations. The 32 unit stucture is a 2x2x2 lattice of the conventional unit cell and is noticeably smaller than the size required for an accurate Gibbs free energy result which suggests that a much larger size is required for an accurate result using the MD method. As the MD method allows movement in 3-dimensions, the thermal expansion results in a more irregular structure in comparison to the structures obtained from the QH model, which only allows movement in 1-direction. This leads to the assumption that at high temperatures, the cell volumes computed using the MD method will be drastically higher than that approximated using the QH model.

Conclusion

The thermal expansion of crystal MgO was studied using the QH model and MD methods. These were both found to be relatively accurate at low temperatures but found deviations at higher temperatures due to the larger anharmonic effects which are not accounted for fully by the two methods. The approximations used allowed a more time-efficient determination of the effects of temperature on various thermodynamic properties such as the Gibbs free energies and cell volumes.

Calculations

Quasi-harmonic approximation thermal coefficient from 0-300K

Linear CTE:

${\alpha }_L=\frac{1}{L}{\left(\frac{dL}{dT}\right)}_p$

${\alpha }_L=\frac{1}{2.986563}{\left(\frac{2.989390-2.986563}{300}\right)}_p$

${\alpha }_L=3.152\times 10^{-6}{K}^{-1}$

Volumetric CTE:

${\alpha }_V=\frac{1}{V}{\left(\frac{\partial V}{\partial T}\right)}_p$

${\alpha }_V=\frac{1}{18.836496}{\left(\frac{18.890024-18.836496}{300}\right)}_p$

${\alpha }_V=9.472\times 10^{-6}{K}^{-1}$

Quasi-harmonic approximation thermal coefficient from 0-1000K

Linear CTE:

${\alpha }_L=\frac{1}{L}{\left(\frac{dL}{dT}\right)}_p$

${\alpha }_L=\frac{1}{2.986563}{\left(\frac{3.008783-2.986463}{1000}\right)}_p$

${\alpha }_L=7.440\times 10^{-6}{K}^{-1}$

Volumetric CTE:

${\alpha }_V=\frac{1}{V}{\left(\frac{\partial V}{\partial T}\right)}_p$

${\alpha }_V=\frac{1}{18.836496}{\left(\frac{19.260042-18.836496}{1000}\right)}_p$

${\alpha }_V=2.249\times 10^{-5}{K}^{-1}$

Molecular dynamics thermal coefficient from 100-1000K

${\alpha }_V=\frac{1}{V}{\left(\frac{\partial V}{\partial T}\right)}_p$

${\alpha }_V=\frac{1}{18.732935}{\left(\frac{(19.238604-18.732935)}{1000}\right)}_p$

${\alpha }_V=2.999\times 10^{-5}{K}^{-1}$