View source for Rep:Modelling MgO RS7913

== Introduction ==
<span lang="EN-US">The fact that solid materials vibrate at all temperatures, including 0 K, is an important and useful phenomenon. </span>

I<span lang="EN-US">n the field of statistical thermodynamics, ‘bulk’ thermodynamic properties of systems or materials such as free energies, heat capacities, expansion coefficients and phase transitions can be calculated from knowledge of the partition function, ‘Q’.</span><ref>P. Atkins and J. De Paula, ''Atkins' Physical Chemistry'', OUP, 2014.</ref>
[[File:Equation1 MgO RS7913.png|thumb|Equation 1|centre]]

<span lang="EN-US">This crudely corresponds to the total number of thermally accessible energy levels in a system, and can be determined via calculation using approximate models, and/or examining the system using various spectroscopic methods.</span>

<span lang="EN-US">In a binary crystal such as magnesium oxide, at temperatures below its melting point (3125 K for MgO</span><ref>W. M. Haynes, ''CRC Handbook of Chemistry and Physics'', CRC Press, Boca Raton FL, 2011.</ref><span lang="EN-US">), there is essentially only population of excited vibrational levels, since population of electronic energy levels above the ground state requires so much energy that this can be considered as negligible. There are no rotations to be considered since we are dealing with an infinite solid lattice. The partition function of the crystal, and therefore thermodynamic properties such as those mentioned above, can thus be reasonably estimated simply by calculation of the vibrational energy levels of the crystal.</span>

=== <span lang="EN-US">Vibrations In Crystals: The Harmonic and Quasi-Harmonic Approximations</span> ===
<span lang="EN-US">In simple crystal structures, vibrations can be modelled using the harmonic approximation. </span>

<span lang="EN-US">When modelling infinite systems such as crystals, the harmonic approximation assumes only nearest-neighbour harmonic vibrations, which are small relative to the interatomic separation. The Buckingham potential (a simplification of the Lennard-Jones potential) is used to model the energy of interaction between the Mg<sup>2+</sup> and O<sup>2-</sup> ions that make up the lattice. This can be expanded about the equilibrium bond distance, r<sub>0</sub>, as follows (Equation 2):</span>
[[File:Equation 2 MgO RS7913.png|centre|thumb|518x518px|Equation 2]]

<span lang="EN-US">The first two terms of this Taylor series serve as the harmonic potential, and the motion of a one-dimensional monoatomic chain can be modelled as a linear superposition of travelling waves, by applying Newton’s second law of motion as well as the Born-von Kármán periodic boundary
condition, which joins the ends of the chain. This ensures each atom has two nearest neighbours, and is appropriate provided the chain is long. This model
yields Equation 3, which can be used to calculate the angular frequency, </span>''ω<sub>k</sub>'', of <span lang="EN-US">a one-dimensional monoatomic chain</span><ref name=":0">M. T. Dove, ''Introduction to Lattice Dynamics'', Cambridge University Press, 1993.</ref><span lang="EN-US">:</span>
[[File:Equation 3 MgO RS7913.png|centre|thumb|200x200px|Equation 3]]
'''<nowiki/>'''

J corresponds to the harmonic force constant
(analogous to that in simple harmonic motion), m to the atomic mass, k is the
wave vector and a is the interatomic separation. This model can be extended to
three dimensions by considering three different wavevectors, k<sub>x</sub>, k<sub>y</sub>
and k<sub>z</sub>, for the x-, y-, and z-directions, respectively. All
combinations of allowed wavevectors in each of the directions gives rise to all
possible vibrations in a 3D crystal, and this also gives rise to degeneracy,
i.e. different vibrations having the same energy – degeneracy can be considered
by looking at a phonon density of states (Figures 2, 3, 4, 5). Of course, a
higher density at a particular frequency (''ω<sub>k</sub>'') corresponds to
greater degeneracy of vibrations of this energy. A plot of ''ω<sub>k</sub>''
against k, or some reduced wavevector component if considering a
three-dimensional system, gives a phonon dispersion, such as that calculated
for MgO (Figure 1). 

