=Thermal Expansion of MgO=

==Introduction==

This investigation involves modelling the thermal expansion of a crystal of MgO, and calculating its thermal expansion coefficient <math>\alpha </math>:

<math>\alpha = \frac{1}{V_0}\left(\frac{\partial V}{\partial T}\right)_P </math> '''Eq. 1'''

The Quasi-Harmonic Approximation is used to calculate the vibrational energy levels and to understand the phonon dispersion of the material and its vibrational density of states. This approach is compared to the Molecular dynamics results used to simulate the vibrations as random motions of atoms within a cell.

To calculate the thermal expansion coefficient the lattice vibrations must be calculated; the main contributor to thermal expansion. To calculate how these vibrations change with temperature and therefore how the separation between the atoms change, giving rise to thermal expansion, the harmonic potential model is insufficient. This is due to the equilibrium distance between atoms being invariant within this model as the displacement of atoms is symmetrical about their equilibrium positions. Instead the quasi-harmonic approximation (QHA), which is a phonon-based model, is used to model the thermal expansion by assuming the frequencies are volume dependent: <math>\omega_j\left(V\right)</math>

In this model the free energy is:

<math>F=</math> (free energy from vibrational theory)

The k within the equations above is the wave vector, essentially a label that can be given to each unique point within the unit cell which corresponds to a unique vibrational energy within the lattice. To find the k-values the crystal must be converted into reciprocal space with lattice parameters <math>a^*</math> found using the real space lattice parameter <math>a</math>:

<math> a^*=\frac{2\pi}{a}</math>

Within the reciprocal cell <math>k</math> is periodic between <math>-\frac{\pi}{a}</math> and <math>\frac{\pi}{a}</math>.

The frequencies of the vibrations <math>\omega</math> vary with <math>k</math> by the following relationship:

<math>k=</math>

The crystal structure of MgO in a conventional cell is a face-centered cubic lattice which consists of a total of 8 atoms '''(Fig 1)''', it can also be modelled as a primitive cell '''(Fig. 2)'''; the simplest unit cell, with one ion in the centre and the opposing ion on each corner giving a total of 2 atoms within the lattice. The conventional cell has 4 times the volume of the primitive unit cell however is a more realistic representation of the symmetry of the crystal therefore is still used in many crystallographic studies. In this investigation the primitive cell will be dealt with initially when calculating using the QHA, but a super cell '''(Fig. 3)''' is used for the molecular dynamics simulations which is made up of 32 conventional cells.

{|style="margin: 0 auto;"
| [[File:Convencell.png|thumb|'''Fig. 1:'''Conventional cell of MgO|upright]]
| [[File:Primcell.png|thumb|'''Fig. 2:''' Primitive cell of MgO|upright]]
| [[File:Supercell.png|thumb|'''Fig. 3:''' Supercell of MgO]]
|}

==Method==

All calculations were carried out using DLVisualize, a graphical interface used for modelling. Within this programme the modelling code GULP (General Utility Lattice Program) was utilised to calculate the theoretical values evaluated in the Results and Discussion below.

The internal energy (<math>U</math>) of the primitive unit cell was initially calculated by modelling the crystal as a purely ionic material which was used in future calculations of the Helmholtz Free Energy (<math>A</math>).

To enable us to visualise how the frequencies of vibrations of this lattice vary with <math>k</math> a phonon dispersion curve was calculated using a sample of 50 k-points and represented graphically '''(FIg. 4)'''. Subsequently a calculation was run to find the density of states for a lattice grid size of 1x1x1. This calculation was then repeated, doubling the shrinking factor each time until a grid size of 64x64x64 was obtained.

'''Initially the thermal expansion coefficient was calculated using the QHA:'''

Using the optimum grid size, found to be 32x32x32 from the above simulations, the Helmholtz Free Energy was found for the MgO crystal at temperatures between 0 K and 1000 K in 100 K intervals. The lattice parameters and the total volume for each cell was recorded at each temperature and plotted [Fig ....]

From the plot of cell volume with respect to temperature, the gradient of the linear section was taken with the initial volume and used in '''Eq. 1''' to give a value of the thermal expansion coefficient, <math>\alpha</math>, as <math>2.64</math>x<math>10^{-5} K^{-1}</math> for a temperature range 400 - 1000 K.

'''The thermal expansion coefficient was then calculated using Molecular Dynamics.'''

Within this computation the system has a timestep of 1 femtosecond <math>(10^{-15} s)</math> and is run for 500 steps after equilibrium has been reached at the selected temperature. The same temperature range from 0 K to 1000 K was used as above. The system used was made up of 32 MgO conventional cells. The average equilibrium volumes were recorded at each temperature for this simulation [Fig ...]. '''Eq. 1''' was used to calculate the thermal expansion coefficient as above using the gradient of the volume as a function of temperature line of best fit and the initial volume found. Giving a value of <math>\alpha</math> as <math>2.86</math>x<math>10^{-5} K^{-1}</math> with a temperature range 100 - 1000 K.

==Results and Discussion==

The internal energy per primitive cell of MgO was found to be -41.0753 eV. MgO was modelled as purely ionic crystal with Mg 2+ ions and O 2- ions having simple Buckingham repulsive potentials operating between the ions. <ref>T.S.Bush, J.D.Gale, C.R.A.Catlow and P.D. Battle J. Mater Chem., 4, 831-837 (1994) </ref>, this energy is also described as the 'binding energy' of the crystal as it is the energy necessary to separate the ions to an infinite distance.

The Phonon dispersion curve is shown in Fig. 4 with k-space plotted against frequency of the vibrations.

[[File:Phonon_Dispersion.png|500x500px|thumb|Fig. 4: Phonon dispersion curve for MgO|centre]]

The labels on the x-axis are points of interest in k-space. For example the origin is labelled <math>\Gamma</math> with lattice coordinates [0, 0, 0]. The bands in the phonon dispersion model correspond to...

To compute the free energy using the QHA, the sum over all vibrations of the lattice is required. To ensure the correct level of accuracy a Density of States (DOS) calculation was run with increasing shrinking factors. The grid sizes used were: 1x1x1 '''(Fig. 5)''', 2x2x2 '''(Fig. 6)''', 4x4x4 '''(Fig. 7)''', 8x8x8 '''(Fig. 8)''', 16x16x16 '''(Fig. 9)''', 32x32x32 '''(Fig. 10)''' and 64x64x64 '''(Fig. 11)'''.

[[File:DOS1ce.png|400x400px|thumb|'''Fig. 5:''' Density of states using 1x1x1 grid size|centre]]

In the initial DOS calculation for the 1x1x1 grid size the graph was computed for a single k-point, corresponding to the point labelled L in k-space with coordinates (0.5, 0.5, 0.5). This can be determined from the examination of the four frequency signals in the 1x1x1 DOS graph '''(Fig. 5)''' corresponding to the frequencies that are observed at k-point L on the phonon dispersion graph '''(Fig. 4)''', and also read off from the Log file.

<div><ul>
<li style="display: inline-block;"> [[File:DOS2Ce.png|400px|thumb|'''Fig. 6:''' 2x2x2 ]]</li>
<li style="display: inline-block;"> [[File:DOS4Ce.png|400px|thumb|'''Fig. 7:''' 4x4x4]]</li>
<li style="display: inline-block;"> [[File:DOS8.png|400px|thumb|'''Fig. 8:''' 8x8x8]] </li>
<li style="display: inline-block;"> [[File:DOS16Ce.png|400px|thumb|'''Fig. 9:''' 16x16x16]] </li>
<li style="display: inline-block;"> [[File:DOS32Ce.png|400px|thumb|'''Fig. 10:''' 32x32x32]] </li>
<li style="display: inline-block;"> [[File:DOS64.png|400px|thumb|'''Fig. 11:''' 64x64x64]] </li>
</ul></div>

It is clear from the above figures that as the shrinking factor is increasing the accuracy in the Density of States graphs in also increasing. This is due to the number of k-points sampled increasing in each grid size, therefore more phonon frequencies are recorded which increases the smoothness of the graphs. However using this computation method the difference between using a shrinking factor of 32 and a shrinking factor of 64 has very little effect on the shape of the curves, this can be seen in Figure 10 and 11. Therefore the grid size of 32x32x32 was chosen for the calculation of the free-energy of MgO using the QHA as this size being the optimum balance of accuracy and CPU time in this case.

For similar crystal such as CaO, this grid size would also be appropriate as it's lattice parameters are very similar to MgO and they have the same unit cell structure. However for a Zeolite such as Faujasite a smaller grid size would be suitable. This is due to Faujasite having a much larger lattice constant, <math>a</math>: 24.638 Å <ref>Wise, W.S. (1982) New occurrence of faujasite in southeastern California. Am. Miner. 67, 794-798.</ref> compared to MgO with lattice parameter of 4.212 Å, therefore in reciprocal space the lattice parameter, <math>a^*</math> is much smaller and a smaller sample of k-points would give a sufficient degree of accuracy for the DOS calculation. When considering a metal such as Li, the environment of the electrons must be taken into account when determining the necessary grid size for the DOS calculation. In the vibrational structure of a metal, the vibrational bands are much less disperse when compared to the bands of an insulator. This can be explained by the behaviour of the electrons shielding the metal ions as they vibrate within the crystal to minimise the repulsive forces, therefore the electrons do not move as freely and are therefore have less freedom. Thus when calculating the DOS for a metal, a smaller grid size will be suitable as fewer k-points are needed.

An optimum grid size of 32x32x32 was further confirmed when computing the Free Energy of the primitive cell of MgO at 300 K shown in Table 1. As the value for the Helmholtz Free Energy can be said to be converging as grid size increases, the value for the shrinking factor of 64 can be said to be the most accurate. Therefore the Free Energy for the 2x2x2 grid size is accurate to 1 meV and 0.5 meV and the 4x4x4 grid is accurate to 0.1 meV which can be seen from the values in the final column in Table 1. However as a further accuracy can be obtained a larger grid size was preferred and 32x32x32 was chosen as it gave the same value for the Helmholtz Free Energy as 64x64x64 to 6 decimal places without taking as much CPU time.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|+ style="text-align: left;" |'''Table 1'''
!Grid Size
!Helmholtz Free Energy (eV)
!Difference in Energy from 64x64x64 grid (eV)
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|3.82x10<sup>-3</sup>
|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|1.26x10<sup>-4</sup>
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|3.30x10<sup>-5</sup>
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|5.00x10<sup>-6</sup>
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|1.00x10<sup>-6</sup>
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|
|}

Using the Quasi-harmonic approximation the Helmholtz Free Energy, lattice constant,<math>a</math>, and cell volume of a primitive lattice of MgO were calculated and recorder at 100 K intervals within the range 0 K to 1000 K '''(Table 2)'''.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|+ style="text-align: left;" |'''Table 2'''
!Temperature (K) !! Helmholtz Free Energy (eV) !! Lattice constant (Å) !! Cell volume (Å<sup>3</sup>)
|-
|0 || -40.901906 || 2.986563 || 18.838260
|-
|100 || -40.902420 || 2.986563 || 18.838260
|-
|200 || -40.909377 || 2.987604 || 18.856201
|-
|300 || -40.928125 || 2.989390 || 18.890025
|-
|400 || -40.958595 || 2.991629 || 18.932507
|-
|500 || -40.999436 || 2.994134 || 18.980110
|-
|600 || -41.049316 || 2.996820 || 19.031221
|-
|700 || -41.107120 || 2.999643 || 19.085057
|-
|800 || -41.171892 || 3.002587 || 19.141316
|-
|900 || -41.243018 || 3.005634 || 19.199638
|-
|1000 || -41.319849 || 3.008783 || 19.260042
|}

