Rep:Comp:jtl14y3t1

The Free Energy and Thermal Expansion of MgO

Introduction

Through the operating system RedHat Linux, DLVisualize can graphically model MgO and the behaviour of MgO can be studied when it is heated. MgO's free energy and thermal expansion coefficient can be calculated through both quasi-harmonic approximation and molecular dynamics calculations with the aid of a modelling code GULP.

It is ideal to run calculations on MgO as it is a simple face-centred cubic crystal. Quasi-harmonic approximation calculations are run on its primitive cell: a = b = c = 2.99 Å, α = β = γ = 60°, N_p = 2, V_p = 37.36 ų; and molecular dynamic calculations are run on a supercell: a = b = c= 8.42 Å, α = β = γ = 90°, N_p = 64, V_p = 597.76 ų.

Reciprocal space is the Fourier transform of direct space. For a crystal in reciprocal space, energy can be plotted as a function of k; as a crystal can be repeated infinitely, there will an infinite number of vibrational modes and thus an infinite amount of k-points. From equation 1., it can be seen that if a crystal has a large lattice constant,$a$ , then its value in reciprocal space $a^{*}$ will be small.

$$\begin{gathered} a^{*}=\frac{\mathrm{2}\mathrm{\pi }}{a} \\ (1) \end{gathered}$$

The Helmholtz free energies that are calculated show numerically the density of states. The well known equation below shows that the free energy can be minimised when temperatures is increased as a function of

$V$

.

$$\begin{gathered} A=U-TS \\ (2) \end{gathered}$$

The quasi-harmonic approximation is used as the equilibrium bond length remains constant in a diatomic molecule with an exactly harmonic potential; however the quasi-harmonic approximation can be applied to materials in their solid state and takes into account that the lattice constant can increase from fluctuations in bond lengths(due to increased internal energy) in the structure; the harmonic approximation can then be applied to the larger lattice constant values, giving more accurate calculations.

Thermal expansion arises from an larger vibrations along harmonic oscillator due to the increased internal energy in the system, and thus the lattice constant and the equilibrium bond length increases, hence an expansion in 3 dimensions occurs. The thermal expansion coefficient, α, can be determined using the equation below:

$$\begin{gathered} \alpha =\frac{\mathrm{1}}{{\mathrm{V}}_0}{\left(\frac{\mathrm{\partial }\mathrm{V}}{\mathrm{\partial }T}\right)}_P \\ (3) \end{gathered}$$

• $\alpha$ is the thermal expansion coefficient
• $V$ is the cell volume per unit formula
• $T$ is the temperature in kelvin

Method

Figure 1. MgO visualised by DLVis

DLVisualize was launched from RedHat Linux, and MgO was graphically modelled (Fig. 1).

The internal energy per unit cell of MgO crystal was determined using a single point calculation, where the force field between Mg(+2e) and O (-2e) can be represented through a Buckingham potential. An output file from the GULP execution provides a break down of the total energy.

The phonon modes were computed for 50 points along the conventional path in k-space and graphed as a dispersion curve (Fig 2). The phonon mode at around 117 was animated to aid in the understanding of plotting the density of states graph.

Initially the DOS (density of states) was calculated for a single point in k-space (1x1x1) which consisted of four peaks, the grid size was then increased so that more k-points could be sampled and a smooth curve for the DOS could be graphed.

The suitable k-space grid was determined by increasing the grid size of the k-points, and identifying where the free energy converges at an appropriate value. The free energy of MgO was computed with respect to the quasi-harmonic approximation, for increasing grid sizes at 300 K and 0 GPa starting from an initial 1x1x1 up to a 32x32x32 grid, whereby the free energy converges.

With the grid size suitable at 16x16x16, the free energy was minimised and lattice constant was calculated for every 100 K between 0 K and 1000 K by optimising the structure of MgO at these temperatures.

MgO with 32 Mg-O units was loaded into DLVis and molecular dynamic simulations were ran.The number of particles in the simulation and the external pressure were kept constant. The molecular dynamic simulation was used to determine a final equilibrium volume for several different temperatures.

Results

Figure 2. Density of states (left) aligned with phonon dispersion (right)

The relationship between the density of states and the dispersion curve can be seen in Figure 2. The peaks from the DOS plotted for a single k-point align with the high symmetry point L from the dispersion curve.

Figure 3. Graph of free energy against temperature