<span lang="EN-US">The simple harmonic model can explain a number of
features that only depend on the average positions of atoms in a structure and
not their motion, often with good accuracy. Calculations are carried out in
reciprocal space, by considering a reciprocal lattice vector, a*:</span>
[[File:Equation 4 MgO RS7913.png|centre|thumb|101x101px|Equation 4]]

<span lang="EN-US">This accounts for the infinite structure that can be
made up of repeating unit cells, and a* defines the dimensions of the first
Brillouin zone – this is the area over which sampling of phonons is carried out,
and is a cube in reciprocal space (k-space).</span>

<span lang="EN-US">However, this approximation cannot be used to model
the thermal expansion of a crystal since, in this approximation, the
equilibrium interatomic distance, r<sub>0</sub>, is not dependent on
temperature. Hence, whilst a number of features that only depend on the average
positions of atoms in a structure and not their motion can be explained, such
as chemical properties, optical properties and hardness, a static model cannot
explain features such as expansion coefficients, heat capacities and
conductivity.</span>

<span lang="EN-US">The quasi-harmonic approximation introduces the
dependence of equilibrium bond distance on temperature and allows for
calculations of volume-dependent thermal effects such as expansion coefficients
or heat capacities. The approximation assumes that, at each value of the
lattice constant (i.e. at each interatomic separation), the harmonic approximation
explained earlier holds.</span>

<span lang="EN-US">The quasi-harmonic approximation uses the harmonic
approximation to calculate vibrations in a crystal, but in addition finds the
structure (i.e. interatomic separation) that is of minimum the free energy at
the temperature being studied, and this structure is considered as the
equilibrium structure. The Helmholtz free energy, being pressure-independent, is
calculated using Equation 5</span><ref name=":0" /><span lang="EN-US">:</span>
[[File:Equation 5 MgO RS7913.png|centre|thumb|373x373px|Equation 5]]

<span lang="EN-US">E<sub>0</sub> is the internal energy contribution, and
the second term is the zero-point energy; these terms can be considered as the
‘U’ term in the </span>‘A = U – TS’ definition of Helmholtz free
energy<span lang="EN-US">. The third term is a statistical thermodynamic term
that corresponds to the vibrational contribution to the entropy of the system,
and can be thought of as the ‘-T</span>S’ term. This term uses the vibrational
partition function (Q<sup>v</sup>) to calculate the entropy.

The harmonic and
quasi-harmonic approximations work especially well at lower temperatures since
they take into account the zero-point vibrational energy of the system being
studied. However, these models do not take into account the intrinsic
temperature-dependence of phonon frequencies; the relationship between energy
and internuclear separation (for the vibrations) is considered as harmonic i.e.
parabolic. Hence, due to neglecting anharmonicity, these models break down at
higher temperatures and cannot be used to model phase transitions since the parabolic
function does not allow for bond dissociation. The molecular dynamics
simulation method works better at higher temperatures, and can sometimes be
used to model phase transitions with reasonable accuracy.    

=== Classical Simulation: Molecular Dynamics ===
Molecular dynamics methods basically simulate the motion of
atoms in a real system. These simulations work by taking the ideal structure of
a system as the initial configuration, and randomly assigning velocities to all
of the atoms such that the overall kinetic energy within the system corresponds
to the target temperature being studied. Newton’s second law is then used to
calculate the forces on the atoms, and their accelerations. After a certain
time dt, known as the timestep, the velocities and positions of the atoms are
updated, and this process is repeated until average properties such as the internal
energy and temperature of the system settle. At this point, certain properties
of the system such as volume and free energy can be deduced. The Buckingham
potential is used to model interatomic interactions in the system, being the
same as that used in the quasi-harmonic approximation.

Unfortunately, molecular dynamics simulations are
considerably more computationally demanding compared to use of the
quasi-harmonic approximation, and a compromise must be struck between the size
of the system to be simulated, the length of the timestep, and the
computational power required to carry out the simulation. Of course, a larger
system size and shorter timestep will result in a more realistic and more accurate
simulation of a real infinite structure. However, this method does take the
anharmonicity of real bonds into account and, as noted earlier, is therefore
capable of producing some reasonably accurate results at higher temperature.