The value for the Helmholtz Free Energy is becoming more negative as the temperature is increased which is inline with the equation used by the QHA to calculate the Free Energy (F):

<math>F\left(T, V\right) = U\left(V\right) - TS\left(T, V\right)</math>

As temperature, <math>T</math>, increases the term <math>TS\left(T, V\right)</math> becomes more negative and dominates the equation causing the free energy, <math>F</math> to decrease'''(Fig. 12)'''. Conversely the lattice constant, <math>a</math>, increases with temperature '''(Fig. 13)''' as does the volume '''(Fig.14)'''. The volume of the cell increases with temperature as the higher the temperature the more kinetic energy which can manifest as vibrational energy within the crystal. As this energy increases the amplitude of the crystal vibrations is increasing causing the equilibrium bond length to become greater as the repulsive forces become stronger at shorter distances therefore the distance between the atoms becomes greater<ref>Ziman, J.M. (1972). Principles of the Theory of Solids, 2nd Ed. Cambridge University Press. ISBN 0521297338 9780521297332</ref>. This is the origin of thermal expansion within this approximation.

<div><ul>
<li style="display: inline-block;"> [[File:FreeEnergyQHA.png|500px|thumb|'''Fig. 12''']]</li>
<li style="display: inline-block;"> [[File:Latticeconst.png|500px|thumb|'''Fig. 13''']]</li>
</ul></div>

[[File:Thermalexpce.png|800px|thumb|'''Fig. 14'''|centre]]
A line of best fit was plotted for the linear region of the data (400 - 1000 K) and its gradient taken for the rate of change of volume with temperature in '''Eq. 1''' and <math>V_0</math> taken as the volume of the cell at 400 K.

<math>\alpha = \frac{1}{18.932507}0.0005</math>

<math>\alpha = 2.640960</math>x<math>10^{-5}</math> <math>K^{-1}</math> (400 - 1000 K)

The thermal expansion coefficient was found using molecular dynamics. This computational approach models the movements of the atoms in the lattice and calculates the thermodynamic properties from the time-averaged behaviour. A supercell made up of 32 conventional lattice cell is used due to the movement of the atoms no longer being periodical therefore more unit cells must be used to model this as accurately as possible. The larger the system in the computation the closer to a real system it becomes and the more accurate the values produced. However due to this computation taking much higher CPU time than when calculating in the QHA only a system of 32 conventional cells were used.

The volume change was plotted against temperature '''(Fig 15)'''.

[[File:MDTherm.png|800px|thumb|'''Fig. 15'''|centre]]

The volume of the system increases linearly, a value for 0 K cannot be calculated for this approximation as MD is based on classical mechanics and does not take into account quantum effects. Therefore at 0 K there would be no kinetic energy and the atom would be at rest; due to the uncertainty principle <ref>W. Heisenberg, Remarks on the origin of the relations of uncertainty, in: W. Price, S. Chissick (Eds.), The Uncertainty Principle and Foundations of Quantum Mechanics. A Fifty Years’ Survey, J. Wiley & Sons, London, 1977, pp. 3–6.</ref> both the position and momentum of a particle cannot be known simultaneously, and hence the computation breaks down here.

As with the QHA a line of best fit was plotted for the data and its gradient used for the rate of change of volume with temperature in '''Eq. 1''' and <math>V_0</math> the volume at 100 K.

<math>\alpha = \frac{1}{599.513295}0.0173</math>

<math>\alpha = 2.885674</math>x<math>10^{-5}</math> <math>K^{-1}</math> (100 - 1000 K)

==Conclusion==

Both the Quasi-harmonic approximation (QHA) and the Molecular Dynamics (MD) computations gave values for the thermal expansion coefficient: 2.640960x10<sup>-5</sup> K<sup>-1</sup> and 2.668831x10<sup>-5</sup> K<sup>-1</sup> respectively.

The results for the change in volume per unit cell with respect to temperature is compared for the QHA and for MD '''Fig. 16'''.

[[File:Comparevol.png|800px|thumb|'''Fig. 16'''|centre]]

When compared with the literature value of ... (teap range)<ref>email ref</ref> both values give relatively good estimations however the MD result is closer to the literature value.

When in lower temperature ranges the QHA can be said to be more accurate due to a value for 0 K being attainable as this model takes into account quantum effects so the zero point energy is used within the calculation for free energy. As the temperature increases the QHA breaks down, with its calculation for volume increasing more rapidly with temperature than in the MD case '''(Fig. 16)'''. The overestimation in volume for the QHA case is due to the neglect of higher-order anharmonic terms <ref> http://onlinelibrary.wiley.com/store/10.1029/94GL01370/asset/grl7676.pdf?v=1&t=iyyt51qe&s=5a4d486609f7da491c838c7fbdd84526802b5647 </ref>. As the temperature increases the volume in the QHA will just continue increasing infinitely as the model does not allow for the breaking of the bonds between the ions, therefore cannot be used near or above the melting point of a crystal. In higher energy ranges the MD approximation is preferred as it can take into account the breaking of bonds.

Overall, in the temperature ranges within these computations either method can be used as both give a relatively accurate result. It is difficult to say which method is preferred overall as the accuracy varies with the temperature range of the experiment; the QHA approximation working better at lower temperatures and the MD more accurate at higher temperatures.

==References==

Rep:Mod:CE1014MgO

2017-02-09T20:32:41Z

Ce1014:

=Thermal Expansion of MgO=

==Introduction==

This investigation involves modelling the thermal expansion of a crystal of MgO, and calculating its thermal expansion coefficient <math>\alpha </math>:

<math>\alpha = \frac{1}{V_0}\left(\frac{\partial V}{\partial T}\right)_P </math> '''Eq. 1'''

The Quasi-Harmonic Approximation is used to calculate the vibrational energy levels and to understand the phonon dispersion of the material and its vibrational density of states. This approach is compared to the Molecular dynamics results used to simulate the vibrations as random motions of atoms within a cell.

To calculate the thermal expansion coefficient the lattice vibrations must be calculated; the main contributor to thermal expansion. To calculate how these vibrations change with temperature and therefore how the separation between the atoms change, giving rise to thermal expansion, the harmonic potential model is insufficient. This is due to the equilibrium distance between atoms being invariant within this model as the displacement of atoms is symmetrical about their equilibrium positions. Instead the quasi-harmonic approximation (QHA), which is a phonon-based model, is used to model the thermal expansion by assuming the frequencies are volume dependent: <math>\omega_j\left(V\right)</math>

In this model the free energy is:

<math>F=</math> (free energy from vibrational theory)

The k within the equations above is the wave vector, essentially a label that can be given to each unique point within the unit cell which corresponds to a unique vibrational energy within the lattice. To find the k-values the crystal must be converted into reciprocal space with lattice parameters <math>a^*</math> found using the real space lattice parameter <math>a</math>:

<math> a^*=\frac{2\pi}{a}</math>

Within the reciprocal cell <math>k</math> is periodic between <math>-\frac{\pi}{a}</math> and <math>\frac{\pi}{a}</math>.

The frequencies of the vibrations <math>\omega</math> vary with <math>k</math> by the following relationship:

<math>k=</math>

The crystal structure of MgO in a conventional cell is a face-centered cubic lattice which consists of a total of 8 atoms '''(Fig 1)''', it can also be modelled as a primitive cell '''(Fig. 2)'''; the simplest unit cell, with one ion in the centre and the opposing ion on each corner giving a total of 2 atoms within the lattice. The conventional cell has 4 times the volume of the primitive unit cell however is a more realistic representation of the symmetry of the crystal therefore is still used in many crystallographic studies. In this investigation the primitive cell will be dealt with initially when calculating using the QHA, but a super cell '''(Fig. 3)''' is used for the molecular dynamics simulations which is made up of 32 conventional cells.

{|style="margin: 0 auto;"
| [[File:Convencell.png|thumb|'''Fig. 1:'''Conventional cell of MgO|upright]]
| [[File:Primcell.png|thumb|'''Fig. 2:''' Primitive cell of MgO|upright]]
| [[File:Supercell.png|thumb|'''Fig. 3:''' Supercell of MgO]]
|}

==Method==

All calculations were carried out using DLVisualize, a graphical interface used for modelling. Within this programme the modelling code GULP (General Utility Lattice Program) was utilised to calculate the theoretical values evaluated in the Results and Discussion below.

The internal energy (<math>U</math>) of the primitive unit cell was initially calculated by modelling the crystal as a purely ionic material which was used in future calculations of the Helmholtz Free Energy (<math>A</math>).

To enable us to visualise how the frequencies of vibrations of this lattice vary with <math>k</math> a phonon dispersion curve was calculated using a sample of 50 k-points and represented graphically '''(FIg. 4)'''. Subsequently a calculation was run to find the density of states for a lattice grid size of 1x1x1. This calculation was then repeated, doubling the shrinking factor each time until a grid size of 64x64x64 was obtained.

'''Initially the thermal expansion coefficient was calculated using the QHA:'''

Using the optimum grid size, found to be 32x32x32 from the above simulations, the Helmholtz Free Energy was found for the MgO crystal at temperatures between 0 K and 1000 K in 100 K intervals. The lattice parameters and the total volume for each cell was recorded at each temperature and plotted [Fig ....]

From the plot of cell volume with respect to temperature, the gradient of the linear section was taken with the initial volume and used in '''Eq. 1''' to give a value of the thermal expansion coefficient, <math>\alpha</math>, as <math>2.64</math>x<math>10^{-5} K^{-1}</math> for a temperature range 400 - 1000 K.

'''The thermal expansion coefficient was then calculated using Molecular Dynamics.'''

Within this computation the system has a timestep of 1 femtosecond <math>(10^{-15} s)</math> and is run for 500 steps after equilibrium has been reached at the selected temperature. The same temperature range from 0 K to 1000 K was used as above. The system used was made up of 32 MgO conventional cells. The average equilibrium volumes were recorded at each temperature for this simulation [Fig ...]. '''Eq. 1''' was used to calculate the thermal expansion coefficient as above using the gradient of the volume as a function of temperature line of best fit and the initial volume found. Giving a value of <math>\alpha</math> as <math>2.86</math>x<math>10^{-5} K^{-1}</math> with a temperature range 100 - 1000 K.

==Results and Discussion==

The internal energy per primitive cell of MgO was found to be -41.0753 eV. MgO was modelled as purely ionic crystal with Mg 2+ ions and O 2- ions having simple Buckingham repulsive potentials operating between the ions. <ref>T.S.Bush, J.D.Gale, C.R.A.Catlow and P.D. Battle J. Mater Chem., 4, 831-837 (1994) </ref>, this energy is also described as the 'binding energy' of the crystal as it is the energy necessary to separate the ions to an infinite distance.

The Phonon dispersion curve is shown in Fig. 4 with k-space plotted against frequency of the vibrations.