The optimal grid size for DOS computations was determined to be 16x16x16 (Fig. 6), it provides a reasonable value for the free energy as it converges. Also, smaller grid sizes did not provide a smooth enough curve, and though larger grids are able to compute accurately (to the nearest 0.1 eV) they take longer to compute. A 16x16x16 grid size would also be suitable for GULP calculations on similar structures to MgO (eg. CaO). However, structures which are much larger than MgO (eg. Zeolites) would not require as large a grid, as they have a larger lattice parameter and thus a smaller reciprocal space needed to be sampled. On the other hand a metal (eg. Lithium) has a much large reciprocal space, and requires more k-points for adequate convergence. The free energy was plotted against temperature (Fig. 3) and it can be seen that when the temperature is increased the Helmholtz free energy can be minimised.This is because, as the temperature is increased, the entropy of the system also increases, and with Equation 2, the Helmholtz free energy is minimised. Unfortunately, as the temperature nears the melting point of MgO, 2825 °C/3098.15 K , the phonon modes poorly represent the actual motions of the ions as the solid structure of the lattice behaves less like a solid and more like a liquid, and it cannot be accurately modelled within the quasi-harmonic approximation.

Figure 4. Graph of lattice constant against temperature

The lattice constant was also plotted against temperature to examine how the volume of the cell changes with respect to temperature. The graph was fitted with a second-degree polynomial and shows that the volume per unit cell increases with increasing temperature at a relatively constant rate,

The thermal volume expansion coefficient of MgO was calculated with the aid of a volume/temperature graph (Fig. 5), the fitted second-degree polynomial [y = 1x 10^-8 T² + 1 x10^-5 T + 1] was differentiated to [2 x10^-8 T + 1 x10^-5] at T=300 K and with the use of Equation 3, was found to be 1.6x10^-5 K^-1 at 300 K. Literature values were often reported as thermal linear expansion which ranged from (1.0-1.4) x10^-5 K^-1 at 273-300 K;which is equivalent to a thermal volume coefficient of (2.0-2.8) x10^-5 K^-1, around double the calculated value.

Figure 5. Graph of normalised cell volume against temperature. LD = quasi-harmonic approximation, MD = molecular dynamics

A graph of normalised cell volume per formula unit vs temperature was plotted for both the quasi-harmonic approximation and the molecular dynamics (Fig. 5). It can be seen that points from the quasi-harmonic approximation follow closely a second-degree polynomial, whereas the molecular dynamics plot follows a linear trend. Molecular dynamics gives the thermal volume expansion coefficient of MgO to be 3x10^-5 K^-1, which is fairly close to the literature values . This is because MD can simulate random motion and allows for the system to develop and equilibriate following Newton's equations of motions.

As the temperature increases, the difference between the predicted thermal expansion coefficients reduces and the value computed within quasi-harmonic approximation tends to that predicted by molecular dynamics. As the quasi-harmonic approximation is fitted with a second-degree polynomial, a higher temperature will also lead to a higher predicted thermal expansion coefficient which becomes more accurate as anharmonic effects are taken into account.

Conclusion

The free energy was successfully minimised by running calculations at higher temperatures. At high temperatures, the entropic contribution in the thermodynamic equation 2 outweighs the increase in internal energy from the increased temperature.

The molecular dynamics calculations predicted the thermal expansion coefficient of MgO more accurately in the temperature range 0 K to 1000 K, however as the thermal expansion coefficient is a function of temperature, an extrapolation of the molecular dynamics plot would only yield a straight line; implying constant rate of expansion.

Additional Information

Table 1. Data from the computed free energies at different temperatures

Grid size	Helmholtz / eV
2	-40.9266
4	-40.9265
8	-40.9265
16	-40.9281
32	-41.3198

Figure 6. 16x16x16 density of states

Table 2. Table showing the convergence of energy as grid size increases

Temperature / K	Hemholtz / eV	Temperature / K	Hemholtz / eV
100	-40.9024	600	-41.0493
200	-40.9094	700	-41.1071
300	-40.9265	800	-41.1719
400	-40.9586	900	-41.243
500	-40.9994	1000	-41.3198

Table 3. Data from the calculated lattice constants at different temperatures:

Temperature / K	α/α_o	Temperature / K	α/α_o
100	1.008473	600	1.018802
200	1.009433	700	1.021684
300	1.011244	800	1.024696
400	1.013518	900	1.027818
500	1.016066	1000	1.031052

Table 4. Data of cell volume and temperature from molecular dynamic calculations

Temperature /K	V/Vo	Temperature / K	V/Vo
100	1.00046	550	1.013306
150	1.001086	650	1.016016
200	1.002291	700	1.018172
250	1.003908	750	1.018399
300	1.006161	800	1.022175
350	1.007372	850	1.023244
400	1.008313	900	1.024081
450	1.010385	950	1.023244
500	1.012382	1000	1.027974