Due to the restricted computational power available for this
investigation, the ‘Include Shells’ option was not selected, which means a ‘rigid-ion’
model was used. This basically treats the ions as non-deformable, charged
spheres, and means that the simulation is not as accurate as ions indeed can of
course be deformed to some degree. Also, whilst this method does take
anharmonicity into account, it does not account for the zero-point energy of
the structure, since it is a classical simulation. Therefore, this method tends
to break down at lower temperatures, and cannot be used to model a system at 0
K.

=== The Coefficient Of Thermal Expansion ===
As can be seen in Equation 6, the coefficient of thermal
expansion is a measure of the rate at which a material expands with respect to temperature,
at constant pressure. It is an intrinsic property of a material that is
normalised by dividing by the initial volume, V<sub>0</sub>. The physical
origin of thermal expansion stems from the fact that temperature corresponds to
the average atomic, or molecular, kinetic energy of a substance. As temperature
is increased, more vibrations are excited and the average interatomic spacing
will increase due to an increase in the occupancy of excited vibrational states.
This can be approached from another angle, by realising that an increase in
interatomic separation (i.e. volume) will result in lower frequency vibrations,
which is entropically advantageous due to smaller separation between
vibrational energy levels leading to an increased population of excited states,
and hence an increase in the vibrational partition function, Q<sup>v</sup>.
Since the ‘-TS’ entropic term in the calculation of free energy becomes more
influential at higher temperature, it makes sense that an increase in
temperature would result in an increase in volume, due to an increase in the
entropy of the system being favoured.
[[File:Equation 6 MgO RS7913.png|centre|thumb|120x120px|Equation 6]]

=== Magnesium Oxide Structure ===
MgO takes up a face-centred cubic (fcc) crystal structure in the solid state. The conventional unit cell is a cube, containing 4 MgO units, and is most often considered as opposed to the primitive unit cell since it better reflects the symmetry of the lattice.
{|style="margin: 0 auto;"
| [[File:CONVENTIONAL CELL RS7913 MGO.gif|left|thumb|300x400px|MgO Conventional Unit Cell]]
| [[File:CONVENTIONAL CELL MGO RS7913.png|thumb|780x650px|MgO Primitive Unit Cell]]
|}

== Experimental ==
<span lang="EN-US">All calculations were carried out using the
DLVisualise GUI on a virtual (VMware Workstation on Windows 10) Linux operating
system. The GULP (General Utility Lattice Program) code was used to run
molecular simulations, and the ‘ionic.lib’ force field was selected for all
calculations in the ‘GULP Potential Parameters’ panel, with the ‘Include
Shells’ option not selected.</span>

== Results & Discussion ==

=== <span lang="EN-US">Calculation of the phonon modes of MgO</span> ===
<span lang="EN-US">The ‘Phonon Dispersion’ method was run from the
execute GULP panel, with Npoints set to 50. This corresponds to the phonons
being computed at 50 points in k-space. A 2D view of the resultant plot is
shown in Figure 1. In addition, a ‘Phonon DOS’ (density of states) was
calculated for a 1x1x1 grid in k-space i.e. a single k-point, and is featured
in Figure 2.
{|style="margin: 0 auto;"
| [[File:MgO Phonon Dispersion Plot RS7913.PNG|left|thumb|500x500px|Figure 1: Phonon dispersion for MgO.]]
| [[File:Exercise 2 - Phonon Density of States 32x32x32 Plot RS7913.PNG|thumb|502x502px|Figure 2: MgO phonon DOS, calculated for a 1x1x1 grid.]]
|}
<span lang="EN-US">Examination of Figure 2, noting the two peaks
(acoustic modes) at around 295 and 350, and two more peaks (optical modes) at
around 685 and 805 that are half of the intensity, it is clear that the Phonon
DOS was calculated for the k-point corresponding to L in the phonon dispersion
shown in Figure 1. The fact that the intensity of the optical modes is double
that of the acoustic modes indicates that the optical modes are doubly
degenerate relative to the acoustic modes. The double degeneracy of the
acoustic modes is reflected in the phonon distribution (Figure 2), since each
of the dispersion curves that cross L at 295 and 350 cm<sup>-1</sup> both branch
off into two separate curves to the left of L, whereas the other two curves do
not branch off either side of L. By examination of the output ‘Log’ file for
the phonon dispersion carried out, one can confirm that the frequencies at
point L match the peaks in the phonon DOS, and that this k-point corresponding
to the point L is at k<sub>x</sub>=k<sub>y</sub>=k<sub>z</sub>=0.5.</span>