[[File:Phonon_Dispersion.png|500x500px|thumb|Fig. 4: Phonon dispersion curve for MgO|centre]]

The labels on the x-axis are points of interest in k-space. For example the origin is labelled <math>\Gamma</math> with lattice coordinates [0, 0, 0]. The bands in the phonon dispersion model correspond to...

To compute the free energy using the QHA, the sum over all vibrations of the lattice is required. To ensure the correct level of accuracy a Density of States (DOS) calculation was run with increasing shrinking factors. The grid sizes used were: 1x1x1 '''(Fig. 5)''', 2x2x2 '''(Fig. 6)''', 4x4x4 '''(Fig. 7)''', 8x8x8 '''(Fig. 8)''', 16x16x16 '''(Fig. 9)''', 32x32x32 '''(Fig. 10)''' and 64x64x64 '''(Fig. 11)'''.

[[File:DOS1ce.png|400x400px|thumb|'''Fig. 5:''' Density of states using 1x1x1 grid size|centre]]

In the initial DOS calculation for the 1x1x1 grid size the graph was computed for a single k-point, corresponding to the point labelled L in k-space with coordinates (0.5, 0.5, 0.5). This can be determined from the examination of the four frequency signals in the 1x1x1 DOS graph '''(Fig. 5)''' corresponding to the frequencies that are observed at k-point L on the phonon dispersion graph '''(Fig. 4)''', and also read off from the Log file.

<div><ul>
<li style="display: inline-block;"> [[File:DOS2Ce.png|400px|thumb|'''Fig. 6:''' 2x2x2 ]]</li>
<li style="display: inline-block;"> [[File:DOS4Ce.png|400px|thumb|'''Fig. 7:''' 4x4x4]]</li>
<li style="display: inline-block;"> [[File:DOS8.png|400px|thumb|'''Fig. 8:''' 8x8x8]] </li>
<li style="display: inline-block;"> [[File:DOS16Ce.png|400px|thumb|'''Fig. 9:''' 16x16x16]] </li>
<li style="display: inline-block;"> [[File:DOS32Ce.png|400px|thumb|'''Fig. 10:''' 32x32x32]] </li>
<li style="display: inline-block;"> [[File:DOS64.png|400px|thumb|'''Fig. 11:''' 64x64x64]] </li>
</ul></div>

It is clear from the above figures that as the shrinking factor is increasing the accuracy in the Density of States graphs in also increasing. This is due to the number of k-points sampled increasing in each grid size, therefore more phonon frequencies are recorded which increases the smoothness of the graphs. However using this computation method the difference between using a shrinking factor of 32 and a shrinking factor of 64 has very little effect on the shape of the curves, this can be seen in Figure 10 and 11. Therefore the grid size of 32x32x32 was chosen for the calculation of the free-energy of MgO using the QHA as this size being the optimum balance of accuracy and CPU time in this case.

For similar crystal such as CaO, this grid size would also be appropriate as it's lattice parameters are very similar to MgO and they have the same unit cell structure. However for a Zeolite such as Faujasite a smaller grid size would be suitable. This is due to Faujasite having a much larger lattice constant, <math>a</math>: 24.638 Å <ref>Wise, W.S. (1982) New occurrence of faujasite in southeastern California. Am. Miner. 67, 794-798.</ref> compared to MgO with lattice parameter of 4.212 Å, therefore in reciprocal space the lattice parameter, <math>a^*</math> is much smaller and a smaller sample of k-points would give a sufficient degree of accuracy for the DOS calculation. When considering a metal such as Li, the environment of the electrons must be taken into account when determining the necessary grid size for the DOS calculation. In the vibrational structure of a metal, the vibrational bands are much less disperse when compared to the bands of an insulator. This can be explained by the behaviour of the electrons shielding the metal ions as they vibrate within the crystal to minimise the repulsive forces, therefore the electrons do not move as freely and are therefore have less freedom. Thus when calculating the DOS for a metal, a smaller grid size will be suitable as fewer k-points are needed.

An optimum grid size of 32x32x32 was further confirmed when computing the Free Energy of the primitive cell of MgO at 300 K shown in Table 1. As the value for the Helmholtz Free Energy can be said to be converging as grid size increases, the value for the shrinking factor of 64 can be said to be the most accurate. Therefore the Free Energy for the 2x2x2 grid size is accurate to 1 meV and 0.5 meV and the 4x4x4 grid is accurate to 0.1 meV which can be seen from the values in the final column in Table 1. However as a further accuracy can be obtained a larger grid size was preferred and 32x32x32 was chosen as it gave the same value for the Helmholtz Free Energy as 64x64x64 to 6 decimal places without taking as much CPU time.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|+ style="text-align: left;" |'''Table 1'''
!Grid Size
!Helmholtz Free Energy (eV)
!Difference in Energy from 64x64x64 grid (eV)
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|3.82x10<sup>-3</sup>
|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|1.26x10<sup>-4</sup>
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|3.30x10<sup>-5</sup>
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|5.00x10<sup>-6</sup>
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|1.00x10<sup>-6</sup>
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|
|}

Using the Quasi-harmonic approximation the Helmholtz Free Energy, lattice constant,<math>a</math>, and cell volume of a primitive lattice of MgO were calculated and recorder at 100 K intervals within the range 0 K to 1000 K '''(Table 2)'''.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|+ style="text-align: left;" |'''Table 2'''
!Temperature (K) !! Helmholtz Free Energy (eV) !! Lattice constant (Å) !! Cell volume (Å<sup>3</sup>)
|-
|0 || -40.901906 || 2.986563 || 18.838260
|-
|100 || -40.902420 || 2.986563 || 18.838260
|-
|200 || -40.909377 || 2.987604 || 18.856201
|-
|300 || -40.928125 || 2.989390 || 18.890025
|-
|400 || -40.958595 || 2.991629 || 18.932507
|-
|500 || -40.999436 || 2.994134 || 18.980110
|-
|600 || -41.049316 || 2.996820 || 19.031221
|-
|700 || -41.107120 || 2.999643 || 19.085057
|-
|800 || -41.171892 || 3.002587 || 19.141316
|-
|900 || -41.243018 || 3.005634 || 19.199638
|-
|1000 || -41.319849 || 3.008783 || 19.260042
|}

The value for the Helmholtz Free Energy is becoming more negative as the temperature is increased which is inline with the equation used by the QHA to calculate the Free Energy (F):

<math>F\left(T, V\right) = U\left(V\right) - TS\left(T, V\right)</math>

As temperature, <math>T</math>, increases the term <math>TS\left(T, V\right)</math> becomes more negative and dominates the equation causing the free energy, <math>F</math> to decrease'''(Fig. 12)'''. Conversely the lattice constant, <math>a</math>, increases with temperature '''(Fig. 13)''' as does the volume '''(Fig.14)'''. The volume of the cell increases with temperature as the higher the temperature the more kinetic energy which can manifest as vibrational energy within the crystal. As this energy increases the amplitude of the crystal vibrations is increasing causing the equilibrium bond length to become greater as the repulsive forces become stronger at shorter distances therefore the distance between the atoms becomes greater<ref>Ziman, J.M. (1972). Principles of the Theory of Solids, 2nd Ed. Cambridge University Press. ISBN 0521297338 9780521297332</ref>. This is the origin of thermal expansion within this approximation.

<div><ul>
<li style="display: inline-block;"> [[File:FreeEnergyQHA.png|500px|thumb|'''Fig. 12''']]</li>
<li style="display: inline-block;"> [[File:Latticeconst.png|500px|thumb|'''Fig. 13''']]</li>
</ul></div>

[[File:Thermalexpce.png|800px|thumb|'''Fig. 14'''|centre]]
A line of best fit was plotted for the linear region of the data (400 - 1000 K) and its gradient taken for the rate of change of volume with temperature in '''Eq. 1''' and <math>V_0</math> taken as the volume of the cell at 400 K.

<math>\alpha = \frac{1}{18.932507}0.0005</math>

<math>\alpha = 2.640960</math>x<math>10^{-5}</math> <math>K^{-1}</math> (400 - 1000 K)

The thermal expansion coefficient was found using molecular dynamics. This computational approach models the movements of the atoms in the lattice and calculates the thermodynamic properties from the time-averaged behaviour. A supercell made up of 32 conventional lattice cell is used due to the movement of the atoms no longer being periodical therefore more unit cells must be used to model this as accurately as possible. The larger the system in the computation the closer to a real system it becomes and the more accurate the values produced. However due to this computation taking much higher CPU time than when calculating in the QHA only a system of 32 conventional cells were used.

The volume change was plotted against temperature '''(Fig 15)'''.

[[File:MDTherm.png|800px|thumb|'''Fig. 15'''|centre]]

The volume of the system increases linearly, a value for 0 K cannot be calculated for this approximation as MD is based on classical mechanics and does not take into account quantum effects. Therefore at 0 K there would be no kinetic energy and the atom would be at rest; due to the uncertainty principle <ref>W. Heisenberg, Remarks on the origin of the relations of uncertainty, in: W. Price, S. Chissick (Eds.), The Uncertainty Principle and Foundations of Quantum Mechanics. A Fifty Years’ Survey, J. Wiley & Sons, London, 1977, pp. 3–6.</ref> both the position and momentum of a particle cannot be known simultaneously, and hence the computation breaks down here.

As with the QHA a line of best fit was plotted for the data and its gradient used for the rate of change of volume with temperature in '''Eq. 1''' and <math>V_0</math> the volume at 100 K.

<math>\alpha = \frac{1}{599.513295}0.0173</math>

<math>\alpha = 2.885674</math>x<math>10^{-5}</math> <math>K^{-1}</math> (100 - 1000 K)

==Conclusion==

Both the Quasi-harmonic approximation (QHA) and the Molecular Dynamics (MD) computations gave values for the thermal expansion coefficient: 2.640960x10<sup>-5</sup> K<sup>-1</sup> and 2.668831x10<sup>-5</sup> K<sup>-1</sup> respectively.

The results for the change in volume per unit cell with respect to temperature is compared for the QHA and for MD '''Fig. 16'''.

[[File:Comparevol.png|800px|thumb|'''Fig. 16'''|centre]]

When compared with the literature value of ... (teap range)<ref>email ref</ref> both values give relatively good estimations but the MD approach gives a value closer to the literature. When in lower temperature ranges, however, the QHA can be said to be more accurate due to a value for 0 K being attainable as the zero point energy is taken into account in this model. As the temperature increases the QHA breaks down, with its calculation for volume increasing more rapidly with temperature than in the MD case '''(Fig. 16)'''. The overestimation in volume for the QHA case is due to the neglect of higher-order anharmonic terms <ref> http://onlinelibrary.wiley.com/store/10.1029/94GL01370/asset/grl7676.pdf?v=1&t=iyyt51qe&s=5a4d486609f7da491c838c7fbdd84526802b5647 </ref>. As the temperature increases the volume in the QHA will just continue increasing infinitely as the model does not allow for the breaking of the bonds between the ions therefore cannot be used near or above the melting point of a crystal. Therefore in higher energy ranges the MD approximation is preferred as it can take into account the breaking of bonds.

==References==

File:MDTherm.png

2017-02-09T20:07:46Z

Ce1014: Ce1014 uploaded a new version of File:MDTherm.png

File:MDTherm.png

2017-02-09T20:07:03Z

Ce1014: Ce1014 uploaded a new version of File:MDTherm.png

Rep:Mod:CE1014MgO

2017-02-09T20:06:43Z

Ce1014:

=Thermal Expansion of MgO=

==Introduction==

This investigation involves modelling the thermal expansion of a crystal of MgO, and calculating its thermal expansion coefficient <math>\alpha </math>:

<math>\alpha = \frac{1}{V_0}\left(\frac{\partial V}{\partial T}\right)_P </math> '''Eq. 1'''

The Quasi-Harmonic Approximation is used to calculate the vibrational energy levels and to understand the phonon dispersion of the material and its vibrational density of states. This approach is compared to the Molecular dynamics results used to simulate the vibrations as random motions of atoms within a cell.

To calculate the thermal expansion coefficient the lattice vibrations must be calculated; the main contributor to thermal expansion. To calculate how these vibrations change with temperature and therefore how the separation between the atoms change, giving rise to thermal expansion, the harmonic potential model is insufficient. This is due to the equilibrium distance between atoms being invariant within this model as the displacement of atoms is symmetrical about their equilibrium positions. Instead the quasi-harmonic approximation (QHA), which is a phonon-based model, is used to model the thermal expansion by assuming the frequencies are volume dependent: <math>\omega_j\left(V\right)</math>

In this model the free energy is:

<math>F=</math> (free energy from vibrational theory)

The k within the equations above is the wave vector, essentially a label that can be given to each unique point within the unit cell which corresponds to a unique vibrational energy within the lattice. To find the k-values the crystal must be converted into reciprocal space with lattice parameters <math>a^*</math> found using the real space lattice parameter <math>a</math>:

<math> a^*=\frac{2\pi}{a}</math>

Within the reciprocal cell <math>k</math> is periodic between <math>-\frac{\pi}{a}</math> and <math>\frac{\pi}{a}</math>.

The frequencies of the vibrations <math>\omega</math> vary with <math>k</math> by the following relationship:

<math>k=</math>

The crystal structure of MgO in a conventional cell is a face-centered cubic lattice which consists of a total of 8 atoms '''(Fig 1)''', it can also be modelled as a primitive cell '''(Fig. 2)'''; the simplest unit cell, with one ion in the centre and the opposing ion on each corner giving a total of 2 atoms within the lattice. The conventional cell has 4 times the volume of the primitive unit cell however is a more realistic representation of the symmetry of the crystal therefore is still used in many crystallographic studies. In this investigation the primitive cell will be dealt with initially when calculating using the QHA, but a super cell '''(Fig. 3)''' is used for the molecular dynamics simulations which is made up of 32 conventional cells.

{|style="margin: 0 auto;"
| [[File:Convencell.png|thumb|'''Fig. 1:'''Conventional cell of MgO|upright]]
| [[File:Primcell.png|thumb|'''Fig. 2:''' Primitive cell of MgO|upright]]
| [[File:Supercell.png|thumb|'''Fig. 3:''' Supercell of MgO]]
|}

==Method==

All calculations were carried out using DLVisualize, a graphical interface used for modelling. Within this programme the modelling code GULP (General Utility Lattice Program) was utilised to calculate the theoretical values evaluated in the Results and Discussion below.

The internal energy (<math>U</math>) of the primitive unit cell was initially calculated by modelling the crystal as a purely ionic material which was used in future calculations of the Helmholtz Free Energy (<math>A</math>).

To enable us to visualise how the frequencies of vibrations of this lattice vary with <math>k</math> a phonon dispersion curve was calculated using a sample of 50 k-points and represented graphically '''(FIg. 4)'''. Subsequently a calculation was run to find the density of states for a lattice grid size of 1x1x1. This calculation was then repeated, doubling the shrinking factor each time until a grid size of 64x64x64 was obtained.

'''Initially the thermal expansion coefficient was calculated using the QHA:'''

Using the optimum grid size, found to be 32x32x32 from the above simulations, the Helmholtz Free Energy was found for the MgO crystal at temperatures between 0 K and 1000 K in 100 K intervals. The lattice parameters and the total volume for each cell was recorded at each temperature and plotted [Fig ....]

From the plot of cell volume with respect to temperature, the gradient of the linear section was taken with the initial volume and used in '''Eq. 1''' to give a value of the thermal expansion coefficient, <math>\alpha</math>, as <math>2.64</math>x<math>10^{-5} K^{-1}</math> for a temperature range 400 - 1000 K.

'''The thermal expansion coefficient was then calculated using Molecular Dynamics.'''

Within this computation the system has a timestep of 1 femtosecond <math>(10^{-15} s)</math> and is run for 500 steps after equilibrium has been reached at the selected temperature. The same temperature range from 0 K to 1000 K was used as above. The system used was made up of 32 MgO conventional cells. The average equilibrium volumes were recorded at each temperature for this simulation [Fig ...]. '''Eq. 1''' was used to calculate the thermal expansion coefficient as above using the gradient of the volume as a function of temperature line of best fit and the initial volume found. Giving a value of <math>\alpha</math> as <math>2.86</math>x<math>10^{-5} K^{-1}</math> with a temperature range 100 - 1000 K.

==Results and Discussion==

The internal energy per primitive cell of MgO was found to be -41.0753 eV. MgO was modelled as purely ionic crystal with Mg 2+ ions and O 2- ions having simple Buckingham repulsive potentials operating between the ions. <ref>T.S.Bush, J.D.Gale, C.R.A.Catlow and P.D. Battle J. Mater Chem., 4, 831-837 (1994) </ref>, this energy is also described as the 'binding energy' of the crystal as it is the energy necessary to separate the ions to an infinite distance.

The Phonon dispersion curve is shown in Fig. 4 with k-space plotted against frequency of the vibrations.

[[File:Phonon_Dispersion.png|500x500px|thumb|Fig. 4: Phonon dispersion curve for MgO|centre]]

The labels on the x-axis are points of interest in k-space. For example the origin is labelled <math>\Gamma</math> with lattice coordinates [0, 0, 0]. The bands in the phonon dispersion model correspond to...

To compute the free energy using the QHA, the sum over all vibrations of the lattice is required. To ensure the correct level of accuracy a Density of States (DOS) calculation was run with increasing shrinking factors. The grid sizes used were: 1x1x1 '''(Fig. 5)''', 2x2x2 '''(Fig. 6)''', 4x4x4 '''(Fig. 7)''', 8x8x8 '''(Fig. 8)''', 16x16x16 '''(Fig. 9)''', 32x32x32 '''(Fig. 10)''' and 64x64x64 '''(Fig. 11)'''.

[[File:DOS1ce.png|400x400px|thumb|'''Fig. 5:''' Density of states using 1x1x1 grid size|centre]]

In the initial DOS calculation for the 1x1x1 grid size the graph was computed for a single k-point, corresponding to the point labelled L in k-space with coordinates (0.5, 0.5, 0.5). This can be determined from the examination of the four frequency signals in the 1x1x1 DOS graph '''(Fig. 5)''' corresponding to the frequencies that are observed at k-point L on the phonon dispersion graph '''(Fig. 4)''', and also read off from the Log file.

<div><ul>
<li style="display: inline-block;"> [[File:DOS2Ce.png|400px|thumb|'''Fig. 6:''' 2x2x2 ]]</li>
<li style="display: inline-block;"> [[File:DOS4Ce.png|400px|thumb|'''Fig. 7:''' 4x4x4]]</li>
<li style="display: inline-block;"> [[File:DOS8.png|400px|thumb|'''Fig. 8:''' 8x8x8]] </li>
<li style="display: inline-block;"> [[File:DOS16Ce.png|400px|thumb|'''Fig. 9:''' 16x16x16]] </li>
<li style="display: inline-block;"> [[File:DOS32Ce.png|400px|thumb|'''Fig. 10:''' 32x32x32]] </li>
<li style="display: inline-block;"> [[File:DOS64.png|400px|thumb|'''Fig. 11:''' 64x64x64]] </li>
</ul></div>

It is clear from the above figures that as the shrinking factor is increasing the accuracy in the Density of States graphs in also increasing. This is due to the number of k-points sampled increasing in each grid size, therefore more phonon frequencies are recorded which increases the smoothness of the graphs. However using this computation method the difference between using a shrinking factor of 32 and a shrinking factor of 64 has very little effect on the shape of the curves, this can be seen in Figure 10 and 11. Therefore the grid size of 32x32x32 was chosen for the calculation of the free-energy of MgO using the QHA as this size being the optimum balance of accuracy and CPU time in this case.

For similar crystal such as CaO, this grid size would also be appropriate as it's lattice parameters are very similar to MgO and they have the same unit cell structure. However for a Zeolite such as Faujasite a smaller grid size would be suitable. This is due to Faujasite having a much larger lattice constant, <math>a</math>: 24.638 Å <ref>Wise, W.S. (1982) New occurrence of faujasite in southeastern California. Am. Miner. 67, 794-798.</ref> compared to MgO with lattice parameter of 4.212 Å, therefore in reciprocal space the lattice parameter, <math>a^*</math> is much smaller and a smaller sample of k-points would give a sufficient degree of accuracy for the DOS calculation. When considering a metal such as Li, the environment of the electrons must be taken into account when determining the necessary grid size for the DOS calculation. In the vibrational structure of a metal, the vibrational bands are much less disperse when compared to the bands of an insulator. This can be explained by the behaviour of the electrons shielding the metal ions as they vibrate within the crystal to minimise the repulsive forces, therefore the electrons do not move as freely and are therefore have less freedom. Thus when calculating the DOS for a metal, a smaller grid size will be suitable as fewer k-points are needed.

An optimum grid size of 32x32x32 was further confirmed when computing the Free Energy of the primitive cell of MgO at 300 K shown in Table 1. As the value for the Helmholtz Free Energy can be said to be converging as grid size increases, the value for the shrinking factor of 64 can be said to be the most accurate. Therefore the Free Energy for the 2x2x2 grid size is accurate to 1 meV and 0.5 meV and the 4x4x4 grid is accurate to 0.1 meV which can be seen from the values in the final column in Table 1. However as a further accuracy can be obtained a larger grid size was preferred and 32x32x32 was chosen as it gave the same value for the Helmholtz Free Energy as 64x64x64 to 6 decimal places without taking as much CPU time.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|+ style="text-align: left;" |'''Table 1'''
!Grid Size
!Helmholtz Free Energy (eV)
!Difference in Energy from 64x64x64 grid (eV)
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|3.82x10<sup>-3</sup>
|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|1.26x10<sup>-4</sup>
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|3.30x10<sup>-5</sup>
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|5.00x10<sup>-6</sup>
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|1.00x10<sup>-6</sup>
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|
|}