{|style="margin: 0 auto;"
| [[File:Exercise 2 - Phonon Density of States 8x8x8 Plot RS7913.PNG|left|thumb|500x500px|Figure 3: MgO phonon DOS, calculated for a 8x8x8 grid.]]
| [[File:Exercise 2 - Phonon Density of States 32x32x32 Plot RS7913.PNG|thumb|502x502px|Figure 4: MgO phonon DOS, calculated for a 32x32x32 grid.]]
|}
[[File:Exercise 2 - Phonon Density of States 50x50x50 Plot RS7913.PNG|none|thumb|511x511px|Figure 5: MgO phonon DOS calculated for a 50x50x50 grid.]]

<span lang="EN-US">The ‘Phonon DOS’ method was run in the same way as
earlier for cubic grid sizes between 2x2x2 and 8x8x8 (Figure 3) inclusive by sequentially increasing the shrinking
factors by 1 unit. Densities of states for larger cubic grid sizes were
calculated by increasing the shrinking factors in larger steps up until a grid
size of 32x32x32 (Figure 4), and a
phonon DOS was also calculated for a grid size of 50x50x50 (Figure 5). Increasing the grid size for
which the phonon DOS is calculated for resulted in a more accurate DOS being calculated
due to a larger number of k-points, i.e. different vibrations, being sampled.
The number of k-points the DOS is computed for is equal to the ‘volume’ of the
grid, i.e. a 32x32x32 grid size corresponds to 32,768 k-points being sampled.
Increasing the grid size resulted in a smoother, more continuous curve, and it
was found that a grid size of 8x8x8 gave a continuous curve, however a grid
size of 32x32x32 gave a reasonable approximation to the density of states. It
can be seen that there is an insignificant difference between the plots of the
density of states for the 32x32x32 and 50x50x50 grid sizes, therefore using a
32x32x32 grid size is a reasonable compromise between computational power (i.e.
the length of time required to carry out a calculation) and accuracy. The
slight increase in accuracy on computing densities of states for larger grid
sizes would not be worth the additional time required to carry out the
calculation for the setup used in this investigation. Of course, this may not
be the case if a very powerful computer were used to carry out the simulations.</span>
[[File:Free energy vs Grid size RS7913.png|none|thumb|560x560px|Figure 6: The variation of Helmholtz free energy calculated with grid size used. All calculations were carried out at a temperature of 300 K.]]
{| class="wikitable"
!Grid Size
!Helmholtz Free Energy / eV
!Difference Vs Converged Value
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|0.003818
|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|0.000126
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|0.000051
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|0.000033
|-
|5x5x5
|<nowiki>-40.926463</nowiki>
|0.000020
|-
|6x6x6
|<nowiki>-40.926471</nowiki>
|0.000012
|-
|7x7x7
|<nowiki>-40.926475</nowiki>
|0.000008
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|0.000005
|-
|12x12x12
|<nowiki>-40.926481</nowiki>
|0.000002
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|0.000001
|-
|22x22x22
|<nowiki>-40.926483</nowiki>
|0
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0
|-
|100x100x100
|<nowiki>-40.926483</nowiki>
|0
|}
'''Table 1: Calculated free energy vs grid size.'''

<span lang="EN-US">Increasing the grid size resulted in convergence in
the free energy calculated as can be seen in Figure 6 and Table 1. As the grid
size was increased from 2x2x2, the free energy calculated decreased until it
converged for a grid size of 20x20x20. A larger grid size corresponds to a
larger number of k-points being sampled, just as for the phonon densities of
states calculated above, and therefore the possibility of calculating a more
accurate vibrational partition function, Q<sup>v</sup>. As explained in the
introduction, Equation 5 can be used to calculate the Helmholtz free energy,
and a sum over all ‘k’ values is required to calculate the zero-point energy
and the vibrational contribution to the free energy of the system. Therefore,
the more k-values sampled, the more accurate the calculation. Of course, as k-points
corresponding to higher and higher coordinates in k-space (i.e. higher
k-values) correspond to phonons of higher energy, their contribution to the
free energy of the system becomes less important since population of these
higher vibrational levels is less likely. This is why the free energy
calculated converges as grid size is increased.</span>