Using the Quasi-harmonic approximation the Helmholtz Free Energy, lattice constant,<math>a</math>, and cell volume of a primitive lattice of MgO were calculated and recorder at 100 K intervals within the range 0 K to 1000 K '''(Table 2)'''.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|+ style="text-align: left;" |'''Table 2'''
!Temperature (K) !! Helmholtz Free Energy (eV) !! Lattice constant (Å) !! Cell volume (Å<sup>3</sup>)
|-
|0 || -40.901906 || 2.986563 || 18.838260
|-
|100 || -40.902420 || 2.986563 || 18.838260
|-
|200 || -40.909377 || 2.987604 || 18.856201
|-
|300 || -40.928125 || 2.989390 || 18.890025
|-
|400 || -40.958595 || 2.991629 || 18.932507
|-
|500 || -40.999436 || 2.994134 || 18.980110
|-
|600 || -41.049316 || 2.996820 || 19.031221
|-
|700 || -41.107120 || 2.999643 || 19.085057
|-
|800 || -41.171892 || 3.002587 || 19.141316
|-
|900 || -41.243018 || 3.005634 || 19.199638
|-
|1000 || -41.319849 || 3.008783 || 19.260042
|}

The value for the Helmholtz Free Energy is becoming more negative as the temperature is increased which is inline with the equation used by the QHA to calculate the Free Energy (F):

<math>F\left(T, V\right) = U\left(V\right) - TS\left(T, V\right)</math>

As temperature, <math>T</math>, increases the term <math>TS\left(T, V\right)</math> becomes more negative and dominates the equation causing the free energy, <math>F</math> to decrease'''(Fig. 12)'''. Conversely the lattice constant, <math>a</math>, increases with temperature '''(Fig. 13)''' as does the volume '''(Fig.14)'''. The volume of the cell increases with temperature as the higher the temperature the more kinetic energy which can manifest as vibrational energy within the crystal. As this energy increases the amplitude of the crystal vibrations is increasing causing the equilibrium bond length to become greater as the repulsive forces become stronger at shorter distances therefore the distance between the atoms becomes greater<ref>Ziman, J.M. (1972). Principles of the Theory of Solids, 2nd Ed. Cambridge University Press. ISBN 0521297338 9780521297332</ref>. This is the origin of thermal expansion within this approximation.

<div><ul>
<li style="display: inline-block;"> [[File:FreeEnergyQHA.png|500px|thumb|'''Fig. 12''']]</li>
<li style="display: inline-block;"> [[File:Latticeconst.png|500px|thumb|'''Fig. 13''']]</li>
</ul></div>

[[File:Thermalexpce.png|800px|thumb|'''Fig. 14'''|centre]]
A line of best fit was plotted for the linear region of the data (400 - 1000 K) and its gradient taken for the rate of change of volume with temperature in '''Eq. 1''' and <math>V_0</math> taken as the volume of the cell at 400 K.

<math>\alpha = \frac{1}{18.932507}0.0005</math>

<math>\alpha = 2.640960</math>x<math>10^{-5}</math> <math>K^{-1}</math> (lit.)<ref>email ref</ref>

The thermal expansion coefficient was found using molecular dynamics. This computational approach models the movements of the atoms in the lattice and calculates the thermodynamic properties from the time-averaged behaviour. A supercell made up of 32 conventional lattice cell is used due to the movement of the atoms no longer being periodical therefore more unit cells must be used to model this as accurately as possible. The larger the system in the computation the closer to a real system it becomes and the more accurate the values produced. However due to this computation taking much higher CPU time than when calculating in the QHA only a system of 32 conventional cells were used.

The volume change was plotted against temperature '''(Fig 15)'''.

[[File:MDTherm.png|800px|thumb|'''Fig. 15'''|centre]]

The volume of the system increases linearly, a value for 0 K cannot be calculated for this approximation as at 0 K there would be 0 kinetic energy therefore the atom would be at rest; due to the uncertainty principle <ref>W. Heisenberg, Remarks on the origin of the relations of uncertainty, in: W. Price, S. Chissick (Eds.), The Uncertainty Principle and Foundations of Quantum Mechanics. A Fifty Years’ Survey, J. Wiley & Sons, London, 1977, pp. 3–6.</ref> both the position and momentum of a particle cannot be known therefore the computation breaks down here.

As with the QHA a line of best fit was plotted for the data and its gradient used for the rate of change of volume with temperature in '''Eq. 1''' and <math>V_0</math> the volume at 100 K.

<math>\alpha = \frac{1}{599.513295}0.0173</math>

<math>\alpha = 2.668831</math>x<math>10^{-5}</math> <math>K^{-1}</math>

==Conclusion==

Both the Quasi-harmonic approximation and the Molecular Dynamics computations gave values for the thermal expansion coefficient: 2.640960x10<sup>-5</sup> K<sup>-1</sup> and 2.668831x10<sup>-5</sup>
The results for the change in volume with respect to temperature is compared for the QHA and for MD '''Fig. 16'''.

[[File:Comparevol.png|800px|thumb|'''Fig. 16'''|centre]]

How does the thermal expansion predicted by MD compare to that predicted by the quasi-harmonic approximation ?
Why do the two methods produce different answers ? - how does the difference depend on temperature ?

What would happen to the cell volume at high temperature in the quasi-harmonic approximation and in the MD ?

As the temperature approaches the melting point of MgO how well do the phonon modes represent the actual motions of the ions ?
In a diatomic molecule with an exactly harmonic potential would you expect the bond length to increase with temperature ? Why does it increase in the solid when we are using an quasi-harmonic approximation ?

==References==

Rep:Mod:CE1014MgO

2017-02-09T19:55:10Z

Ce1014:

=Thermal Expansion of MgO=

==Introduction==

This investigation involves modelling the thermal expansion of a crystal of MgO, and calculating its thermal expansion coefficient <math>\alpha </math>:

<math>\alpha = \frac{1}{V_0}\left(\frac{\partial V}{\partial T}\right)_P </math> '''Eq. 1'''

The Quasi-Harmonic Approximation is used to calculate the vibrational energy levels and to understand the phonon dispersion of the material and its vibrational density of states. This approach is compared to the Molecular dynamics results used to simulate the vibrations as random motions of atoms within a cell.

To calculate the thermal expansion coefficient the lattice vibrations must be calculated; the main contributor to thermal expansion. To calculate how these vibrations change with temperature and therefore how the separation between the atoms change, giving rise to thermal expansion, the harmonic potential model is insufficient. This is due to the equilibrium distance between atoms being invariant within this model as the displacement of atoms is symmetrical about their equilibrium positions. Instead the quasi-harmonic approximation (QHA), which is a phonon-based model, is used to model the thermal expansion by assuming the frequencies are volume dependent: <math>\omega_j\left(V\right)</math>

In this model the free energy is:

<math>F=</math> (free energy from vibrational theory)

The k within the equations above is the wave vector, essentially a label that can be given to each unique point within the unit cell which corresponds to a unique vibrational energy within the lattice. To find the k-values the crystal must be converted into reciprocal space with lattice parameters <math>a^*</math> found using the real space lattice parameter <math>a</math>:

<math> a^*=\frac{2\pi}{a}</math>

Within the reciprocal cell <math>k</math> is periodic between <math>-\frac{\pi}{a}</math> and <math>\frac{\pi}{a}</math>.

The frequencies of the vibrations <math>\omega</math> vary with <math>k</math> by the following relationship:

<math>k=</math>

The crystal structure of MgO in a conventional cell is a face-centered cubic lattice which consists of a total of 8 atoms '''(Fig 1)''', it can also be modelled as a primitive cell '''(Fig. 2)'''; the simplest unit cell, with one ion in the centre and the opposing ion on each corner giving a total of 2 atoms within the lattice. The conventional cell has 4 times the volume of the primitive unit cell however is a more realistic representation of the symmetry of the crystal therefore is still used in many crystallographic studies. In this investigation the primitive cell will be dealt with initially when calculating using the QHA, but a super cell '''(Fig. 3)''' is used for the molecular dynamics simulations which is made up of 32 conventional cells.

{|style="margin: 0 auto;"
| [[File:Convencell.png|thumb|'''Fig. 1:'''Conventional cell of MgO|upright]]
| [[File:Primcell.png|thumb|'''Fig. 2:''' Primitive cell of MgO|upright]]
| [[File:Supercell.png|thumb|'''Fig. 3:''' Supercell of MgO]]
|}

==Method==

All calculations were carried out using DLVisualize, a graphical interface used for modelling. Within this programme the modelling code GULP (General Utility Lattice Program) was utilised to calculate the theoretical values evaluated in the Results and Discussion below.

The internal energy (<math>U</math>) of the primitive unit cell was initially calculated by modelling the crystal as a purely ionic material which was used in future calculations of the Helmholtz Free Energy (<math>A</math>).

To enable us to visualise how the frequencies of vibrations of this lattice vary with <math>k</math> a phonon dispersion curve was calculated using a sample of 50 k-points and represented graphically '''(FIg. 4)'''. Subsequently a calculation was run to find the density of states for a lattice grid size of 1x1x1. This calculation was then repeated, doubling the shrinking factor each time until a grid size of 64x64x64 was obtained.

'''Initially the thermal expansion coefficient was calculated using the QHA:'''

Using the optimum grid size, found to be 32x32x32 from the above simulations, the Helmholtz Free Energy was found for the MgO crystal at temperatures between 0 K and 1000 K in 100 K intervals. The lattice parameters and the total volume for each cell was recorded at each temperature and plotted [Fig ....]

From the plot of cell volume with respect to temperature, the gradient of the linear section was taken with the initial volume and used in '''Eq. 1''' to give a value of the thermal expansion coefficient, <math>\alpha</math>, as <math>2.64</math>x<math>10^{-5} K^{-1}</math> for a temperature range 400 - 1000 K.

'''The thermal expansion coefficient was then calculated using Molecular Dynamics.'''

Within this computation the system has a timestep of 1 femtosecond <math>(10^{-15} s)</math> and is run for 500 steps after equilibrium has been reached at the selected temperature. The same temperature range from 0 K to 1000 K was used as above. The system used was made up of 32 MgO conventional cells. The average equilibrium volumes were recorded at each temperature for this simulation [Fig ...]. '''Eq. 1''' was used to calculate the thermal expansion coefficient as above using the gradient of the volume as a function of temperature line of best fit and the initial volume found. Giving a value of <math>\alpha</math> as <math>2.86</math>x<math>10^{-5} K^{-1}</math> with a temperature range 100 - 1000 K.

==Results and Discussion==

The internal energy per primitive cell of MgO was found to be -41.0753 eV. MgO was modelled as purely ionic crystal with Mg 2+ ions and O 2- ions having simple Buckingham repulsive potentials operating between the ions. <ref>T.S.Bush, J.D.Gale, C.R.A.Catlow and P.D. Battle J. Mater Chem., 4, 831-837 (1994) </ref>, this energy is also described as the 'binding energy' of the crystal as it is the energy necessary to separate the ions to an infinite distance.

The Phonon dispersion curve is shown in Fig. 4 with k-space plotted against frequency of the vibrations.

[[File:Phonon_Dispersion.png|500x500px|thumb|Fig. 4: Phonon dispersion curve for MgO|centre]]

The labels on the x-axis are points of interest in k-space. For example the origin is labelled <math>\Gamma</math> with lattice coordinates [0, 0, 0]. The bands in the phonon dispersion model correspond to...

To compute the free energy using the QHA, the sum over all vibrations of the lattice is required. To ensure the correct level of accuracy a Density of States (DOS) calculation was run with increasing shrinking factors. The grid sizes used were: 1x1x1 '''(Fig. 5)''', 2x2x2 '''(Fig. 6)''', 4x4x4 '''(Fig. 7)''', 8x8x8 '''(Fig. 8)''', 16x16x16 '''(Fig. 9)''', 32x32x32 '''(Fig. 10)''' and 64x64x64 '''(Fig. 11)'''.

[[File:DOS1ce.png|400x400px|thumb|'''Fig. 5:''' Density of states using 1x1x1 grid size|centre]]

In the initial DOS calculation for the 1x1x1 grid size the graph was computed for a single k-point, corresponding to the point labelled L in k-space with coordinates (0.5, 0.5, 0.5). This can be determined from the examination of the four frequency signals in the 1x1x1 DOS graph '''(Fig. 5)''' corresponding to the frequencies that are observed at k-point L on the phonon dispersion graph '''(Fig. 4)''', and also read off from the Log file.

<div><ul>
<li style="display: inline-block;"> [[File:DOS2Ce.png|400px|thumb|'''Fig. 6:''' 2x2x2 ]]</li>
<li style="display: inline-block;"> [[File:DOS4Ce.png|400px|thumb|'''Fig. 7:''' 4x4x4]]</li>
<li style="display: inline-block;"> [[File:DOS8.png|400px|thumb|'''Fig. 8:''' 8x8x8]] </li>
<li style="display: inline-block;"> [[File:DOS16Ce.png|400px|thumb|'''Fig. 9:''' 16x16x16]] </li>
<li style="display: inline-block;"> [[File:DOS32Ce.png|400px|thumb|'''Fig. 10:''' 32x32x32]] </li>
<li style="display: inline-block;"> [[File:DOS64.png|400px|thumb|'''Fig. 11:''' 64x64x64]] </li>
</ul></div>

It is clear from the above figures that as the shrinking factor is increasing the accuracy in the Density of States graphs in also increasing. This is due to the number of k-points sampled increasing in each grid size, therefore more phonon frequencies are recorded which increases the smoothness of the graphs. However using this computation method the difference between using a shrinking factor of 32 and a shrinking factor of 64 has very little effect on the shape of the curves, this can be seen in Figure 10 and 11. Therefore the grid size of 32x32x32 was chosen for the calculation of the free-energy of MgO using the QHA as this size being the optimum balance of accuracy and CPU time in this case.

For similar crystal such as CaO, this grid size would also be appropriate as it's lattice parameters are very similar to MgO and they have the same unit cell structure. However for a Zeolite such as Faujasite a smaller grid size would be suitable. This is due to Faujasite having a much larger lattice constant, <math>a</math>: 24.638 Å <ref>Wise, W.S. (1982) New occurrence of faujasite in southeastern California. Am. Miner. 67, 794-798.</ref> compared to MgO with lattice parameter of 4.212 Å, therefore in reciprocal space the lattice parameter, <math>a^*</math> is much smaller and a smaller sample of k-points would give a sufficient degree of accuracy for the DOS calculation. When considering a metal such as Li, the environment of the electrons must be taken into account when determining the necessary grid size for the DOS calculation. In the vibrational structure of a metal, the vibrational bands are much less disperse when compared to the bands of an insulator. This can be explained by the behaviour of the electrons shielding the metal ions as they vibrate within the crystal to minimise the repulsive forces, therefore the electrons do not move as freely and are therefore have less freedom. Thus when calculating the DOS for a metal, a smaller grid size will be suitable as fewer k-points are needed.

An optimum grid size of 32x32x32 was further confirmed when computing the Free Energy of the primitive cell of MgO at 300 K shown in Table 1. As the value for the Helmholtz Free Energy can be said to be converging as grid size increases, the value for the shrinking factor of 64 can be said to be the most accurate. Therefore the Free Energy for the 2x2x2 grid size is accurate to 1 meV and 0.5 meV and the 4x4x4 grid is accurate to 0.1 meV which can be seen from the values in the final column in Table 1. However as a further accuracy can be obtained a larger grid size was preferred and 32x32x32 was chosen as it gave the same value for the Helmholtz Free Energy as 64x64x64 to 6 decimal places without taking as much CPU time.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|+ style="text-align: left;" |'''Table 1'''
!Grid Size
!Helmholtz Free Energy (eV)
!Difference in Energy from 64x64x64 grid (eV)
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|3.82x10<sup>-3</sup>
|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|1.26x10<sup>-4</sup>
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|3.30x10<sup>-5</sup>
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|5.00x10<sup>-6</sup>
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|1.00x10<sup>-6</sup>
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|
|}

Using the Quasi-harmonic approximation the Helmholtz Free Energy, lattice constant,<math>a</math>, and cell volume of a primitive lattice of MgO were calculated and recorder at 100 K intervals within the range 0 K to 1000 K '''(Table 2)'''.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|+ style="text-align: left;" |'''Table 2'''
!Temperature (K) !! Helmholtz Free Energy (eV) !! Lattice constant (Å) !! Cell volume (Å<sup>3</sup>)
|-
|0 || -40.901906 || 2.986563 || 18.838260
|-
|100 || -40.902420 || 2.986563 || 18.838260
|-
|200 || -40.909377 || 2.987604 || 18.856201
|-
|300 || -40.928125 || 2.989390 || 18.890025
|-
|400 || -40.958595 || 2.991629 || 18.932507
|-
|500 || -40.999436 || 2.994134 || 18.980110
|-
|600 || -41.049316 || 2.996820 || 19.031221
|-
|700 || -41.107120 || 2.999643 || 19.085057
|-
|800 || -41.171892 || 3.002587 || 19.141316
|-
|900 || -41.243018 || 3.005634 || 19.199638
|-
|1000 || -41.319849 || 3.008783 || 19.260042
|}

The value for the Helmholtz Free Energy is becoming more negative as the temperature is increased which is inline with the equation used by the QHA to calculate the Free Energy (F):

<math>F\left(T, V\right) = U\left(V\right) - TS\left(T, V\right)</math>

As temperature, <math>T</math>, increases the term <math>TS\left(T, V\right)</math> becomes more negative and dominates the equation causing the free energy, <math>F</math> to decrease'''(Fig. 12)'''. Conversely the lattice constant, <math>a</math>, increases with temperature '''(Fig. 13)''' as does the volume '''(Fig.14)'''. The volume of the cell increases with temperature as the higher the temperature the more kinetic energy which can manifest as vibrational energy within the crystal. As this energy increases the amplitude of the crystal vibrations is increasing causing the equilibrium bond length to become greater as the repulsive forces become stronger at shorter distances therefore the distance between the atoms becomes greater<ref>Ziman, J.M. (1972). Principles of the Theory of Solids, 2nd Ed. Cambridge University Press. ISBN 0521297338 9780521297332</ref>. This is the origin of thermal expansion within this approximation.

<div><ul>
<li style="display: inline-block;"> [[File:FreeEnergyQHA.png|500px|thumb|'''Fig. 12''']]</li>
<li style="display: inline-block;"> [[File:Latticeconst.png|500px|thumb|'''Fig. 13''']]</li>
</ul></div>

[[File:Thermalexpce.png|800px|thumb|'''Fig. 14'''|centre]]
A line of best fit was plotted for the linear region of the data (400 - 1000 K) and its gradient taken for the rate of change of volume with temperature in '''Eq. 1''' and <math>V_0</math> taken as the volume of the cell at 400 K.

<math>\alpha = \frac{1}{18.932507}0.0005</math>

<math>\alpha = 2.640960</math>x<math>10^{-5}</math> <math>K^{-1}</math> (lit.)<ref>email ref</ref>

The thermal expansion coefficient was found using molecular dynamics. This computational approach models the movements of the atoms in the lattice and calculates the thermodynamic properties from the time-averaged behaviour. A supercell made up of 32 conventional lattice cell is used due to the movement of the atoms no longer being periodical therefore more unit cells must be used to model this as accurately as possible. The larger the system in the computation the closer to a real system it becomes and the more accurate the values produced. However due to this computation taking much higher CPU time than when calculating in the QHA only a system of 32 conventional cells were used.

The volume change was plotted against temperature '''(Fig 15)'''.

[[File:MDTherm.png|800px|thumb|'''Fig. 15'''|centre]]

The volume of the system increases linearly, a value for 0 K cannot be calculated for this approximation as at 0 K there would be 0 kinetic energy therefore the atom would be at rest; due to the uncertainty principle <ref>W. Heisenberg, Remarks on the origin of the relations of uncertainty, in: W. Price, S. Chissick (Eds.), The Uncertainty Principle and Foundations of Quantum Mechanics. A Fifty Years’ Survey, J. Wiley & Sons, London, 1977, pp. 3–6.</ref> both the position and momentum of a particle cannot be known therefore the computation breaks down here.

As with the QHA a line of best fit was plotted for the data and its gradient used for the rate of change of volume with temperature in '''Eq. 1''' and <math>V_0</math> the volume at 100 K.

<math>\alpha = \frac{1}{599.513295}0.0173</math>

<math>\alpha = 2.668831</math>x<math>10^{-5}</math> <math>K^{-1}</math>

==Conclusion==

Both the Quasi-harmonic approximation
The results for the change in volume with respect to temperature is compared for the QHA and for MD '''Fig. 16'''.

[[File:Comparevol.png|800px|thumb|'''Fig. 16'''|centre]]

How does the thermal expansion predicted by MD compare to that predicted by the quasi-harmonic approximation ?
Why do the two methods produce different answers ? - how does the difference depend on temperature ?

What would happen to the cell volume at high temperature in the quasi-harmonic approximation and in the MD ?

As the temperature approaches the melting point of MgO how well do the phonon modes represent the actual motions of the ions ?
In a diatomic molecule with an exactly harmonic potential would you expect the bond length to increase with temperature ? Why does it increase in the solid when we are using an quasi-harmonic approximation ?

==References==

File:Comparevol.png

2017-02-09T19:49:14Z

Ce1014:

2017-02-08T13:06:54Z

Ce1014:

=Thermal Expansion of MgO=

==Introduction==

[[File:Primcell.png|265x265px|thumb|Primitive Cell of MgO]][[File:Convencell.png|265x265px|thumb|Conventional Cell of MgO]][[File:Supercell.png|265x265px|thumb|Supercell of MgO]]

This investigation involves modelling the thermal expansion of a crystal of MgO, and calculating its thermal expansion coefficient <math>\alpha </math>:

<math>\alpha = \frac{1}{V_0}\left(\frac{\partial V}{\partial T}\right)_P </math> '''Eq. 1'''

The Quasi-Harmonic Approximation is used to calculate the vibrational energy levels and to understand the phonon dispersion of the material and its vibrational density of states. This approach is compared to the Molecular dynamics results used to simulate the vibrations as random motions of atoms within a cell.

To calculate the thermal expansion coefficient the lattice vibrations must be calculated; the main contributor to thermal expansion. To calculate how these vibrations change with temperature and therefore how the separation between the atoms change, giving rise to thermal expansion, the harmonic potential model is insufficient. This is due to the equilibrium distance between atoms being invariant within this model as the displacement of atoms is symmetrical about their equilibrium positions. Instead the quasi-harmonic approximation (QHA), which is a phonon-based model, is used to model the thermal expansion by assuming the frequencies are volume dependent: <math>\omega_j\left(V\right)</math>