<span lang="EN-US">As can be seen in the ‘Difference Vs Converged Value’
column in Table 1, the grid size 2x2x2 would be appropriate for calculations
accurate to 1 meV and 0.5 meV, and the 3x3x3 grid size would be appropriate for
calculations accurate to 0.1 meV. It is also noted that the free energy
calculated converges to </span>-40.926483 eV for a grid size of 22x22x22, with
the accuracy<span lang="EN-US"> of 0.000001 eV being the upper limit for the software
used. A further calculation was carried out using a grid size of 100x100x100 to
ensure true convergence had been achieved. Even though a grid size of 22x22x22
would in theory be large enough for further calculations of free energy, a grid
size of 32x32x32 was chosen to ensure accuracy, since the phonon DOS plot
calculated using this grid size also gave a reasonable approximation to the
density of states, as can be seen in Figure 4. Calculations carried out using
this grid size still took under two minutes, hence this was considered a good
compromise between accuracy and time efficiency.</span>

=== <span lang="EN-US">Optimal Grid Size Vs Solid Structure For The Quasi-Harmonic Approximation</span> ===
<span lang="EN-US">The optimal grid size to be used will vary with the
solid structure being modelled and, namely, the interatomic spacing in this
solid structure. This corresponds to the lattice constant, a, of the structure,
which in turn corresponds to the size of the first Bruillon zone in reciprocal
space. As noted in the introduction, to calculate phonon densities of states
and free energies using the quasi-harmonic approximation, sampling of phonons
is done over the first Brillouin zone, which has dimensions of 2p/a. The smaller the lattice constant, a, of the structure, the larger
the size of the first Brillouin zone in k-space. Thus, a larger grid size,
corresponding to a greater number of phonons sampled within this first
Brillouin zone, would be required to gain a reasonable approximation to the DOS
and the calculation of accurate free energies. Conversely, for a structure that
has a larger lattice constant, a smaller grid size would be sufficient to gain
a reasonable approximation to the DOS and the calculation of accurate free
energies, and a smaller grid size would therefore be appropriate to ensure
computational power is not wasted.</span>

<span lang="EN-US">CaO has a fcc crystal structure that is analogous to
that of MgO, since Ca is directly below Mg in the periodic table and therefore
has very similar properties to Mg. However, Ca<sup>2+</sup> does have a
slightly larger ionic radius in comparison to Mg<sup>2+</sup> due to having an
extra shell of electrons, and this results in CaO having a slightly larger
lattice constant of 4.80 Å compared to 4.21 Å in the case of MgO. Whilst a
slightly smaller grid size would indeed be sufficient to model CaO, this would
not spare a great amount of computational power, thus it can be considered that
the optimal grid size of 32x32x32 would be appropriate for a calculation on
CaO, or another similar oxide.</span>

<span lang="EN-US">Zeolites, such as Faujasite, have complex porous
structures with much larger interatomic spacings. This corresponds to far
larger lattice constants in comparison to an oxide like MgO, and hence a
smaller first Brillouin zone in reciprocal space. Therefore, a smaller grid
size would be appropriate for a calculation on a Zeolite.</span>

<span lang="EN-US">The structure of metals such as lithium can be thought
about in terms of the ‘lattice of cations in a sea of delocalised electrons’
picture of metallic bonding, where favourable electrostatic forces of
attraction between the cations and electrons hold the structure together. Metals
generally take up close-packed structures with the atoms being held very close
to one another, giving rise to smaller lattice constants in comparison to an
ionic solid like MgO. Therefore, a larger grid size would be appropriate when
carrying out calculations on metals, such as lithium, to ensure the accuracy of
these calculations.</span>

=== <span lang="EN-US">Modelling The Thermal Expansion Of MgO Using The Quasi-Harmonic and Molecular Dynamics Simulations</span> ===
<span lang="EN-US">All calculations using the quasi-harmonic approximation
were carried out using the ‘Optimisation’ function in the ‘Execute GULP’ panel,
with ‘Optimise Gibbs free energy’ selected in the ‘Optimisation opts’ panel. In
the ‘General opts’ panel, the ‘Phonon DOS’ option was selected, all shrinking
factors were set to 32, and calculations were carried out for temperatures
between 0 and 2900 K, inclusive, by increasing in 100 K increments. The ‘Final
free energy’ and ‘Primitive cell volume’ in the output ‘Log’ files were
recorded for each of these calculations.</span>

{|style="margin: 0 auto;"
| [[File:Free Energy Vs Temp RS7913.png|left|thumb|500x500px|Figure 7: Plot of Helmholtz free energy against grid size as calculated for MgO, using the quasi-harmonic approximation.]]
| [[File:Lattice Constant Vs Grid Size RS7913.png|thumb|520x520px|Figure 8: Plot of lattice constant, a, against temperature, as calculated for MgO using a quasi-harmonic model.]]
|}

<span lang="EN-US">The variation of free energy against temperature for
these calculations carried out using the quasi-harmonic approximation was
plotted (Figure 7), and it can of
course be seen that the Helmholtz free energy decreases as the temperature is
increased. To begin with this decrease follows a curve, however, when the
temperature is greater than around 1000 K, this decrease essentially becomes
linear. By examination of Equation 5, and comparing this to the standard A = U
– TS definition of Helmholtz free energy, it can be seen that the ‘U’ terms
consisting of the E<sub>0</sub> and the zero-point energy term are
temperature-independent, and correspond to the Helmholtz free energy at 0 K,
i.e. the y-intercept in Figure 7.
The final term that corresponds to ‘-TS’ contains two ‘T’ terms, and has two
different kinds of temperature-dependence. The ‘T’ term before the summation
would lead to a linear decrease in the Helmholtz free energy with increasing
temperature, providing the summation – corresponding to ‘-S’ – is a negative,
temperature-independent term. At lower temperatures, the effect of increasing
the temperature resulting in an increase in the partition function is
significant, and results in the vibrational contribution to the ‘-S’ term (i.e.
the summation) becoming more negative as the temperature is increased. As the
temperature is increased to over 1000 K, the vibrational contribution to
entropy can be considered to be temperature-independent, and the
temperature-dependence of the Helmholtz free energy is essentially linear.</span>

<span lang="EN-US">The volumes in the output ‘Log’ files for these
calculations carried out correspond to the primitive unit cell. Thus, to
calculate the lattice constant, a, these volumes were multiplied by 4 to obtain
the volume of the conventional unit cell, and then the cube root was taken to
give a. A plot of lattice constant against temperature (Figure 8) shows that, as expected, an increase in temperature
results in an increase in lattice constant, i.e. greater interatomic
separation.''' '''</span>

Molecular dynamics simulations were carried out
using a supercell containing 32 MgO units. The ‘Molecular Dynamics’ method was
run from the ‘Execute GULP’ panel, with ‘ensemble’ set to ‘NPT’ in the ‘MD
Opts’ panel. The ‘Time Step’ was set to 1.0 femtoseconds, ‘Equilibration steps’ and ‘Production steps’ set to 500, and ‘Sampling steps’ and ‘Trajectory write steps’ set to 5. Calculations <span lang="EN-US">were carried out for temperatures between 100 and 3000 K, inclusive, by increasing in 100 K increments. The ‘Averaged Cell volume’ for the final timestep, 0.5 ps, in the output ‘Log’ files were recorded for each of these calculations.

The ensemble being set to ‘NPT’ corresponds to the number of particles (N), external pressure (P) and temperature (T) being fixed. Of course, the volume of the system is allowed to vary. It is essential that the correct timestep is chosen, and a compromise must be drawn between computational power and accuracy. The timestep must be long enough to ensure the calculation is efficient, but short enough to ensure all possible vibrations are reasonably well sampled; 10 steps along a given vibration is typically adequate. In this investigation the 1 femtosecond timestep gave reasonable results, and can be deemed acceptable for MgO. In addition, the size of the system to be simulated is important. The supercell used in this investigation, containing 32 MgO units, which also chosen as a compromise between computational power and accuracy of the calculations carried out. The use of a larger supercell would result in greater accuracy, since the increased flexibility would allow for the simulation of a greater number of vibrational modes. Indeed, altering the timestep and size of the supercell could certainly be a subject of further investigation, provided one had the computational resources, or time, required to carry out more complex calcuations involving a shorter timestep and/or larger supercell.[[File:Volume vs Temperature QHA and MD RS7913.png|none|thumb|600x600px|Figure 9: Plot of MgO conventional unit cell volume against temperature, for both calculations carried out using the quasi-harmonic and molecular dynamics simulations.]]The conventional unit cell volumes calculated – using both the quasi-harmonic and molecular dynamics models – at each of the temperatures studied – are plotted in Figure 9. Equation 6 was used to calculate the coefficient of thermal expansion, where dv/dt was taken as the gradient of the linear
region of the graph. By examination of Figure 9, for the quasi-harmonic
approximation, the linear region was taken as the points between 500 K and 1900
K, inclusive, and for the molecular dynamics simulation the points between 100
K and 2100 K, inclusive. The plots of these linear regions and the
corresponding linear fits can be seen in figures 10 and 11 for the
quasi-harmonic and molecular dynamics simulations, respectively. Calculation of
the expansion coefficients gave rise to values of 36.6 x 10<sup>-6</sup> K<sup>-1</sup>
and 29.8 x 10<sup>-6</sup> K<sup>-1</sup> for the quasi-harmonic approximation
and molecular dynamics simulation, respectively. The expansion coefficient does
have a considerable temperature dependence, thus the literature values which
the calculated expansion coefficients were compared to were averaged over the
temperature range being studied, and these values can be seen in Table 2.
{| style="margin: 0 auto;"
| [[File:QHA VOL VS TEMP RS7913.png|left|thumb|516x500px|Figure 10: Plot of MgO conventional unit cell volume against temperature between 500 K and 1900 K, calculated using the quasi-harmonic approximation.]]
| [[File:Conv unit cell vs temp MD RS7913.png|thumb|530x500px|Figure 11: Plot of MgO conventional unit cell volume against temperature between 100 K and 2100 K, calculated using the molecular dynamics model.]]
|}
{| class="wikitable"
!
!Quasi-Harmonic Approximation
!Molecular Dynamics Simulation
|-
|Temperature range of the linear region / K
|500-1900
|100-2100
|-
|Computed coefficient of thermal expansion / 10<sup>-6</sup> K<sup>-1</sup>
|36.6
|29.8
|-
|Literature value corresponding to average coefficient of thermal expansion over the temperature range being studied / 10<sup>-6</sup> K<sup>-1</sup>
|45.3
|42.3
|}
'''Table 2: Coefficients of thermal expansion for MgO calculated using the quasi-harmonic and molecular dynamics models, and corresponding literature'''<ref name=":1">O. L. Anderson and K. Zou, ''J. Phys. Chem. Ref. Data'', 1990, '''19''', 69-81.</ref>''' values.'''

As can be seen in Table 2, the computed values of the
coefficient of thermal expansion do agree reasonably well with the literature
values, and the quasi-harmonic approximation provides a more accurate value for
the temperature range studied. In addition, the expansion coefficient
calculated over the temperature range 500-1900 K using the quasi-harmonic
approximation was greater than that calculated over the temperature range
100-2100 K, which agrees with the cited literature values.

By examination of Figure 9,
it can be seen that the quasi-harmonic approximation best models the system at
lower temperatures. As noted in the introduction, this is because the
quasi-harmonic approximation takes into account the zero-point vibrational
energy of the system. This is neglected by the molecular dynamics model, and
results in the model predicting a volume that is excessively small. However, at
temperatures above around 2300 K, the quasi-harmonic approximation begins to
break down. As the melting point of MgO, 3125 K, is reached, the volume
computed using the quasi-harmonic approximation starts to increase at a rapidly
steeper and steeper rate. This is because the phonon modes are harmonic, and
therefore do not represent the actual motion of the ions at higher temperature,
since the population of higher vibrational levels means that anharmonic effects
and bond dissociation become important. The quasi-harmonic continues to treat
MgO as a fcc lattice even at these higher temperatures, and this ultimately
results in an error when attempting calculations for temperatures above 3000 K.
It can be seen that, on the other hand, the molecular dynamics simulation does
better model the system as the temperature is increased to near to the boiling
point of MgO. This model does not neglect anharmonicity, and the atoms are
allowed to move apart from one another since they are not restricted by a
harmonic potential. As noted in the introduction, the molecular dynamics
simulation method can sometimes be used to model phase transitions with
reasonable accuracy. However, the molecular dynamics simulations carried out in
this investigation do start to break down at temperatures above roughly 2500 K,
where the cell volume calculated starts to hover up and down. This is because the
timestep of 1.0 femtoseconds used was not short enough to model the motion of
atoms accurately at these higher temperatures, due to the greater kinetic
energy the particles have. This could in theory be resolved by using a shorter
timestep and, as mentioned in the introduction, the molecular dynamics
simulation method can sometimes be used to model phase transitions with
reasonable accuracy. Using a shorter timestep for these calculations could be a subject for further research.        

== Conclusion ==
An investigation was carried out to compare the quasi-harmonic approximation and molecular dynamics model in approximating the thermal expansion of magnesium oxide. To begin with, it was discovered that the optimal grid size to use for quasi-harmonic calculations in this investigation was 32x32x32, since this grid size gave rise to a sufficiently realistic phonon DOS as well as convergence of the Helmholtz free energy calculated to the maximum level of precision possible using the GULP program used. This grid size was used for all subsequent calculations.

The quasi-harmonic approximation was used to calculate the Helmholtz free energy of MgO at temperatures between 0 K and 2900 K, inclusive, and it was found that an increase in temperature resulted in a decrease in the free energy. 

Finally, the thermal expansion of MgO was modelled using both the simple harmonic approximation and a molecular dynamics simulation over a temperature range of 0 K and 3000 K. It was found that the quasi-harmonic approximation was good at modelling MgO at low temperatures since it takes the zero-point vibrational energy of the system into account. However, the approximation broke down at higher temperatures since it does not take anharmonic effects into account, a calculation at 3000 K resulted in an error output file. The molecular dynamics method did not model MgO well at lower temperatures since this approximation does not take the zero-point vibrational energy of the system into account, and thus a calculation at 0 K resulted in an error output file. It did fare better than the quasi-harmonic approximation at higher temperatures, however the molecular dynamics model did start to break down at temperatures approaching 3000 K. This is likely to be due to the 1.0 femtosecond timestep used not being sufficiently short to model the motion of atoms with such high kinetic energies. The coefficient of thermal expansion of MgO was calculated using the quasi-harmonic and molecular dynamics models, over the temperature range 500-1900 K and 100-2100 K, respectively. This yielded values of 36.6 x 10<sup>-6</sup> K<sup>-1</sup> and 29.8 x 10<sup>-6</sup> K<sup>-1</sup> for the quasi-harmonic and molecular dynamics approximations, which correlate reasonably well with literature values<ref name=":1" /> of 45.3 x 10<sup>-6</sup> K<sup>-1</sup> and 42.3 x 10<sup>-6</sup> K<sup>-1</sup> for the corresponding temperature ranges. Discrepancy in the results in comparison to the literature could in theory be minimised by using a larger grid size and including shell modes in the quasi-harmonic model, and by using a larger supercell and shorter timestep for molecular dynamics simulations. To determine whether this is indeed the case could be the subject of further investigation. However, one must not forget that these computational models provide a convenient and fast way of calculating thermodynamic properties such as free energy and volume without having to carry out experiments on the systems in reality, in the laboratory. Trying to gain more accurate results using computational methods will certainly require more computational power, and this requirement may finally begin to take away from some of the advantages of using a computational method.  

==References==
-