In this model the free energy is:

<math>F=</math> (free energy from vibrational theory)

The '''k''' within the equations above is the wave vector, essentially a label that can be given to each unique point within the unit cell which corresponds to a unique vibrational energy within the lattice. To find the k-values the crystal must be converted into reciprocal space with lattice parameters '''a*''' found using the real space lattice parameter '''a''':

<math> a^*=\frac{2\pi}{a}</math>

Within the reciprocal cell <math>k</math> is periodic between <math>-\frac{\pi}{a}</math> and <math>\frac{\pi}{a}</math>.

The frequencies of the vibrations <math>\omega</math> vary with <math>k</math> by the following relationship:

<math>k=</math>

The crystal structure of MgO in a conventional cell is a face-centered cubic lattice which consists of a total of 8 atoms [Fig 1], it can also be modelled as a primitive cell; the simplest unit cell, with one ion in the centre and the opposing ion on each corner giving a total of 2 atoms within the lattice, which can be repeated to give the overall crystal structure. The conventional cell has 4 times the volume of the primitive unit cell however is a more realistic representation of the symmetry of the crystal therefore is still used in many crystallographic studies. In this investigation the primitive cell will be dealt with initially when calculating using the QHA but a super cell is used for the molecular dynamics simulations which is made up of 32 conventional cells.

==Method==

All calculations were carried out using DLVisualize, a graphical interface used for modelling. Within this programme the modelling code GULP (General Utility Lattice Program) was utilised to calculate the theoretical values evaluated in the Results and Discussion below.

The internal energy (<math>U</math>) of the primitive unit cell was initially calculated by modelling the crystal as a purely ionic material which was used in future calculations of the Helmholtz Free Energy (<math>A</math>).

To enable us to visualise how the frequencies of vibrations of this lattice vary with <math>k</math> a phonon dispersion curve was calculated using a sample of 50 k-points and represented graphically '''(FIg. 4)'''. Subsequently a calculation was run to find the density of states for a lattice grid size of 1x1x1. This calculation was then repeated, doubling the shrinking factor each time until a grid size of 64x64x64 was obtained.

'''Initially the thermal expansion coefficient was calculated using the QHA:'''

Using the optimum grid size, found to be 32x32x32 from the above simulations, the Helmholtz Free Energy was found for the MgO crystal at temperatures between 0 K and 1000 K in 100 K intervals. The lattice parameters and the total volume for each cell was recorded at each temperature and plotted [Fig ....]

From the plot of cell volume with respect to temperature, the gradient of the linear section was taken with the initial volume and used in '''Eq. 1''' to give a value of the thermal expansion coefficient, <math>\alpha</math>, as <math>2.64</math>x<math>10^{-5} K^{-1}</math> for a temperature range 400 - 1000 K.

The thermal expansion coefficient was then calculated using Molecular Dynamics.

Within this computation the system has a timestep of 1 femtosecond <math>(10^{-15} s)</math> and is run for 500 steps after equilibrium has been reached at the selected temperature. The same temperature range from 0 K to 1000 K was used as above. The system used was made up of 32 MgO conventional cells. The average equilibrium volumes were recorded at each temperature for this simulation [Fig ...]. '''Eq. 1''' was used to calculate the thermal expansion coefficient as above using the gradient of the volume as a function of temperature line of best fit and the initial volume found. Giving a value of <math>\alpha</math> as <math>2.86</math>x<math>10^{-5} K^{-1}</math> with a temperature range 100 - 1000 K.

==Results and Discussion==

The internal energy per primitive cell of MgO was found to be -41.0753 eV. MgO was modelled as purely ionic crystal with Mg 2+ ions and O 2- ions having simple Buckingham repulsive potentials operating between the ions. <ref>T.S.Bush, J.D.Gale, C.R.A.Catlow and P.D. Battle J. Mater Chem., 4, 831-837 (1994) </ref>, this energy is also described as the 'binding energy' of the crystal as it is the energy necessary to separate the ions to an infinite distance.

The Phonon dispersion curve is shown in Fig. 4 with k-space plotted against frequency of the vibrations.

[[File:Phonon_Dispersion.png|500x500px|thumb|Fig. 4: Phonon dispersion curve for MgO|centre]]

The labels on the x-axis are points of interest in k-space. For example the origin is labelled <math>\Gamma</math> with lattice coordinates [0, 0, 0].

To compute the free energy using the QHA, the sum over all vibrations of the lattice is required. To ensure the correct level of accuracy a Density of States (DOS) calculation was run with increasing shrinking factors. The grid sizes used were: 1x1x1 '''(Fig. 5)''', 2x2x2 '''(Fig. 6)''', 4x4x4 '''(Fig. 7)''', 8x8x8 '''(Fig. 8)''', 16x16x16 '''(Fig. 9)''', 32x32x32 '''(Fig. 10)''' and 64x64x64 '''(Fig. 11)'''.

[[File:DOS1ce.png|400x400px|thumb|'''Fig. 5:''' DOS for MgO with 1x1x1 grid]]

For the original 1x1x1 DOS graph, the k-point corresponds to the point labelled L on the phonon dispersion curve at coordinates 0.5, 0.5, 0.5 which can be seen as the wavelengths detected for the DOS calculation correspond the wavelengths seen at point L on the phonon dispersion graph (wavelengths).
As the shrinking factors increase the DOS graphs are smoothing off as more phonon vibrations are taken into account. However the difference between the 32x32x32 grid and the 64x64x64 grid is minimal, with the 64x grid taking a greater time to run for little more accuaracy.
The corresponding free energy at 300K for each of these grid sizes can be seen in Table 1. The difference between 32x and 64x is ….. which is a close enough accuracy to continue the free energy calculations with 32x.

The internal energy of a crystal of MgO

'''Calculating phonon modes of MgO'''

• The density of states for 1x1x1 grid was computed for a single k-point.

• Can you work out which k-point by examining the dispersion curves ? look in the Log File to check.

The k-point used to calculate the Density of states was at the cell coordinate 0.5, 0.5, 0.5 which is also the k-point label 'L' seen on the phonon dispersion graph below.

[[File:DOS2Ce.png|500px]]

[[File:DOS4Ce.png|500px]]

[[File:DOS8.png|500px]]

[[File:DOS16Ce.png|500px]]

[[File:DOS32Ce.png|500px]]

[[File:DOS64.png|500px]]

• How does the density of states vary with grid size ?

As the grid size increases, therefore the k points analysed increases the density of states graphs vary as seen below:

• As the grid size increases more and more of the possible vibrations are sampled and more features appear - which grid size is the minimum for a reasonable approximation to the density of states ?

grid size 16x16x16

• How does the density of states computed at this optimal grid size compare to that computed for smaller and larger grids ?

Smaller grids don't have enough resolution, above 16 graph smooths off but has no further resolution.

• How is the density of states related to the dispersion curves ?

==References==

Rep:Mod:CE1014MgO

2017-02-08T12:49:11Z

Ce1014: /* Results and Discussion */

=Thermal Expansion of MgO=

==Introduction==

[[File:Primcell.png|265x265px|thumb|Primitive Cell of MgO]][[File:Convencell.png|265x265px|thumb|Conventional Cell of MgO]][[File:Supercell.png|265x265px|thumb|Supercell of MgO]]

This investigation involves modelling the thermal expansion of a crystal of MgO, and calculating its thermal expansion coefficient <math>\alpha </math> <ref>Paper on Blackboard</ref>:

<math>\alpha = \frac{1}{V_0}\left(\frac{\partial V}{\partial T}\right)_P </math> '''Eq. 1'''

The Quasi-Harmonic Approximation is used to calculate the vibrational energy levels and to understand the phonon dispersion of the material and its vibrational density of states. This approach is compared to the Molecular dynamics results used to simulate the vibrations as random motions of atoms within a cell.

To calculate the thermal expansion coefficient the lattice vibrations must be calculated; the main contributor to thermal expansion. To calculate how these vibrations change with temperature and therefore how the separation between the atoms change, giving rise to thermal expansion, the harmonic potential model is insufficient. This is due to the equilibrium distance between atoms being invariant within this model as the displacement of atoms is symmetrical about their equilibrium positions. Instead the quasi-harmonic approximation (QHA), which is a phonon-based model, is used to model the thermal expansion by assuming the frequencies are volume dependent: <math>\omega_j\left(V\right)</math>

In this model the free energy is:

<math>F=</math> (free energy from vibrational theory)

The '''k''' within the equations above is the wave vector, essentially a label that can be given to each unique point within the unit cell which corresponds to a unique vibrational energy within the lattice. To find the k-values the crystal must be converted into reciprocal space with lattice parameters '''a*''' found using the real space lattice parameter '''a''':

<math> a^*=\frac{2\pi}{a}</math>

Within the reciprocal cell <math>k</math> is periodic between <math>-\frac{\pi}{a}</math> and <math>\frac{\pi}{a}</math>.

The frequencies of the vibrations <math>\omega</math> vary with <math>k</math> by the following relationship:

<math>k=</math>

The crystal structure of MgO in a conventional cell is a face-centered cubic lattice which consists of a total of 8 atoms [Fig 1], it can also be modelled as a primitive cell; the simplest unit cell, with one ion in the centre and the opposing ion on each corner giving a total of 2 atoms within the lattice, which can be repeated to give the overall crystal structure. The conventional cell has 4 times the volume of the primitive unit cell however is a more realistic representation of the symmetry of the crystal therefore is still used in many crystallographic studies. In this investigation the primitive cell will be dealt with initially when calculating using the QHA but a super cell is used for the molecular dynamics simulations which is made up of 32 conventional cells.

==Method==

All calculations were carried out using DLVisualize, a graphical interface used for modelling. Within this programme the modelling code GULP (General Utility Lattice Program) was utilised to calculate the theoretical values evaluated in the Results and Discussion below.

Initially the thermal expansion coefficient was calculated using the QHA.
The internal energy (<math>U</math>) of the primitive unit cell was initially calculated by modelling the crystal as a purely ionic material which was used in future calculations of the Helmholtz Free Energy (<math>A</math>).

To enable us to visualise how the frequencies of vibrations of this lattice vary with <math>k</math> a phonon dispersion curve was calculated using a sample of 50 k points and represented graphically [Fig ...]. Subsequently a calculation was run to find the density of states for a lattice grid size of 1x1x1 then repeated, doubling the shrinking factor each time until a grid size of 64x64x64 was obtained.

Using the optimum grid size, found to be 32x32x32 from the above simulations, the Helmholtz Free Energy was found for the MgO crystal at temperatures between 0 K and 1000 K in 100 K intervals. The lattice parameters and the total volume for each cell was recorded at each temperature and plotted [Fig ....]

From the plot of cell volume with respect to temperature, the gradient of the linear section was taken with the initial volume and used in '''Eq. 1''' to give a value of the thermal expansion coefficient, <math>\alpha</math>, as <math>2.64</math>x<math>10^{-5} K^{-1}</math> for a temperature range 400 - 1000 K.

The thermal expansion coefficient was then calculated using Molecular Dynamics.

Within this computation the system has a timestep of 1 femtosecond <math>(10^{-15} s)</math> and is run for 500 steps after equilibrium has been reached at the selected temperature. The same temperature range from 0 K to 1000 K was used as above. The system used was made up of 32 MgO conventional cells. The average equilibrium volumes were recorded at each temperature for this simulation [Fig ...]. '''Eq. 1''' was used to calculate the thermal expansion coefficient as above using the gradient of the volume as a function of temperature line of best fit and the initial volume found. Giving a value of <math>\alpha</math> as <math>2.86</math>x<math>10^{-5} K^{-1}</math> with a temperature range 100 - 1000 K.

==Results and Discussion==

The internal energy per primitive cell of MgO was found to be -41.0753 eV. MgO was modelled as purely ionic crystal with Mg 2+ ions and O 2- ions having simple Buckingham repulsive potentials operating between the ions. <ref>T.S.Bush, J.D.Gale, C.R.A.Catlow and P.D. Battle J. Mater Chem., 4, 831-837 (1994) </ref>.

On the X axis are points of interest along the conventional k-space for example () corresponds to centre of the lattice with coordinates (0, 0, 0) and L corresponding to the corner of the primitive cell with coordinates (0.5, 0.5, 0.5).
To compute the free energy, the sum over all vibrations of the lattice is required. To ensure the correct level of accuracy a Density of States (DOS) calculation was run with increasing shrinking factors from 1x1x1 doubling up until 64x64x64.
For the original 1x1x1 DOS graph, the k-point corresponds to the point labelled L on the phonon dispersion curve at coordinates 0.5, 0.5, 0.5 which can be seen as the wavelengths detected for the DOS calculation correspond the wavelengths seen at point L on the phonon dispersion graph (wavelengths).
As the shrinking factors increase the DOS graphs are smoothing off as more phonon vibrations are taken into account. However the difference between the 32x32x32 grid and the 64x64x64 grid is minimal, with the 64x grid taking a greater time to run for little more accuaracy.
The corresponding free energy at 300K for each of these grid sizes can be seen in Table 1. The difference between 32x and 64x is ….. which is a close enough accuracy to continue the free energy calculations with 32x.

The internal energy of a crystal of MgO

'''Calculating phonon modes of MgO'''

• The density of states for 1x1x1 grid was computed for a single k-point.

[[File:DOS1ce.png|500px]]

• Can you work out which k-point by examining the dispersion curves ? look in the Log File to check.

The k-point used to calculate the Density of states was at the cell coordinate 0.5, 0.5, 0.5 which is also the k-point label 'L' seen on the phonon dispersion graph below.

[[File:Phonon_Dispersion.png|500px]]

[[File:DOS2Ce.png|500px]]

[[File:DOS4Ce.png|500px]]

[[File:DOS8.png|500px]]

[[File:DOS16Ce.png|500px]]

[[File:DOS32Ce.png|500px]]

[[File:DOS64.png|500px]]

• How does the density of states vary with grid size ?

As the grid size increases, therefore the k points analysed increases the density of states graphs vary as seen below:

• As the grid size increases more and more of the possible vibrations are sampled and more features appear - which grid size is the minimum for a reasonable approximation to the density of states ?

grid size 16x16x16

• How does the density of states computed at this optimal grid size compare to that computed for smaller and larger grids ?

Smaller grids don't have enough resolution, above 16 graph smooths off but has no further resolution.

• How is the density of states related to the dispersion curves ?

==References==

Rep:Mod:CE1014MgO

2017-02-08T11:53:16Z

Ce1014:

=Thermal Expansion of MgO=

==Introduction==

[[File:Primcell.png|265x265px|thumb|Primitive Cell of MgO]][[File:Convencell.png|265x265px|thumb|Conventional Cell of MgO]][[File:Supercell.png|265x265px|thumb|Supercell of MgO]]

This investigation involves modelling the thermal expansion of a crystal of MgO, and calculating its thermal expansion coefficient <math>\alpha </math> <ref>Paper on Blackboard</ref>:

<math>\alpha = \frac{1}{V_0}\left(\frac{\partial V}{\partial T}\right)_P </math> '''Eq. 1'''

The Quasi-Harmonic Approximation is used to calculate the vibrational energy levels and to understand the phonon dispersion of the material and its vibrational density of states. This approach is compared to the Molecular dynamics results used to simulate the vibrations as random motions of atoms within a cell.

To calculate the thermal expansion coefficient the lattice vibrations must be calculated; the main contributor to thermal expansion. To calculate how these vibrations change with temperature and therefore how the separation between the atoms change, giving rise to thermal expansion, the harmonic potential model is insufficient. This is due to the equilibrium distance between atoms being invariant within this model as the displacement of atoms is symmetrical about their equilibrium positions. Instead the quasi-harmonic approximation (QHA), which is a phonon-based model, is used to model the thermal expansion by assuming the frequencies are volume dependent: <math>\omega_j\left(V\right)</math>

In this model the free energy is:

<math>F=</math> (free energy from vibrational theory)

The '''k''' within the equations above is the wave vector, essentially a label that can be given to each unique point within the unit cell which corresponds to a unique vibrational energy within the lattice. To find the k-values the crystal must be converted into reciprocal space with lattice parameters '''a*''' found using the real space lattice parameter '''a''':

<math> a^*=\frac{2\pi}{a}</math>

Within the reciprocal cell <math>k</math> is periodic between <math>-\frac{\pi}{a}</math> and <math>\frac{\pi}{a}</math>.

The frequencies of the vibrations <math>\omega</math> vary with <math>k</math> by the following relationship:

<math>k=</math>

The crystal structure of MgO in a conventional cell is a face-centered cubic lattice which consists of a total of 8 atoms [Fig 1], it can also be modelled as a primitive cell; the simplest unit cell, with one ion in the centre and the opposing ion on each corner giving a total of 2 atoms within the lattice, which can be repeated to give the overall crystal structure. The conventional cell has 4 times the volume of the primitive unit cell however is a more realistic representation of the symmetry of the crystal therefore is still used in many crystallographic studies. In this investigation the primitive cell will be dealt with initially when calculating using the QHA but a super cell is used for the molecular dynamics simulations which is made up of 32 conventional cells.

==Method==

All calculations were carried out using DLVisualize, a graphical interface used for modelling. Within this programme the modelling code GULP (General Utility Lattice Program) was utilised to calculate the theoretical values evaluated in the Results and Discussion below.

Initially the thermal expansion coefficient was calculated using the QHA.
The internal energy (<math>U</math>) of the primitive unit cell was initially calculated by modelling the crystal as a purely ionic material which was used in future calculations of the Helmholtz Free Energy (<math>A</math>).

To enable us to visualise how the frequencies of vibrations of this lattice vary with <math>k</math> a phonon dispersion curve was calculated using a sample of 50 k points and represented graphically [Fig ...]. Subsequently a calculation was run to find the density of states for a lattice grid size of 1x1x1 then repeated, doubling the shrinking factor each time until a grid size of 64x64x64 was obtained.

Using the optimum grid size, found to be 32x32x32 from the above simulations, the Helmholtz Free Energy was found for the MgO crystal at temperatures between 0 K and 1000 K in 100 K intervals. The lattice parameters and the total volume for each cell was recorded at each temperature and plotted [Fig ....]

From the plot of cell volume with respect to temperature, the gradient of the linear section was taken with the initial volume and used in '''Eq. 1''' to give a value of the thermal expansion coefficient, <math>\alpha</math>, as <math>2.64</math>x<math>10^{-5} K^{-1}</math> for a temperature range 400 - 1000 K.

The thermal expansion coefficient was then calculated using Molecular Dynamics.

Within this computation the system has a timestep of 1 femtosecond <math>(10^{-15} s)</math> and is run for 500 steps after equilibrium has been reached at the selected temperature. The same temperature range from 0 K to 1000 K was used as above. The system used was made up of 32 MgO conventional cells. The average equilibrium volumes were recorded at each temperature for this simulation [Fig ...]. '''Eq. 1''' was used to calculate the thermal expansion coefficient as above using the gradient of the volume as a function of temperature line of best fit and the initial volume found. Giving a value of <math>\alpha</math> as <math>2.86</math>x<math>10^{-5} K^{-1}</math> with a temperature range 100 - 1000 K.

==Results and Discussion==

The internal energy per primitive cell of MgO was found to be -41.0753 eV by choosing to model MgO as an ionic crystal with the Mg ion having a charge of +2 and the O ion -2 and having the simple Buckingham repulsive potentials operating between the ions . This is also known as the ‘binding energy’ of the crystal as it is the amount of energy that must be supplied to pull apart the ions to infinite separation.

On the X axis are points of interest along the conventional k-space for example () corresponds to centre of the lattice with coordinates (0, 0, 0) and L corresponding to the corner of the primitive cell with coordinates (0.5, 0.5, 0.5).
To compute the free energy, the sum over all vibrations of the lattice is required. To ensure the correct level of accuracy a Density of States (DOS) calculation was run with increasing shrinking factors from 1x1x1 doubling up until 64x64x64.
For the original 1x1x1 DOS graph, the k-point corresponds to the point labelled L on the phonon dispersion curve at coordinates 0.5, 0.5, 0.5 which can be seen as the wavelengths detected for the DOS calculation correspond the wavelengths seen at point L on the phonon dispersion graph (wavelengths).
As the shrinking factors increase the DOS graphs are smoothing off as more phonon vibrations are taken into account. However the difference between the 32x32x32 grid and the 64x64x64 grid is minimal, with the 64x grid taking a greater time to run for little more accuaracy.
The corresponding free energy at 300K for each of these grid sizes can be seen in Table 1. The difference between 32x and 64x is ….. which is a close enough accuracy to continue the free energy calculations with 32x.

The internal energy of a crystal of MgO

'''Calculating phonon modes of MgO'''

• The density of states for 1x1x1 grid was computed for a single k-point.

[[File:DOS1ce.png|500px]]

• Can you work out which k-point by examining the dispersion curves ? look in the Log File to check.

The k-point used to calculate the Density of states was at the cell coordinate 0.5, 0.5, 0.5 which is also the k-point label 'L' seen on the phonon dispersion graph below.

[[File:Phonon_Dispersion.png|500px]]

[[File:DOS2Ce.png|500px]]

[[File:DOS4Ce.png|500px]]

[[File:DOS8.png|500px]]

[[File:DOS16Ce.png|500px]]

[[File:DOS32Ce.png|500px]]

[[File:DOS64.png|500px]]

• How does the density of states vary with grid size ?

As the grid size increases, therefore the k points analysed increases the density of states graphs vary as seen below:

• As the grid size increases more and more of the possible vibrations are sampled and more features appear - which grid size is the minimum for a reasonable approximation to the density of states ?

grid size 16x16x16

• How does the density of states computed at this optimal grid size compare to that computed for smaller and larger grids ?

Smaller grids don't have enough resolution, above 16 graph smooths off but has no further resolution.

• How is the density of states related to the dispersion curves ?

==References==

Rep:Mod:CE1014MgO

2017-02-07T17:33:01Z

Ce1014:

Rep:Mod:CE1014MgO

2017-02-07T17:32:19Z

Ce1014:

File:Supercell.png

2017-02-07T17:30:23Z

Ce1014:

Rep:Mod:CE1014MgO

2017-02-07T17:29:00Z

Ce1014:

File:Convencell.png

2017-02-07T17:28:29Z

Ce1014: