== Abstract ==
== Introduction ==
Thermal properties of materials are described by their vibrational free energies, which can be described in terms of the relative motion of atoms or the motion of their centre-of-mass. <ref name='phonons'>G. Srivastava, The physics of phonons, A. Hilger, Bristol, 1990.</ref> These concepts give rise to different approaches in calculating vibrational free energy, and both methodologies will be explored in greater detail.

=== Quasi-harmonic Approximation (QHA) ===
Fundamentally, QHA invokes the description of a crystalline solid as a primitive unit cell. This is essential due to the impracticality of calculating all the vibrational degrees of freedom in a crystal—for a crystal of size <math>N</math>, there are <math>3N</math> degrees of vibrational freedom, and in an infinitely large crystal lattice, 3N --> infinity.
Nonetheless, the translational periodicity of the crystal lattice, where <math display="inline">f(r + T) = f(r)</math>, simplifies the dynamics of all atoms in the lattice into that of a unit cell. For such a simplification to be appropriate, the following assumptions are made.

==== The Adiabatic Approximation ====
The adiabatic approximation separates the motion of the ion cores from that of the electrons since former are much more massive than the latter. Hence, the ion cores can be assumed to be in their equilibrium positions and that their motion is dependent on the potential field generated from the average motion of electrons.<ref name='phonons' />

==== The Harmonic Approximation ====
The total potential energy of a crystal can be expressed as the sum of all interatomic potentials. A two-body system typically has an anharmonic potential energy surface <math display="inline">U\left( r\right) </math>, where <math display="inline">r</math> is the interatomic separation. By considering a small displacement <math display="inline">x</math>,
<math display="block">x = r-r_0 \left( 1 \right)</math>

Where <math display="inline">r_0</math> is the equilibrium distance between the first and second atoms and is a minimum on the potential energy surface.

<math display="inline">U\left( r\right) </math> can be expanded in a Taylor series about <math display="inline">r_0</math>.
<math display="block">U\left( r\right) = U(r_0) + \frac{\partial U}{\partial x} x + \frac{\partial^2 U}{\partial x^2} x^2 + ... </math>
Since <math display="inline">U(r_0) </math> is unimportant in dynamics, <math display = "inline">\frac{\partial U}{\partial x}</math> is a force term and must be 0 for an equilibrium configuration, and all higher order terms <math display = "inline">x_n</math> , where <math display = "inline">n \ge 3</math> are assumed to be close to 0. As such, only the quadratic term is considered in the harmonic approximation. The solutions are the normal modes of vibrations for a system of independent quantum oscillators.

A phonon is a quantum of vibrational energy, hw, associated with a wave vector k.

Hence, for a crystal, its potential energy is given in the following equation.
[First principles, pdf, page 2]
Where l and k are the labels of the unit cells and atoms in each unit cell respectively. (http://dx.doi.org/10.1016/j.scriptamat.2015.07.021)

===== Limitations of Harmonic Approximation =====
The harmonic approximation predicts symmetric atomic vibrations about r0 at all temperatures, and is therefore incongruent with observed phenomena such as thermal expansion and heat conductivity.<ref>G. Peckham, PhD, Trinity College, Cambridge, 1964.</ref> The QHA causes renormalisation of the phonon frequencies and atomic force constants as is appropriate for the thermal equation of state.<ref>G. Leibfried and W. Ludwig, Solid State Physics, 1961, 275-444.</ref>

=== MD Simulation ===
MD considers the forces exerted on each atom and provides a classical description of an <math display="inline">N</math>-atom system. This is given by
<math display="block">M \left( \frac{\partial^2 \mathbf{r_i}}{\partial t^2}\right) = \sum_{j=1, j \ne i}^N \mathbf{F_{ij}} </math>

== Methodology ==
Unless otherwise stated, all calculations were performed on a primitive unit cell of MgO with lattice parameters <math display="inline">a = 2.9783 \AA</math>, <math display = "inline">a = b = c</math>; <math display="inline">\alpha = 60</math>°, and <math display="inline">\alpha = \beta = \gamma </math> with GULP version 1.4.43 and crystals visualised with DLV interface.

A phonon dispersion curve was computed by sampling 100 points within the first Brillouin zone. The phonon density of states (DOS) was calculated with various shrinking factors, and the graphs subsequently plotted with matplotlib. The free energy of MgO was calculated with different shrinking factors at 300 K, and a suitable shrinking factor selected for the subsequent investigation of the thermal expansion of MgO. For every run, the Gibbs free energy was optimised, and calculations were performed from 0 to 2960 K in temperature steps of 20 K.

All MD simulations were performed on an isothermal-isobaric ensemble of MgO supercell of 32 formula units, with the following cell parameters:
<math display = "block">a = 8.4239 \AA</math>
<math display = "block">a = b = c</math>
<math display = "block">\alpha = 90^o</math>
<math display = "block>\alpha = \beta = \gamma</math>

MD was performed over a temperature range of 20 K to 4000 K, with temperature steps of 20 K. All calculations were performed with a time step of 1 fs. From 20 K to 1680 K, the system was allowed to first equilibrate for 1 ps; this was increased to 5 ps from 1700 K to 4000 K. Following which, MD production was allowed to run for 5 ps for all temperatures.

All data was analysed with Python on Jupyter notebook, and all graphs plotted with matplotlib.

== Results and Discussion ==
The lattice energy of MgO was calculated to be -41.0753 eV per primitive unit cell.
=== Phonon Modes of MgO ===
Figure 1 illustrates the phonon dispersion curve computed at 100 points for the primitive unit cell.

[[File:Syl815Phonon_Dispersion_Graph_100.PNG|thumb|center|600px|'''Figure 1''' Phonon Dispersion Curve of MgO]]

A salient feature is the presence of 6 branches in the dispersion diagram. Assuming that the Born-von Karman boundary condition is satisfied, the edge effects of cells on dynamics can be ignored and <math display = "inline">u_{N+1} = u_1</math>, where <math display = "inline">u</math> is the displacement and <math display = "inline">N</math> is the number of unit cells. This also implies the translational symmetry in <math display = "inline">k</math>-space, such that all information of phonon dispersion can be derived by sampling in the first Brillouin zone (FBZ).

By considering a linear diatomic chain satisfying the periodic boundary condition, the solutions to the vibrational frequency can be expressed in the form
<math display = "block">\omega^2 = \Lambda \left( \frac{1}{m} + \frac{1}{M} \right) \pm \left[ \left( \frac{1}{m} + \frac{1}{M} \right) ^2 - \frac{4}{mM} \sin^2 ka \right]^{\frac{1}{2}}</math>

where <math display = "inline">\omega</math> is the angular frequency, <math display = "inline">\Lambda</math> is the force constant of the bond, <math display = "inline">\left( \frac{1}{m} + \frac{1}{M} \right)</math> is the reduced mass of the system, where <math display = "inline">m</math> is the mass of the lighter atom (O) and <math display = "inline">M</math> is the mass of the more massive atom (Mg), and <math display = "inline">a</math> is the length of the unit cell.

The equation highlights two possible solutions for each <math display = "inline">k</math>-value in a linear chain. Moreover, when <math display = "inline">m \ne M</math>, a gap is observed at<math display = "inline">k = \frac{\pi}{2a}</math>, which is observed in Figure 1.<ref>R. Hornreich, M. Kugler, S. Shtrikman and C. Sommers, Journal de Physique I, 1997, 7, 509-519.</ref>

By extending the logic to a 3D crystal lattice, the number of branches observed is given by <math display = "inline">3x</math>, where <math display = "inline">x</math> is the number of atoms per unit cell. This is in agreement with the observation in Figure 1.

By appraising the solutions for <math display = "inline">k=0</math> (long wavelength limit),

<math display = "block">\omega_1 = 2\Lambda \left( \frac{1}{m} + \frac{1}{M} \right)</math>

<math display = "block">\omega_2 = 0 </math>

<math display = "inline">\omega_1</math> corresponds to a high energy mode in which the atoms in the unit cell are moving out-of-phase, where frequency values are within the visible electromagnetic spectrum. The atoms are able to interact with an electric field of appropriate frequency due to the presence of both a positive and negative charge within the unit cell. It is hence naturally termed the optical mode.<ref>M. Dove, Introduction to Lattice Dynamics, Cambridge University Press, Cambridge, 1993.</ref>

On the other hand, <math display = "inline">\omega_2</math> corresponds to a low energy mode with the atoms moving in phase and the wave pattern is similar to sound waves—hence the term acoustic mode. For any crystal with <math display="inline">N</math> atoms in the unit cell, there are only 3 acoustic—2 transverse and 1 longitudinal—and <math display="inline">3N-3</math> optical branches. The transverse modes are perpendicular to <math display="inline">k</math>, while the longitudinal mode is parallel.

=== Computing Density of States (DOS) ===
The impracticality of sampling all <math display="inline">k</math>-points within the FBZ can be circumvented by the use of a commensurate grid of <math display="inline">k</math>-points. To determine this set of <math display="inline">k</math>-points, the Pack-Monkhorst (PM) shrinking factor was used to specify the number of equidistant <math display="inline">k</math>-points taken along each direction of <math display="inline">b_1</math>, <math display="inline">b_2</math> and <math display="inline">b_3</math> in one reciprocal lattice PUC.<ref>A. Parrill and K. Lipkowitz, Reviews in Computational Chemistry, Volume 29, John Wiley & Sons, 2016.</ref> The Cartesian coordinates of the <math display="inline">k</math>-points calculated are given by the equation

[EQUATION FROM https://journals.aps.org/prb/pdf/10.1103/PhysRevB.93.155109]

A major advantage is its computational efficiency by restricting the number of <math display="inline">k</math>-points calculated to a finite value. Moreover, the accuracy obtained from calculations with a PUC can be comparable to that of a supercell as long as the shrinking factor is appropriate.

Table 1 illustrates the effect of modifying the PM shrinking factor on the number of k-points calculated.

'''Table 1'''. Grid size against number of <math display="inline">k</math>-points
{| class="wikitable" style="text-align: center; width: 85%;margin: auto;"
! Grid Size (n x n x n)
! Number of k-points
|-
| 1
| 1
|-
| 2
| 4
|-
| 3
| 18
|-
| 4
| 32
|-
| 5
| 75
|-
| 6
| 108
|-
| 8
| 256
|-
| 10
| 500
|-
| 16
| 2048
|-
| 20
| 4000
|-
| 64
| >99 999
|}

As the mesh of <math display="inline">k</math>-points increases, the number of <math display="inline">k</math>-points calculated increases as well. This is contrary to the prediction from the above equation, where we would expect <math display="inline">k_x \times k_y \times k_z</math> number of points. This can be attributed to the mapping of equivalent <math display="inline">k</math>-points onto each other and thus the number of <math display="inline">k</math>-points calculated is reduced.

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of <math display="inline">k</math>-points is increased.

[[File:Syl815DOS.png|thumb|center|1000px|DOS plots against wavenumber when PM shrinking factor was varied.]]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication, and will be further discussed in section 3.

An initial plot of the density of states was obtained from a 1x1x1 grid yielding six resultant modes. Sharp and distinct peaks are observed in the plot, since only one k-point was sampled.

Notably, only four unique peaks are observed even though we should observe 6 modes of vibrations. The final two modes are degenerate at _____ and _____ wavenumbers. Compared to the non-degenerate acoustic and optical peaks (___ and ____ respectively), the degenerate acoustic modes are higher in energy whereas the degenerate optical modes are lower in energy correspondingly. It can therefore be deduced that the degenerate acoustic and optical modes are transverse in nature.

The <math display="inline">k</math>-point used in the DOS calculation could be identified by comparing with the dispersion curve. Since point M contains all of the frequency values in Figure ___, it can be determined that the point represented in the DOS curve is M, where <math display="inline">k_x</math> = 0.5, <math display="inline">k_y</math> = 0.5 and <math display="inline">k_z</math> = 0.5.

=== Relationship between the Dispersion Curve and DOS ===
The DOS curve illustrates the number of energy states per unit energy, demonstrating a mode at <math display = "inline">414 cm^{-1}</math>. This correlates well with Figure ____. By constructing a horizontal line at frequency = 414 cm-1, it can be observed that the branches intersect this line frequently. This implies that a significant proportion of k-points have vibrational modes of frequency <math display = "inline">414 cm^{-1}</math>. The DOS curve can thus be interpreted as the orthogonal of the dispersion curve.

The dispersion diagram is useful in locating the band gaps of the acoustic and optical modes - for electronic dispersion diagram, this is useful in identifying whether a material has a direct bandgap or an indirect one, which affects the properties of the material and its use.

However, the dispersion diagram only illustrates the energy values calculated at the special points chosen, interpolating the energies of the vibrational modes for the k-points which are not calculated. The DOS plot is in this respect more meaningful, the energy states for all of <math display = "inline">k</math> values are accounted in this representation.

=== Computing the Free Energy Using the Harmonic Approximation ===
The figure below demonstrates the relationship between the PM shrinking factor used and the computed Helmholtz free energy of the system.

[[File:Syl815FEshinking.png|thumb|center|600px|Free energy vs. PM shrinking factor]]

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[[file:Syl815FEshinking.png|thumb|center|600px|DOS against wavenumber for various shrinking factors.]]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication.

From the above figure, the free energy of MgO is observed to increase and converge to a value of -40.926 483 eV, and it is observed that this occurs for a grid size of 8x8x8.

A 2x2x2 grid is sufficient for calculating the free energy of MgO to 1 meV. A 4x4x4 grid is necessary for a precision to 0.5 meV and 0.1 meV.

=== Thermal Expansion ===
The Helmholtz free energy of a crystal is given by the sum of the energies of independent vibrational waves. The energy level <math display="inline">n</math> of a quantum harmonic oscillator are given by:

<math display="block">E_n = \left( n+ \frac{1}{2} \right) h \nu </math>

where <math display="inline">h</math> is Planck's constant and <math display="inline">\nu</math> is the frequency of energy level <math display="inline">n</math>. For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators, the vibrational energy is given by:

<math display="block">E_{vib} = \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu</math>

For a canonical <math display="inline">NVT</math> ensemble, the partition function is given by

<math display="block">Z = \sum_n e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

where <math display="inline">\beta = \frac{1}{kT}</math> and <math display="inline">E_n</math> enumerates all vibrational energy states.

For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators,

<math display="block">Z_N = \prod_n^{3N} e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

The phonon entropy can then be expressed in terms of the partition function:
<math display="block">S=-k_B \ln Z_N</math> where <math display="inline">k_B</math> is the Boltzmann constant.

Given the relation
<math display="block">F=U+TS</math>
where <math display="inline">U</math> is the internal energy of the system— for a crystal this is its electric potential energy
<math display="block">U_E = \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}}</math>

where <math display="inline">i</math> and <math display="inline">j</math> are the indices of the ions, <math display="inline">r_{ij}</math> is the distance between <math display="inline">i</math> and <math display="inline">j</math> and <math display="inline">\varepsilon_0 = 8.8542 \times 10^{-12} F\cdot m^{-1} </math>

The Helmholtz free energy of a crystal is thus given by

<math display="block">F= \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}} + \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu + k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math>

This equation could be used to qualitatively rationalise the free energy dependence on temperature. The data obtained is plotted in Figure ____.
[FIGURE]
Particularly, there are two salient regimes of interest. At low temperatures, T < 100 K, the graph is flat. However, at high temperatures, the behaviour is approximately linear. These observations are in agreement with the above equation, which highlights the temperature dependence of entropy <math display="inline">S</math>. At low temperatures, the term <math> k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> is extremely small, and hence the free energy term is dominated by the internal energy of the crystal. At high temperatures, the term <math>-k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> dominates and therefore the free energy of the system appears to have a dependence in temperature.

==== Variation of Lattice Parameter with Temperature ====
[[File:Syl815CellPvsT.png|thumb|center|600px|This figure illustrates the variation in cell parameter of MgO with temperature]]

As the temperature increases, the lattice parameter increases. It can thus be observed that the cell volume <math display="inline">V </math> has a dependence on temperature <math display="inline">T </math>, and the thermal expansion coefficient <math display="inline">\alpha </math> is given by

<math display="block">\alpha = \frac{1}{3V} \left( \frac{\partial V}{\partial T}\right)_P = \frac{1}{3B} \left( \frac{\partial P}{\partial T}\right)_V </math>
Where <math display="inline">B </math> is the bulk modulus and <math display="inline">P </math> is the pressure.
At 300 K, <math display="inline">\alpha = 2.06\times 10^{-5} K^{-1}</math>, compared to a literature value of <math display="inline">\alpha = 3.06\times 10^{-5} K^{-1}</math>.<ref>N. Corsepius, T. DeVore, B. Reisner and D. Warnaar, Journal of Chemical Education, 2007, 84, 818</ref>

=== MD Simulation ===
The cell volume per formula unit of MgO was plotted against temperatures between 20 K to 4000 K.
[[File:Syl815MD.png|thumb|center|600px|"Experimental Data MD"]]
Under MD, the cell volume generally increases linearly with temperature throughout. By considering the mean kinetic energy of the crystal
<math display="block">\left \langle E_k \right \rangle = \frac{1}{2} M \sum_{i=1}^N v_i^2 = \frac{3}{2} Nk_BT_{MD} </math>
Where <math display="inline">\left \langle E_k \right \rangle </math> is the average kinetic energy of the atoms, <math display="inline">M</math> is the mass of the crystal lattice, and <math display="inline">v_i</math> represents the velocity of the atom <math display="inline">i</math>. It can be observed that the cell energy is linearly dependent on temperature. In a constant pressure system, this would result in volume expansion as temperature increases.

It can be observed that at high temperatures when <math display="inline">T\ge 2000 K </math>, more noise is present in the data due to the large cell volume and the large kinetic energy of the atoms.

[[File:Syl_MDvsQHA.png|thumb|center|600px|"This figure compares the data obtained for the thermal expansion of MgO under QHA and under MD."]]

At extremely low temperatures of <math display="inline">T\le 200 K </math>, QHA predicts a larger cell volume than MD. This can be attributed to the significant quantum effects at such low temperatures. Since MD only accounts for the kinetic energy of the atoms and neglects zero point vibrations, it predicts a smaller cell volume with the atoms closer together.

The data obtained for MD and QHA demonstrate strong agreement for temperatures between 200 to 1000 K. At these temperatures, the thermal energy of the system is sufficiently large such that the motion of the particles can be described classically.

== References ==

2018-02-14T11:59:36Z

Syl815: /* Phonon Modes of MgO */

== Abstract ==
== Introduction ==
Thermal properties of materials are described by their vibrational free energies, which can be described in terms of the relative motion of atoms or the motion of their centre-of-mass. <ref name='phonons'>G. Srivastava, The physics of phonons, A. Hilger, Bristol, 1990.</ref> These concepts give rise to different approaches in calculating vibrational free energy, and both methodologies will be explored in greater detail.

=== Quasi-harmonic Approximation (QHA) ===
Fundamentally, QHA invokes the description of a crystalline solid as a primitive unit cell. This is essential due to the impracticality of calculating all the vibrational degrees of freedom in a crystal—for a crystal of size <math>N</math>, there are <math>3N</math> degrees of vibrational freedom, and in an infinitely large crystal lattice, 3N --> infinity.
Nonetheless, the translational periodicity of the crystal lattice, where <math display="inline">f(r + T) = f(r)</math>, simplifies the dynamics of all atoms in the lattice into that of a unit cell. For such a simplification to be appropriate, the following assumptions are made.

==== The Adiabatic Approximation ====
The adiabatic approximation separates the motion of the ion cores from that of the electrons since former are much more massive than the latter. Hence, the ion cores can be assumed to be in their equilibrium positions and that their motion is dependent on the potential field generated from the average motion of electrons.<ref name='phonons' />

==== The Harmonic Approximation ====
The total potential energy of a crystal can be expressed as the sum of all interatomic potentials. A two-body system typically has an anharmonic potential energy surface <math display="inline">U\left( r\right) </math>, where <math display="inline">r</math> is the interatomic separation. By considering a small displacement <math display="inline">x</math>,
<math display="block">x = r-r_0 \left( 1 \right)</math>

Where <math display="inline">r_0</math> is the equilibrium distance between the first and second atoms and is a minimum on the potential energy surface.

<math display="inline">U\left( r\right) </math> can be expanded in a Taylor series about <math display="inline">r_0</math>.
<math display="block">U\left( r\right) = U(r_0) + \frac{\partial U}{\partial x} x + \frac{\partial^2 U}{\partial x^2} x^2 + ... </math>
Since <math display="inline">U(r_0) </math> is unimportant in dynamics, <math display = "inline">\frac{\partial U}{\partial x}</math> is a force term and must be 0 for an equilibrium configuration, and all higher order terms <math display = "inline">x_n</math> , where <math display = "inline">n \ge 3</math> are assumed to be close to 0. As such, only the quadratic term is considered in the harmonic approximation. The solutions are the normal modes of vibrations for a system of independent quantum oscillators.

A phonon is a quantum of vibrational energy, hw, associated with a wave vector k.

Hence, for a crystal, its potential energy is given in the following equation.
[First principles, pdf, page 2]
Where l and k are the labels of the unit cells and atoms in each unit cell respectively. (http://dx.doi.org/10.1016/j.scriptamat.2015.07.021)

===== Limitations of Harmonic Approximation =====
The harmonic approximation predicts symmetric atomic vibrations about r0 at all temperatures, and is therefore incongruent with observed phenomena such as thermal expansion and heat conductivity. (https://core.ac.uk/download/pdf/42335965.pdf) The QHA causes renormalisation of the phonon frequencies and atomic force constants as is appropriate for the thermal equation of state. (https://ac.els-cdn.com/S0081194708606566/1-s2.0-S0081194708606566-main.pdf?_tid=f9617aac-104c-11e8-87bc-00000aab0f27&acdnat=1518478522_8ee2607e5ac970b2f00844b9ef59e8dc)

=== MD Simulation ===

=== System considered ===
MgO has a face-centred cubic lattice; more specifically, an alkali halide type structure with a Fm3m space group. (Citation)
MgO: insulator
Lattice parameters, symmetry, space group, etc.

=== Aims of the Exercise ===

== Methodology ==
Unless otherwise stated, all calculations were performed on a primitive unit cell of MgO with lattice parameters a = 2.9783 Angstrom, a = b = c; alpha = 60 degrees, and alpha = beta = gamma with GULP version 1.4.43 and crystals visualised with DLV interface.

A phonon dispersion curve was computed by sampling 100 points within the first Brillouin zone. The phonon density of states (DOS) was calculated with various shrinking factors, and the graphs subsequently plotted with matplotlib. The free energy of MgO was calculated with different shrinking factors at 300 K, and a suitable shrinking factor selected for the subsequent investigation of the thermal expansion of MgO. For every run, the Gibbs free energy was optimised, and calculations were performed from 0 to 2960 K in temperature steps of 20 K.

All MD simulations were performed on an isothermal-isobaric ensemble of MgO supercell of 32 formula units, with the following cell parameters:
A = 8.4239 Angstrom
A = b = c
Alpha = 90 degrees
Alpha = beta = gamma

MD was performed over a temperature range of 20 K to 4000 K, with temperature steps of 20 K. All calculations were performed with a time step of 1 fs. From 20 K to 1680 K, the system was allowed to first equilibrate for 1 ps; this was increased to 5 ps from 1700 K to 4000 K. Following which, MD production was allowed to run for 5 ps for all temperatures.

All data was analysed with Python on Jupyter notebook, and all graphs plotted with matplotlib.

== Results and Discussion ==
The lattice energy of MgO was calculated to be -41.0753 eV per primitive unit cell.
=== Phonon Modes of MgO ===
Figure 1 illustrates the phonon dispersion curve computed at 100 points for the primitive unit cell.

[[File:Syl815Phonon_Dispersion_Graph_100.PNG|thumb|center|600px|Phonon Dispersion Curve.]]

A salient feature is the presence of 6 branches in the dispersion diagram. Assuming that the Born-von Karman boundary condition is satisfied, the edge effects of cells on dynamics can be ignored and u(N+1) = u1, where u is the displacement and N is the number of unit cells. This also implies the translational symmetry in k-space, such that all information of phonon dispersion can be derived by sampling in the first Brillouin zone.

By considering a linear diatomic chain satisfying the periodic boundary condition, the solutions to the vibrational frequency can be expressed in the form
[EQUATION, pg. 27 of book]

Highlighting two possible solutions for each k-value in a linear chain. Moreover, when m1 /= m2, a gap is observed at k = pi/2a, which is observed in Figure 1. (https://hal.archives-ouvertes.fr/file/index/docid/247340/filename/ajp-jp1v7p509.pdf)

By extending the logic to a 3D crystal lattice, the number of branches observed is given by 3x, where x is the number of atoms per unit cell. This is in agreement with the observation in Figure 1.

By appraising the solutions for k = 0 (long wavelength limit),

w1 = 2lambda(1/m + 1/M)
w2 = 0

w1 corresponds to a high energy mode where the atoms in the unit cell are moving out-of-phase, where frequency values are within the visible electromagnetic spectrum. The atoms are able to interact with an electric field of appropriate frequency due to the presence of both a positive and negative charge within the unit cell. It is hence naturally termed the optical mode. (Introduction to lattice dynamics)

On the other hand, w2 corresponds to a low energy mode with the atoms moving in phase and the wave pattern is similar to sound waves—hence the term acoustic mode. For any crystal with N atoms in the unit cell, there are only 3 acoustic—2 transverse and 1 longitudinal—and 3N-3 optical branches. The transverse modes are perpendicular to k, while the longitudinal mode is parallel.

For a cubic crystal, the highly symmetric nature indicates the possibility for some vibrations to be degenerate. Typically, transverse modes are lower in energy due to the weaker interaction between atoms in the unit cell.

This is encapsulated in Bloch's theorem and thus the wavefunction psi can be expressed in the form as given in equation 1.

[EQUATION HERE]

While typically applied to electrons in crystals, Bloch's theorem is broadly applicable in describing periodic wave phenomena, such as in phononic crystals. The branches relate to the modes of vibration, as given by equation for a one-dimensional system.

[EQUATION HERE]

Since the system is a three-dimensional one, vibrations can occur independently in perpendicular planes, therefore giving rise to additional normal coordinates and vibrations.

A dispersion diagram enables clear visualisation of the nature of the band gap, as illustrated in Figure ___. Figure __ clearly illustrates acoustic and optical phonons, depicting the in-phase and out-of-phase movement of the ions respectively.

Nonetheless, it is limited since energy values are unknown for points which are not sampled. Calculation of the density of states has been proven to be more useful.

=== Computing Density of States (DOS) ===
The impracticality of sampling all k-points within the FBZ can be circumvented by the use of a commensurate grid of k-points. To determine this set of k-points, the Pack-Monkhorst (PM) shrinking factor was used to specify the number of equidistant k-points taken along each direction of b1, b2 and b3 in one reciprocal lattice PUC. (https://books.google.co.uk/books?id=nX_wG7WaDJsC&pg=PA38&lpg=PA38&dq=Pack-Monkhorst+shrinking+factor&source=bl&ots=vL_-nToT5e&sig=SOp4EsY7oG-ki9tlvsMSRTJ-eTY&hl=en&sa=X&ved=0ahUKEwiG-PHajZfZAhXmLsAKHefcBDAQ6AEINzAC#v=onepage&q=Pack-Monkhorst%20shrinking%20factor&f=false) The Cartesian coordinates of the k-points calculated are given by the equation

[EQUATION FROM https://journals.aps.org/prb/pdf/10.1103/PhysRevB.93.155109]

A major advantage is its computational efficiency by restricting the number of k-points calculated to a finite value. Moreover, the accuracy obtained from calculations with a PUC can be comparable to that of a supercell as long as the shrinking factor is appropriate.

Table 1 illustrates the effect of modifying the PM shrinking factor on the number of k-points calculated.

'''Table 1'''. Grid size against number of k-points
{| class="wikitable" style="text-align: center; width: 85%;margin: auto;"
! Grid Size (n x n x n)
! Number of k-points
|-
| 1
| 1
|-
| 2
| 4
|-
| 3
| 18
|-
| 4
| 32
|-
| 5
| 75
|-
| 6
| 108
|-
| 8
| 256
|-
| 10
| 500
|-
| 16
| 2048
|-
| 20
| 4000
|-
| 64
| >99 999
|}

As the mesh of k-points increases, the number of k-points calculated increases as well. This is contrary to the prediction from the above equation, where we would expect kx * ky* kz number of points. This can be attributed to the mapping of equivalent k-points onto each other and thus the number of k-points calculated is reduced.

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased.

[[File:Syl815FEshinking.png|thumb|center|600px|Total energy vs. grid size]]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication, and will be further discussed in section 3.

An initial plot of the density of states was obtained from a 1x1x1 grid yielding six resultant modes. Sharp and distinct peaks are observed in the plot, since only one k-point was sampled.

Notably, only four unique peaks are observed even though we should observe 6 modes of vibrations. The final two modes are degenerate at _____ and _____ wavenumbers. Compared to the non-degenerate acoustic and optical peaks (___ and ____ respectively), the degenerate acoustic modes are higher in energy whereas the degenerate optical modes are lower in energy correspondingly. It can therefore be deduced that the degenerate acoustic and optical modes are transverse in nature.

The k-point used in the DOS calculation could be identified by comparing with the dispersion curve. Since point M contains all of the frequency values in Figure ___, it can be determined that the point represented in the DOS curve is M, where kx = 0.5, ky = 0.5 and kz = 0.5.

=== Relationship between the Dispersion Curve and DOS ===
The DOS curve illustrates the number of energy states per unit energy, demonstrating a mode at 414 cm-1. This correlates well with Figure ____. By constructing a horizontal line at frequency = 414 cm-1, it can be observed that the branches intersect this line frequently. This implies that a significant proportion of k-points have vibrational modes of frequency 414 cm-1.

Loss of information in dispersion curve: how is that so?
Orthogonal to the dispersion curve: gives the number of energy states.
Information provided for the dispersion curve vs the density of states.
=== Computing the Free Energy Using the Harmonic Approximation ===
The figure below demonstrates the relationship between the PM shrinking factor used and the computed Helmholtz free energy of the system.

[[file:Syl815FEvsT.png|thumb|center|600px|Free Energy vs. PM shrinking factor.]]

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[[file:Syl815FEshinking.png|thumb|center|600px|DOS against wavenumber for various shrinking factors.]]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication.

From Figure ____, the free energy of MgO is observed to increase and converge to a value of -40.926 483 eV, and it is observed that this occurs for a grid size of 8x8x8. (Why does it do so? Explain here)

A 2x2x2 grid is sufficient for calculating the free energy of MgO to 1 meV. A 4x4x4 grid is necessary for a precision to 0.5 meV and 0.1 meV.

=== Thermal Expansion ===
The Helmholtz free energy of a crystal is given by the sum of the energies of independent vibrational waves. The energy level <math display="inline">n</math> of a quantum harmonic oscillator are given by:

<math display="block">E_n = \left( n+ \frac{1}{2} \right) h \nu </math>

where <math display="inline">h</math> is Planck's constant and <math display="inline">\nu</math> is the frequency of energy level <math display="inline">n</math>. For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators, the vibrational energy is given by:

<math display="block">E_{vib} = \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu</math>

For a canonical <math display="inline">NVT</math> ensemble, the partition function is given by

<math display="block">Z = \sum_n e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

where <math display="inline">\beta = \frac{1}{kT}</math> and <math display="inline">E_n</math> enumerates all vibrational energy states.

For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators,

<math display="block">Z_N = \prod_n^{3N} e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

The phonon entropy can then be expressed in terms of the partition function:
<math display="block">S=-k_B \ln Z_N</math> where <math display="inline">k_B</math> is the Boltzmann constant.

Given the relation
<math display="block">F=U+TS</math>
where <math display="inline">U</math> is the internal energy of the system— for a crystal this is its electric potential energy
<math display="block">U_E = \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}}</math>

where <math display="inline">i</math> and <math display="inline">j</math> are the indices of the ions, <math display="inline">r_{ij}</math> is the distance between <math display="inline">i</math> and <math display="inline">j</math> and <math display="inline">\varepsilon_0 = 8.8542 \times 10^{-12} F\cdot m^{-1} </math>

The Helmholtz free energy of a crystal is thus given by

<math display="block">F= \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}} + \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu + k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math>

This equation could be used to qualitatively rationalise the free energy dependence on temperature. The data obtained is plotted in Figure ____.
[FIGURE]
Particularly, there are two salient regimes of interest. At low temperatures, T < 100 K, the graph is flat. However, at high temperatures, the behaviour is approximately linear. These observations are in agreement with the above equation, which highlights the temperature dependence of entropy <math display="inline">S</math>. At low temperatures, the term <math> k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> is extremely small, and hence the free energy term is dominated by the internal energy of the crystal. At high temperatures, the term <math>-k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> dominates and therefore the free energy of the system appears to have a dependence in temperature.

==== Variation of Lattice Parameter with Temperature ====
[[File:Syl815CellPvsT.png|thumb|center|600px|This figure illustrates the variation in cell parameter of MgO with temperature]]

As the temperature increases, the lattice parameter increases. It can thus be observed that the cell volume <math display="inline">V </math> has a dependence on temperature <math display="inline">T </math>, and the thermal expansion coefficient <math display="inline">\alpha </math> is given by

<math display="block">\alpha = \frac{1}{3V} \left( \frac{\partial V}{\partial T}\right)_P = \frac{1}{3B} \left( \frac{\partial P}{\partial T}\right)_V </math>
Where <math display="inline">B </math> is the bulk modulus and <math display="inline">P </math> is the pressure.
At 300 K, <math display="inline">\alpha = 2.06\times 10^{-5} K^{-1}</math>, compared to a literature value of <math display="inline">\alpha = 3.06\times 10^{-5} K^{-1}</math>

=== MD Simulation ===
The cell volume per formula unit of MgO was plotted against temperatures between 20 K to 4000 K.
[[File:Syl815MD.png|thumb|center|600px|"Experimental Data MD"]]
Under MD, the cell volume generally increases linearly with temperature throughout. By considering the mean kinetic energy of the crystal
<math display="block">\left \langle E_k \right \rangle = \frac{1}{2} M \sum_{i=1}^N v_i^2 = \frac{3}{2} Nk_BT_{MD} </math>
Where <math display="inline">\left \langle E_k \right \rangle </math> is the average kinetic energy of the atoms, <math display="inline">M</math> is the mass of the crystal lattice, and <math display="inline">v_i</math> represents the velocity of the atom <math display="inline">i</math>. It can be observed that the cell energy is linearly dependent on temperature. In a constant pressure system, this would result in volume expansion as temperature increases.

It can be observed that at high temperatures when <math display="inline">T\ge 2000 K </math>, more noise is present in the data due to the large cell volume and the large kinetic energy of the atoms.

[[File:Syl_MDvsQHA.png|thumb|center|600px|"This figure compares the data obtained for the thermal expansion of MgO under QHA and under MD."]]

At extremely low temperatures of <math display="inline">T\le 200 K </math>, QHA predicts a larger cell volume than MD. This can be attributed to the significant quantum effects at such low temperatures. Since MD only accounts for the kinetic energy of the atoms and neglects zero point vibrations, it predicts a smaller cell volume with the atoms closer together.

The data obtained for MD and QHA demonstrate strong agreement for temperatures between 200 to 1000 K. At these temperatures, the thermal energy of the system is sufficiently large such that the motion of the particles can be described classically.

Links:
https://pubs.acs.org/doi/pdf/10.1021/ct0500904
https://journals.aps.org/prb/pdf/10.1103/PhysRevB.43.5024

Conclusion

Between MD and Thermal expansion model
Compare quantitatively with values
Why you think it shouldn’t be linear

== References ==

File:Syl815Phonon Dispersion Graph 100.PNG

2018-02-14T11:58:51Z

Syl815:

Rep:MgO:syl815

2018-02-14T11:58:07Z

Syl815: /* Determining Optimal Grid Size for MgO */

== Abstract ==
== Introduction ==
Thermal properties of materials are described by their vibrational free energies, which can be described in terms of the relative motion of atoms or the motion of their centre-of-mass. <ref name='phonons'>G. Srivastava, The physics of phonons, A. Hilger, Bristol, 1990.</ref> These concepts give rise to different approaches in calculating vibrational free energy, and both methodologies will be explored in greater detail.

=== Quasi-harmonic Approximation (QHA) ===
Fundamentally, QHA invokes the description of a crystalline solid as a primitive unit cell. This is essential due to the impracticality of calculating all the vibrational degrees of freedom in a crystal—for a crystal of size <math>N</math>, there are <math>3N</math> degrees of vibrational freedom, and in an infinitely large crystal lattice, 3N --> infinity.
Nonetheless, the translational periodicity of the crystal lattice, where <math display="inline">f(r + T) = f(r)</math>, simplifies the dynamics of all atoms in the lattice into that of a unit cell. For such a simplification to be appropriate, the following assumptions are made.

==== The Adiabatic Approximation ====
The adiabatic approximation separates the motion of the ion cores from that of the electrons since former are much more massive than the latter. Hence, the ion cores can be assumed to be in their equilibrium positions and that their motion is dependent on the potential field generated from the average motion of electrons.<ref name='phonons' />

==== The Harmonic Approximation ====
The total potential energy of a crystal can be expressed as the sum of all interatomic potentials. A two-body system typically has an anharmonic potential energy surface <math display="inline">U\left( r\right) </math>, where <math display="inline">r</math> is the interatomic separation. By considering a small displacement <math display="inline">x</math>,
<math display="block">x = r-r_0 \left( 1 \right)</math>

Where <math display="inline">r_0</math> is the equilibrium distance between the first and second atoms and is a minimum on the potential energy surface.

<math display="inline">U\left( r\right) </math> can be expanded in a Taylor series about <math display="inline">r_0</math>.
<math display="block">U\left( r\right) = U(r_0) + \frac{\partial U}{\partial x} x + \frac{\partial^2 U}{\partial x^2} x^2 + ... </math>
Since <math display="inline">U(r_0) </math> is unimportant in dynamics, <math display = "inline">\frac{\partial U}{\partial x}</math> is a force term and must be 0 for an equilibrium configuration, and all higher order terms <math display = "inline">x_n</math> , where <math display = "inline">n \ge 3</math> are assumed to be close to 0. As such, only the quadratic term is considered in the harmonic approximation. The solutions are the normal modes of vibrations for a system of independent quantum oscillators.

A phonon is a quantum of vibrational energy, hw, associated with a wave vector k.

Hence, for a crystal, its potential energy is given in the following equation.
[First principles, pdf, page 2]
Where l and k are the labels of the unit cells and atoms in each unit cell respectively. (http://dx.doi.org/10.1016/j.scriptamat.2015.07.021)

===== Limitations of Harmonic Approximation =====
The harmonic approximation predicts symmetric atomic vibrations about r0 at all temperatures, and is therefore incongruent with observed phenomena such as thermal expansion and heat conductivity. (https://core.ac.uk/download/pdf/42335965.pdf) The QHA causes renormalisation of the phonon frequencies and atomic force constants as is appropriate for the thermal equation of state. (https://ac.els-cdn.com/S0081194708606566/1-s2.0-S0081194708606566-main.pdf?_tid=f9617aac-104c-11e8-87bc-00000aab0f27&acdnat=1518478522_8ee2607e5ac970b2f00844b9ef59e8dc)

=== MD Simulation ===

=== System considered ===
MgO has a face-centred cubic lattice; more specifically, an alkali halide type structure with a Fm3m space group. (Citation)
MgO: insulator
Lattice parameters, symmetry, space group, etc.

=== Aims of the Exercise ===

== Methodology ==
Unless otherwise stated, all calculations were performed on a primitive unit cell of MgO with lattice parameters a = 2.9783 Angstrom, a = b = c; alpha = 60 degrees, and alpha = beta = gamma with GULP version 1.4.43 and crystals visualised with DLV interface.

A phonon dispersion curve was computed by sampling 100 points within the first Brillouin zone. The phonon density of states (DOS) was calculated with various shrinking factors, and the graphs subsequently plotted with matplotlib. The free energy of MgO was calculated with different shrinking factors at 300 K, and a suitable shrinking factor selected for the subsequent investigation of the thermal expansion of MgO. For every run, the Gibbs free energy was optimised, and calculations were performed from 0 to 2960 K in temperature steps of 20 K.

All MD simulations were performed on an isothermal-isobaric ensemble of MgO supercell of 32 formula units, with the following cell parameters:
A = 8.4239 Angstrom
A = b = c
Alpha = 90 degrees
Alpha = beta = gamma

MD was performed over a temperature range of 20 K to 4000 K, with temperature steps of 20 K. All calculations were performed with a time step of 1 fs. From 20 K to 1680 K, the system was allowed to first equilibrate for 1 ps; this was increased to 5 ps from 1700 K to 4000 K. Following which, MD production was allowed to run for 5 ps for all temperatures.

All data was analysed with Python on Jupyter notebook, and all graphs plotted with matplotlib.

== Results and Discussion ==
The lattice energy of MgO was calculated to be -41.0753 eV per primitive unit cell.
=== Phonon Modes of MgO ===
Figure 1 illustrates the phonon dispersion curve computed at 100 points for the primitive unit cell.

[FIGURE HERE]

Computation of the

A salient feature is the presence of 6 branches in the dispersion diagram. Assuming that the Born-von Karman boundary condition is satisfied, the edge effects of cells on dynamics can be ignored and u(N+1) = u1, where u is the displacement and N is the number of unit cells. This also implies the translational symmetry in k-space, such that all information of phonon dispersion can be derived by sampling in the first Brillouin zone.

By considering a linear diatomic chain satisfying the periodic boundary condition, the solutions to the vibrational frequency can be expressed in the form
[EQUATION, pg. 27 of book]

Highlighting two possible solutions for each k-value in a linear chain. Moreover, when m1 /= m2, a gap is observed at k = pi/2a, which is observed in Figure 1. (https://hal.archives-ouvertes.fr/file/index/docid/247340/filename/ajp-jp1v7p509.pdf)

By extending the logic to a 3D crystal lattice, the number of branches observed is given by 3x, where x is the number of atoms per unit cell. This is in agreement with the observation in Figure 1.

By appraising the solutions for k = 0 (long wavelength limit),

w1 = 2lambda(1/m + 1/M)
w2 = 0

w1 corresponds to a high energy mode where the atoms in the unit cell are moving out-of-phase, where frequency values are within the visible electromagnetic spectrum. The atoms are able to interact with an electric field of appropriate frequency due to the presence of both a positive and negative charge within the unit cell. It is hence naturally termed the optical mode. (Introduction to lattice dynamics)

On the other hand, w2 corresponds to a low energy mode with the atoms moving in phase and the wave pattern is similar to sound waves—hence the term acoustic mode. For any crystal with N atoms in the unit cell, there are only 3 acoustic—2 transverse and 1 longitudinal—and 3N-3 optical branches. The transverse modes are perpendicular to k, while the longitudinal mode is parallel.

For a cubic crystal, the highly symmetric nature indicates the possibility for some vibrations to be degenerate. Typically, transverse modes are lower in energy due to the weaker interaction between atoms in the unit cell.

This is encapsulated in Bloch's theorem and thus the wavefunction psi can be expressed in the form as given in equation 1.

[EQUATION HERE]

While typically applied to electrons in crystals, Bloch's theorem is broadly applicable in describing periodic wave phenomena, such as in phononic crystals. The branches relate to the modes of vibration, as given by equation for a one-dimensional system.

[EQUATION HERE]

Since the system is a three-dimensional one, vibrations can occur independently in perpendicular planes, therefore giving rise to additional normal coordinates and vibrations.

A dispersion diagram enables clear visualisation of the nature of the band gap, as illustrated in Figure ___. Figure __ clearly illustrates acoustic and optical phonons, depicting the in-phase and out-of-phase movement of the ions respectively.

Nonetheless, it is limited since energy values are unknown for points which are not sampled. Calculation of the density of states has been proven to be more useful.

=== Computing Density of States (DOS) ===
The impracticality of sampling all k-points within the FBZ can be circumvented by the use of a commensurate grid of k-points. To determine this set of k-points, the Pack-Monkhorst (PM) shrinking factor was used to specify the number of equidistant k-points taken along each direction of b1, b2 and b3 in one reciprocal lattice PUC. (https://books.google.co.uk/books?id=nX_wG7WaDJsC&pg=PA38&lpg=PA38&dq=Pack-Monkhorst+shrinking+factor&source=bl&ots=vL_-nToT5e&sig=SOp4EsY7oG-ki9tlvsMSRTJ-eTY&hl=en&sa=X&ved=0ahUKEwiG-PHajZfZAhXmLsAKHefcBDAQ6AEINzAC#v=onepage&q=Pack-Monkhorst%20shrinking%20factor&f=false) The Cartesian coordinates of the k-points calculated are given by the equation

[EQUATION FROM https://journals.aps.org/prb/pdf/10.1103/PhysRevB.93.155109]

A major advantage is its computational efficiency by restricting the number of k-points calculated to a finite value. Moreover, the accuracy obtained from calculations with a PUC can be comparable to that of a supercell as long as the shrinking factor is appropriate.

Table 1 illustrates the effect of modifying the PM shrinking factor on the number of k-points calculated.

'''Table 1'''. Grid size against number of k-points
{| class="wikitable" style="text-align: center; width: 85%;margin: auto;"
! Grid Size (n x n x n)
! Number of k-points
|-
| 1
| 1
|-
| 2
| 4
|-
| 3
| 18
|-
| 4
| 32
|-
| 5
| 75
|-
| 6
| 108
|-
| 8
| 256
|-
| 10
| 500
|-
| 16
| 2048
|-
| 20
| 4000
|-
| 64
| >99 999
|}

As the mesh of k-points increases, the number of k-points calculated increases as well. This is contrary to the prediction from the above equation, where we would expect kx * ky* kz number of points. This can be attributed to the mapping of equivalent k-points onto each other and thus the number of k-points calculated is reduced.

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased.

[[File:Syl815FEshinking.png|thumb|center|600px|Total energy vs. grid size]]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication, and will be further discussed in section 3.

An initial plot of the density of states was obtained from a 1x1x1 grid yielding six resultant modes. Sharp and distinct peaks are observed in the plot, since only one k-point was sampled.

Notably, only four unique peaks are observed even though we should observe 6 modes of vibrations. The final two modes are degenerate at _____ and _____ wavenumbers. Compared to the non-degenerate acoustic and optical peaks (___ and ____ respectively), the degenerate acoustic modes are higher in energy whereas the degenerate optical modes are lower in energy correspondingly. It can therefore be deduced that the degenerate acoustic and optical modes are transverse in nature.

The k-point used in the DOS calculation could be identified by comparing with the dispersion curve. Since point M contains all of the frequency values in Figure ___, it can be determined that the point represented in the DOS curve is M, where kx = 0.5, ky = 0.5 and kz = 0.5.

=== Relationship between the Dispersion Curve and DOS ===
The DOS curve illustrates the number of energy states per unit energy, demonstrating a mode at 414 cm-1. This correlates well with Figure ____. By constructing a horizontal line at frequency = 414 cm-1, it can be observed that the branches intersect this line frequently. This implies that a significant proportion of k-points have vibrational modes of frequency 414 cm-1.

Loss of information in dispersion curve: how is that so?
Orthogonal to the dispersion curve: gives the number of energy states.
Information provided for the dispersion curve vs the density of states.
=== Computing the Free Energy Using the Harmonic Approximation ===
The figure below demonstrates the relationship between the PM shrinking factor used and the computed Helmholtz free energy of the system.

[[file:Syl815FEvsT.png|thumb|center|600px|Free Energy vs. PM shrinking factor.]]

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[[file:Syl815FEshinking.png|thumb|center|600px|DOS against wavenumber for various shrinking factors.]]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication.

From Figure ____, the free energy of MgO is observed to increase and converge to a value of -40.926 483 eV, and it is observed that this occurs for a grid size of 8x8x8. (Why does it do so? Explain here)

A 2x2x2 grid is sufficient for calculating the free energy of MgO to 1 meV. A 4x4x4 grid is necessary for a precision to 0.5 meV and 0.1 meV.

=== Thermal Expansion ===
The Helmholtz free energy of a crystal is given by the sum of the energies of independent vibrational waves. The energy level <math display="inline">n</math> of a quantum harmonic oscillator are given by:

<math display="block">E_n = \left( n+ \frac{1}{2} \right) h \nu </math>

where <math display="inline">h</math> is Planck's constant and <math display="inline">\nu</math> is the frequency of energy level <math display="inline">n</math>. For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators, the vibrational energy is given by:

<math display="block">E_{vib} = \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu</math>

For a canonical <math display="inline">NVT</math> ensemble, the partition function is given by

<math display="block">Z = \sum_n e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

where <math display="inline">\beta = \frac{1}{kT}</math> and <math display="inline">E_n</math> enumerates all vibrational energy states.

For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators,

<math display="block">Z_N = \prod_n^{3N} e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

The phonon entropy can then be expressed in terms of the partition function:
<math display="block">S=-k_B \ln Z_N</math> where <math display="inline">k_B</math> is the Boltzmann constant.

Given the relation
<math display="block">F=U+TS</math>
where <math display="inline">U</math> is the internal energy of the system— for a crystal this is its electric potential energy
<math display="block">U_E = \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}}</math>

where <math display="inline">i</math> and <math display="inline">j</math> are the indices of the ions, <math display="inline">r_{ij}</math> is the distance between <math display="inline">i</math> and <math display="inline">j</math> and <math display="inline">\varepsilon_0 = 8.8542 \times 10^{-12} F\cdot m^{-1} </math>

The Helmholtz free energy of a crystal is thus given by

<math display="block">F= \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}} + \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu + k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math>

This equation could be used to qualitatively rationalise the free energy dependence on temperature. The data obtained is plotted in Figure ____.
[FIGURE]
Particularly, there are two salient regimes of interest. At low temperatures, T < 100 K, the graph is flat. However, at high temperatures, the behaviour is approximately linear. These observations are in agreement with the above equation, which highlights the temperature dependence of entropy <math display="inline">S</math>. At low temperatures, the term <math> k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> is extremely small, and hence the free energy term is dominated by the internal energy of the crystal. At high temperatures, the term <math>-k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> dominates and therefore the free energy of the system appears to have a dependence in temperature.

==== Variation of Lattice Parameter with Temperature ====
[[File:Syl815CellPvsT.png|thumb|center|600px|This figure illustrates the variation in cell parameter of MgO with temperature]]

As the temperature increases, the lattice parameter increases. It can thus be observed that the cell volume <math display="inline">V </math> has a dependence on temperature <math display="inline">T </math>, and the thermal expansion coefficient <math display="inline">\alpha </math> is given by

<math display="block">\alpha = \frac{1}{3V} \left( \frac{\partial V}{\partial T}\right)_P = \frac{1}{3B} \left( \frac{\partial P}{\partial T}\right)_V </math>
Where <math display="inline">B </math> is the bulk modulus and <math display="inline">P </math> is the pressure.
At 300 K, <math display="inline">\alpha = 2.06\times 10^{-5} K^{-1}</math>, compared to a literature value of <math display="inline">\alpha = 3.06\times 10^{-5} K^{-1}</math>

=== MD Simulation ===
The cell volume per formula unit of MgO was plotted against temperatures between 20 K to 4000 K.
[[File:Syl815MD.png|thumb|center|600px|"Experimental Data MD"]]
Under MD, the cell volume generally increases linearly with temperature throughout. By considering the mean kinetic energy of the crystal
<math display="block">\left \langle E_k \right \rangle = \frac{1}{2} M \sum_{i=1}^N v_i^2 = \frac{3}{2} Nk_BT_{MD} </math>
Where <math display="inline">\left \langle E_k \right \rangle </math> is the average kinetic energy of the atoms, <math display="inline">M</math> is the mass of the crystal lattice, and <math display="inline">v_i</math> represents the velocity of the atom <math display="inline">i</math>. It can be observed that the cell energy is linearly dependent on temperature. In a constant pressure system, this would result in volume expansion as temperature increases.

It can be observed that at high temperatures when <math display="inline">T\ge 2000 K </math>, more noise is present in the data due to the large cell volume and the large kinetic energy of the atoms.

[[File:Syl_MDvsQHA.png|thumb|center|600px|"This figure compares the data obtained for the thermal expansion of MgO under QHA and under MD."]]

At extremely low temperatures of <math display="inline">T\le 200 K </math>, QHA predicts a larger cell volume than MD. This can be attributed to the significant quantum effects at such low temperatures. Since MD only accounts for the kinetic energy of the atoms and neglects zero point vibrations, it predicts a smaller cell volume with the atoms closer together.

The data obtained for MD and QHA demonstrate strong agreement for temperatures between 200 to 1000 K. At these temperatures, the thermal energy of the system is sufficiently large such that the motion of the particles can be described classically.

Links:
https://pubs.acs.org/doi/pdf/10.1021/ct0500904
https://journals.aps.org/prb/pdf/10.1103/PhysRevB.43.5024

Conclusion

Between MD and Thermal expansion model
Compare quantitatively with values
Why you think it shouldn’t be linear

== References ==

Rep:MgO:syl815

2018-02-14T11:56:32Z

Syl815: /* Computing Density of States (DOS) */

== Abstract ==
== Introduction ==
Thermal properties of materials are described by their vibrational free energies, which can be described in terms of the relative motion of atoms or the motion of their centre-of-mass. <ref name='phonons'>G. Srivastava, The physics of phonons, A. Hilger, Bristol, 1990.</ref> These concepts give rise to different approaches in calculating vibrational free energy, and both methodologies will be explored in greater detail.

=== Quasi-harmonic Approximation (QHA) ===
Fundamentally, QHA invokes the description of a crystalline solid as a primitive unit cell. This is essential due to the impracticality of calculating all the vibrational degrees of freedom in a crystal—for a crystal of size <math>N</math>, there are <math>3N</math> degrees of vibrational freedom, and in an infinitely large crystal lattice, 3N --> infinity.
Nonetheless, the translational periodicity of the crystal lattice, where <math display="inline">f(r + T) = f(r)</math>, simplifies the dynamics of all atoms in the lattice into that of a unit cell. For such a simplification to be appropriate, the following assumptions are made.

==== The Adiabatic Approximation ====
The adiabatic approximation separates the motion of the ion cores from that of the electrons since former are much more massive than the latter. Hence, the ion cores can be assumed to be in their equilibrium positions and that their motion is dependent on the potential field generated from the average motion of electrons.<ref name='phonons' />

==== The Harmonic Approximation ====
The total potential energy of a crystal can be expressed as the sum of all interatomic potentials. A two-body system typically has an anharmonic potential energy surface <math display="inline">U\left( r\right) </math>, where <math display="inline">r</math> is the interatomic separation. By considering a small displacement <math display="inline">x</math>,
<math display="block">x = r-r_0 \left( 1 \right)</math>

Where <math display="inline">r_0</math> is the equilibrium distance between the first and second atoms and is a minimum on the potential energy surface.

<math display="inline">U\left( r\right) </math> can be expanded in a Taylor series about <math display="inline">r_0</math>.
<math display="block">U\left( r\right) = U(r_0) + \frac{\partial U}{\partial x} x + \frac{\partial^2 U}{\partial x^2} x^2 + ... </math>
Since <math display="inline">U(r_0) </math> is unimportant in dynamics, <math display = "inline">\frac{\partial U}{\partial x}</math> is a force term and must be 0 for an equilibrium configuration, and all higher order terms <math display = "inline">x_n</math> , where <math display = "inline">n \ge 3</math> are assumed to be close to 0. As such, only the quadratic term is considered in the harmonic approximation. The solutions are the normal modes of vibrations for a system of independent quantum oscillators.

A phonon is a quantum of vibrational energy, hw, associated with a wave vector k.

Hence, for a crystal, its potential energy is given in the following equation.
[First principles, pdf, page 2]
Where l and k are the labels of the unit cells and atoms in each unit cell respectively. (http://dx.doi.org/10.1016/j.scriptamat.2015.07.021)

===== Limitations of Harmonic Approximation =====
The harmonic approximation predicts symmetric atomic vibrations about r0 at all temperatures, and is therefore incongruent with observed phenomena such as thermal expansion and heat conductivity. (https://core.ac.uk/download/pdf/42335965.pdf) The QHA causes renormalisation of the phonon frequencies and atomic force constants as is appropriate for the thermal equation of state. (https://ac.els-cdn.com/S0081194708606566/1-s2.0-S0081194708606566-main.pdf?_tid=f9617aac-104c-11e8-87bc-00000aab0f27&acdnat=1518478522_8ee2607e5ac970b2f00844b9ef59e8dc)

=== MD Simulation ===

=== System considered ===
MgO has a face-centred cubic lattice; more specifically, an alkali halide type structure with a Fm3m space group. (Citation)
MgO: insulator
Lattice parameters, symmetry, space group, etc.

=== Aims of the Exercise ===

== Methodology ==
Unless otherwise stated, all calculations were performed on a primitive unit cell of MgO with lattice parameters a = 2.9783 Angstrom, a = b = c; alpha = 60 degrees, and alpha = beta = gamma with GULP version 1.4.43 and crystals visualised with DLV interface.

A phonon dispersion curve was computed by sampling 100 points within the first Brillouin zone. The phonon density of states (DOS) was calculated with various shrinking factors, and the graphs subsequently plotted with matplotlib. The free energy of MgO was calculated with different shrinking factors at 300 K, and a suitable shrinking factor selected for the subsequent investigation of the thermal expansion of MgO. For every run, the Gibbs free energy was optimised, and calculations were performed from 0 to 2960 K in temperature steps of 20 K.

All MD simulations were performed on an isothermal-isobaric ensemble of MgO supercell of 32 formula units, with the following cell parameters:
A = 8.4239 Angstrom
A = b = c
Alpha = 90 degrees
Alpha = beta = gamma

MD was performed over a temperature range of 20 K to 4000 K, with temperature steps of 20 K. All calculations were performed with a time step of 1 fs. From 20 K to 1680 K, the system was allowed to first equilibrate for 1 ps; this was increased to 5 ps from 1700 K to 4000 K. Following which, MD production was allowed to run for 5 ps for all temperatures.

All data was analysed with Python on Jupyter notebook, and all graphs plotted with matplotlib.

== Results and Discussion ==
The lattice energy of MgO was calculated to be -41.0753 eV per primitive unit cell.
=== Phonon Modes of MgO ===
Figure 1 illustrates the phonon dispersion curve computed at 100 points for the primitive unit cell.

[FIGURE HERE]

Computation of the

A salient feature is the presence of 6 branches in the dispersion diagram. Assuming that the Born-von Karman boundary condition is satisfied, the edge effects of cells on dynamics can be ignored and u(N+1) = u1, where u is the displacement and N is the number of unit cells. This also implies the translational symmetry in k-space, such that all information of phonon dispersion can be derived by sampling in the first Brillouin zone.

By considering a linear diatomic chain satisfying the periodic boundary condition, the solutions to the vibrational frequency can be expressed in the form
[EQUATION, pg. 27 of book]

Highlighting two possible solutions for each k-value in a linear chain. Moreover, when m1 /= m2, a gap is observed at k = pi/2a, which is observed in Figure 1. (https://hal.archives-ouvertes.fr/file/index/docid/247340/filename/ajp-jp1v7p509.pdf)

By extending the logic to a 3D crystal lattice, the number of branches observed is given by 3x, where x is the number of atoms per unit cell. This is in agreement with the observation in Figure 1.

By appraising the solutions for k = 0 (long wavelength limit),

w1 = 2lambda(1/m + 1/M)
w2 = 0

w1 corresponds to a high energy mode where the atoms in the unit cell are moving out-of-phase, where frequency values are within the visible electromagnetic spectrum. The atoms are able to interact with an electric field of appropriate frequency due to the presence of both a positive and negative charge within the unit cell. It is hence naturally termed the optical mode. (Introduction to lattice dynamics)

On the other hand, w2 corresponds to a low energy mode with the atoms moving in phase and the wave pattern is similar to sound waves—hence the term acoustic mode. For any crystal with N atoms in the unit cell, there are only 3 acoustic—2 transverse and 1 longitudinal—and 3N-3 optical branches. The transverse modes are perpendicular to k, while the longitudinal mode is parallel.

For a cubic crystal, the highly symmetric nature indicates the possibility for some vibrations to be degenerate. Typically, transverse modes are lower in energy due to the weaker interaction between atoms in the unit cell.

This is encapsulated in Bloch's theorem and thus the wavefunction psi can be expressed in the form as given in equation 1.

[EQUATION HERE]

While typically applied to electrons in crystals, Bloch's theorem is broadly applicable in describing periodic wave phenomena, such as in phononic crystals. The branches relate to the modes of vibration, as given by equation for a one-dimensional system.

[EQUATION HERE]

Since the system is a three-dimensional one, vibrations can occur independently in perpendicular planes, therefore giving rise to additional normal coordinates and vibrations.

A dispersion diagram enables clear visualisation of the nature of the band gap, as illustrated in Figure ___. Figure __ clearly illustrates acoustic and optical phonons, depicting the in-phase and out-of-phase movement of the ions respectively.

Nonetheless, it is limited since energy values are unknown for points which are not sampled. Calculation of the density of states has been proven to be more useful.

=== Computing Density of States (DOS) ===
The impracticality of sampling all k-points within the FBZ can be circumvented by the use of a commensurate grid of k-points. To determine this set of k-points, the Pack-Monkhorst (PM) shrinking factor was used to specify the number of equidistant k-points taken along each direction of b1, b2 and b3 in one reciprocal lattice PUC. (https://books.google.co.uk/books?id=nX_wG7WaDJsC&pg=PA38&lpg=PA38&dq=Pack-Monkhorst+shrinking+factor&source=bl&ots=vL_-nToT5e&sig=SOp4EsY7oG-ki9tlvsMSRTJ-eTY&hl=en&sa=X&ved=0ahUKEwiG-PHajZfZAhXmLsAKHefcBDAQ6AEINzAC#v=onepage&q=Pack-Monkhorst%20shrinking%20factor&f=false) The Cartesian coordinates of the k-points calculated are given by the equation

[EQUATION FROM https://journals.aps.org/prb/pdf/10.1103/PhysRevB.93.155109]

A major advantage is its computational efficiency by restricting the number of k-points calculated to a finite value. Moreover, the accuracy obtained from calculations with a PUC can be comparable to that of a supercell as long as the shrinking factor is appropriate.

Table 1 illustrates the effect of modifying the PM shrinking factor on the number of k-points calculated.

'''Table 1'''. Grid size against number of k-points
{| class="wikitable" style="text-align: center; width: 85%;margin: auto;"
! Grid Size (n x n x n)
! Number of k-points
|-
| 1
| 1
|-
| 2
| 4
|-
| 3
| 18
|-
| 4
| 32
|-
| 5
| 75
|-
| 6
| 108
|-
| 8
| 256
|-
| 10
| 500
|-
| 16
| 2048
|-
| 20
| 4000
|-
| 64
| >99 999
|}

As the mesh of k-points increases, the number of k-points calculated increases as well. This is contrary to the prediction from the above equation, where we would expect kx * ky* kz number of points. This can be attributed to the mapping of equivalent k-points onto each other and thus the number of k-points calculated is reduced.

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[GRAPH OF ENERGY AGAINST K-POINTS]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication, and will be further discussed in section 3.

An initial plot of the density of states was obtained from a 1x1x1 grid yielding six resultant modes. Sharp and distinct peaks are observed in the plot, since only one k-point was sampled.

Notably, only four unique peaks are observed even though we should observe 6 modes of vibrations. The final two modes are degenerate at _____ and _____ wavenumbers. Compared to the non-degenerate acoustic and optical peaks (___ and ____ respectively), the degenerate acoustic modes are higher in energy whereas the degenerate optical modes are lower in energy correspondingly. It can therefore be deduced that the degenerate acoustic and optical modes are transverse in nature.

The k-point used in the DOS calculation could be identified by comparing with the dispersion curve. Since point M contains all of the frequency values in Figure ___, it can be determined that the point represented in the DOS curve is M, where kx = 0.5, ky = 0.5 and kz = 0.5.

=== Relationship between the Dispersion Curve and DOS ===
The DOS curve illustrates the number of energy states per unit energy, demonstrating a mode at 414 cm-1. This correlates well with Figure ____. By constructing a horizontal line at frequency = 414 cm-1, it can be observed that the branches intersect this line frequently. This implies that a significant proportion of k-points have vibrational modes of frequency 414 cm-1.

Loss of information in dispersion curve: how is that so?
Orthogonal to the dispersion curve: gives the number of energy states.
Information provided for the dispersion curve vs the density of states.
=== Computing the Free Energy Using the Harmonic Approximation ===
The figure below demonstrates the relationship between the PM shrinking factor used and the computed Helmholtz free energy of the system.

[[file:Syl815FEvsT.png|thumb|center|600px|Free Energy vs. PM shrinking factor.]]

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[[file:Syl815FEshinking.png|thumb|center|600px|DOS against wavenumber for various shrinking factors.]]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication.

From Figure ____, the free energy of MgO is observed to increase and converge to a value of -40.926 483 eV, and it is observed that this occurs for a grid size of 8x8x8. (Why does it do so? Explain here)

A 2x2x2 grid is sufficient for calculating the free energy of MgO to 1 meV. A 4x4x4 grid is necessary for a precision to 0.5 meV and 0.1 meV.

=== Thermal Expansion ===
The Helmholtz free energy of a crystal is given by the sum of the energies of independent vibrational waves. The energy level <math display="inline">n</math> of a quantum harmonic oscillator are given by:

<math display="block">E_n = \left( n+ \frac{1}{2} \right) h \nu </math>

where <math display="inline">h</math> is Planck's constant and <math display="inline">\nu</math> is the frequency of energy level <math display="inline">n</math>. For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators, the vibrational energy is given by:

<math display="block">E_{vib} = \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu</math>

For a canonical <math display="inline">NVT</math> ensemble, the partition function is given by

<math display="block">Z = \sum_n e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

where <math display="inline">\beta = \frac{1}{kT}</math> and <math display="inline">E_n</math> enumerates all vibrational energy states.

For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators,

<math display="block">Z_N = \prod_n^{3N} e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

The phonon entropy can then be expressed in terms of the partition function:
<math display="block">S=-k_B \ln Z_N</math> where <math display="inline">k_B</math> is the Boltzmann constant.

Given the relation
<math display="block">F=U+TS</math>
where <math display="inline">U</math> is the internal energy of the system— for a crystal this is its electric potential energy
<math display="block">U_E = \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}}</math>

where <math display="inline">i</math> and <math display="inline">j</math> are the indices of the ions, <math display="inline">r_{ij}</math> is the distance between <math display="inline">i</math> and <math display="inline">j</math> and <math display="inline">\varepsilon_0 = 8.8542 \times 10^{-12} F\cdot m^{-1} </math>

The Helmholtz free energy of a crystal is thus given by

<math display="block">F= \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}} + \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu + k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math>

This equation could be used to qualitatively rationalise the free energy dependence on temperature. The data obtained is plotted in Figure ____.
[FIGURE]
Particularly, there are two salient regimes of interest. At low temperatures, T < 100 K, the graph is flat. However, at high temperatures, the behaviour is approximately linear. These observations are in agreement with the above equation, which highlights the temperature dependence of entropy <math display="inline">S</math>. At low temperatures, the term <math> k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> is extremely small, and hence the free energy term is dominated by the internal energy of the crystal. At high temperatures, the term <math>-k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> dominates and therefore the free energy of the system appears to have a dependence in temperature.

==== Variation of Lattice Parameter with Temperature ====
[[File:Syl815CellPvsT.png|thumb|center|600px|This figure illustrates the variation in cell parameter of MgO with temperature]]

As the temperature increases, the lattice parameter increases. It can thus be observed that the cell volume <math display="inline">V </math> has a dependence on temperature <math display="inline">T </math>, and the thermal expansion coefficient <math display="inline">\alpha </math> is given by

<math display="block">\alpha = \frac{1}{3V} \left( \frac{\partial V}{\partial T}\right)_P = \frac{1}{3B} \left( \frac{\partial P}{\partial T}\right)_V </math>
Where <math display="inline">B </math> is the bulk modulus and <math display="inline">P </math> is the pressure.
At 300 K, <math display="inline">\alpha = 2.06\times 10^{-5} K^{-1}</math>, compared to a literature value of <math display="inline">\alpha = 3.06\times 10^{-5} K^{-1}</math>

=== MD Simulation ===
The cell volume per formula unit of MgO was plotted against temperatures between 20 K to 4000 K.
[[File:Syl815MD.png|thumb|center|600px|"Experimental Data MD"]]
Under MD, the cell volume generally increases linearly with temperature throughout. By considering the mean kinetic energy of the crystal
<math display="block">\left \langle E_k \right \rangle = \frac{1}{2} M \sum_{i=1}^N v_i^2 = \frac{3}{2} Nk_BT_{MD} </math>
Where <math display="inline">\left \langle E_k \right \rangle </math> is the average kinetic energy of the atoms, <math display="inline">M</math> is the mass of the crystal lattice, and <math display="inline">v_i</math> represents the velocity of the atom <math display="inline">i</math>. It can be observed that the cell energy is linearly dependent on temperature. In a constant pressure system, this would result in volume expansion as temperature increases.

It can be observed that at high temperatures when <math display="inline">T\ge 2000 K </math>, more noise is present in the data due to the large cell volume and the large kinetic energy of the atoms.

[[File:Syl_MDvsQHA.png|thumb|center|600px|"This figure compares the data obtained for the thermal expansion of MgO under QHA and under MD."]]

At extremely low temperatures of <math display="inline">T\le 200 K </math>, QHA predicts a larger cell volume than MD. This can be attributed to the significant quantum effects at such low temperatures. Since MD only accounts for the kinetic energy of the atoms and neglects zero point vibrations, it predicts a smaller cell volume with the atoms closer together.

The data obtained for MD and QHA demonstrate strong agreement for temperatures between 200 to 1000 K. At these temperatures, the thermal energy of the system is sufficiently large such that the motion of the particles can be described classically.

Links:
https://pubs.acs.org/doi/pdf/10.1021/ct0500904
https://journals.aps.org/prb/pdf/10.1103/PhysRevB.43.5024

Conclusion

Between MD and Thermal expansion model
Compare quantitatively with values
Why you think it shouldn’t be linear

== References ==

Rep:MgO:syl815

2018-02-14T11:52:02Z

Syl815: /* Computing the Free Energy Using the Harmonic Approximation */

== Abstract ==
== Introduction ==
Thermal properties of materials are described by their vibrational free energies, which can be described in terms of the relative motion of atoms or the motion of their centre-of-mass. <ref name='phonons'>G. Srivastava, The physics of phonons, A. Hilger, Bristol, 1990.</ref> These concepts give rise to different approaches in calculating vibrational free energy, and both methodologies will be explored in greater detail.

=== Quasi-harmonic Approximation (QHA) ===
Fundamentally, QHA invokes the description of a crystalline solid as a primitive unit cell. This is essential due to the impracticality of calculating all the vibrational degrees of freedom in a crystal—for a crystal of size <math>N</math>, there are <math>3N</math> degrees of vibrational freedom, and in an infinitely large crystal lattice, 3N --> infinity.
Nonetheless, the translational periodicity of the crystal lattice, where <math display="inline">f(r + T) = f(r)</math>, simplifies the dynamics of all atoms in the lattice into that of a unit cell. For such a simplification to be appropriate, the following assumptions are made.

==== The Adiabatic Approximation ====
The adiabatic approximation separates the motion of the ion cores from that of the electrons since former are much more massive than the latter. Hence, the ion cores can be assumed to be in their equilibrium positions and that their motion is dependent on the potential field generated from the average motion of electrons.<ref name='phonons' />

==== The Harmonic Approximation ====
The total potential energy of a crystal can be expressed as the sum of all interatomic potentials. A two-body system typically has an anharmonic potential energy surface <math display="inline">U\left( r\right) </math>, where <math display="inline">r</math> is the interatomic separation. By considering a small displacement <math display="inline">x</math>,
<math display="block">x = r-r_0 \left( 1 \right)</math>

Where <math display="inline">r_0</math> is the equilibrium distance between the first and second atoms and is a minimum on the potential energy surface.

<math display="inline">U\left( r\right) </math> can be expanded in a Taylor series about <math display="inline">r_0</math>.
<math display="block">U\left( r\right) = U(r_0) + \frac{\partial U}{\partial x} x + \frac{\partial^2 U}{\partial x^2} x^2 + ... </math>
Since <math display="inline">U(r_0) </math> is unimportant in dynamics, <math display = "inline">\frac{\partial U}{\partial x}</math> is a force term and must be 0 for an equilibrium configuration, and all higher order terms <math display = "inline">x_n</math> , where <math display = "inline">n \ge 3</math> are assumed to be close to 0. As such, only the quadratic term is considered in the harmonic approximation. The solutions are the normal modes of vibrations for a system of independent quantum oscillators.

A phonon is a quantum of vibrational energy, hw, associated with a wave vector k.

Hence, for a crystal, its potential energy is given in the following equation.
[First principles, pdf, page 2]
Where l and k are the labels of the unit cells and atoms in each unit cell respectively. (http://dx.doi.org/10.1016/j.scriptamat.2015.07.021)

===== Limitations of Harmonic Approximation =====
The harmonic approximation predicts symmetric atomic vibrations about r0 at all temperatures, and is therefore incongruent with observed phenomena such as thermal expansion and heat conductivity. (https://core.ac.uk/download/pdf/42335965.pdf) The QHA causes renormalisation of the phonon frequencies and atomic force constants as is appropriate for the thermal equation of state. (https://ac.els-cdn.com/S0081194708606566/1-s2.0-S0081194708606566-main.pdf?_tid=f9617aac-104c-11e8-87bc-00000aab0f27&acdnat=1518478522_8ee2607e5ac970b2f00844b9ef59e8dc)

=== MD Simulation ===

=== System considered ===
MgO has a face-centred cubic lattice; more specifically, an alkali halide type structure with a Fm3m space group. (Citation)
MgO: insulator
Lattice parameters, symmetry, space group, etc.

=== Aims of the Exercise ===

== Methodology ==
Unless otherwise stated, all calculations were performed on a primitive unit cell of MgO with lattice parameters a = 2.9783 Angstrom, a = b = c; alpha = 60 degrees, and alpha = beta = gamma with GULP version 1.4.43 and crystals visualised with DLV interface.

A phonon dispersion curve was computed by sampling 100 points within the first Brillouin zone. The phonon density of states (DOS) was calculated with various shrinking factors, and the graphs subsequently plotted with matplotlib. The free energy of MgO was calculated with different shrinking factors at 300 K, and a suitable shrinking factor selected for the subsequent investigation of the thermal expansion of MgO. For every run, the Gibbs free energy was optimised, and calculations were performed from 0 to 2960 K in temperature steps of 20 K.

All MD simulations were performed on an isothermal-isobaric ensemble of MgO supercell of 32 formula units, with the following cell parameters:
A = 8.4239 Angstrom
A = b = c
Alpha = 90 degrees
Alpha = beta = gamma

MD was performed over a temperature range of 20 K to 4000 K, with temperature steps of 20 K. All calculations were performed with a time step of 1 fs. From 20 K to 1680 K, the system was allowed to first equilibrate for 1 ps; this was increased to 5 ps from 1700 K to 4000 K. Following which, MD production was allowed to run for 5 ps for all temperatures.

All data was analysed with Python on Jupyter notebook, and all graphs plotted with matplotlib.

== Results and Discussion ==
The lattice energy of MgO was calculated to be -41.0753 eV per primitive unit cell.
=== Phonon Modes of MgO ===
Figure 1 illustrates the phonon dispersion curve computed at 100 points for the primitive unit cell.

[FIGURE HERE]

Computation of the

A salient feature is the presence of 6 branches in the dispersion diagram. Assuming that the Born-von Karman boundary condition is satisfied, the edge effects of cells on dynamics can be ignored and u(N+1) = u1, where u is the displacement and N is the number of unit cells. This also implies the translational symmetry in k-space, such that all information of phonon dispersion can be derived by sampling in the first Brillouin zone.

By considering a linear diatomic chain satisfying the periodic boundary condition, the solutions to the vibrational frequency can be expressed in the form
[EQUATION, pg. 27 of book]

Highlighting two possible solutions for each k-value in a linear chain. Moreover, when m1 /= m2, a gap is observed at k = pi/2a, which is observed in Figure 1. (https://hal.archives-ouvertes.fr/file/index/docid/247340/filename/ajp-jp1v7p509.pdf)

By extending the logic to a 3D crystal lattice, the number of branches observed is given by 3x, where x is the number of atoms per unit cell. This is in agreement with the observation in Figure 1.

By appraising the solutions for k = 0 (long wavelength limit),

w1 = 2lambda(1/m + 1/M)
w2 = 0

w1 corresponds to a high energy mode where the atoms in the unit cell are moving out-of-phase, where frequency values are within the visible electromagnetic spectrum. The atoms are able to interact with an electric field of appropriate frequency due to the presence of both a positive and negative charge within the unit cell. It is hence naturally termed the optical mode. (Introduction to lattice dynamics)

On the other hand, w2 corresponds to a low energy mode with the atoms moving in phase and the wave pattern is similar to sound waves—hence the term acoustic mode. For any crystal with N atoms in the unit cell, there are only 3 acoustic—2 transverse and 1 longitudinal—and 3N-3 optical branches. The transverse modes are perpendicular to k, while the longitudinal mode is parallel.

For a cubic crystal, the highly symmetric nature indicates the possibility for some vibrations to be degenerate. Typically, transverse modes are lower in energy due to the weaker interaction between atoms in the unit cell.

This is encapsulated in Bloch's theorem and thus the wavefunction psi can be expressed in the form as given in equation 1.

[EQUATION HERE]

While typically applied to electrons in crystals, Bloch's theorem is broadly applicable in describing periodic wave phenomena, such as in phononic crystals. The branches relate to the modes of vibration, as given by equation for a one-dimensional system.

[EQUATION HERE]

Since the system is a three-dimensional one, vibrations can occur independently in perpendicular planes, therefore giving rise to additional normal coordinates and vibrations.

A dispersion diagram enables clear visualisation of the nature of the band gap, as illustrated in Figure ___. Figure __ clearly illustrates acoustic and optical phonons, depicting the in-phase and out-of-phase movement of the ions respectively.

Nonetheless, it is limited since energy values are unknown for points which are not sampled. Calculation of the density of states has been proven to be more useful.

=== Computing Density of States (DOS) ===
The impracticality of sampling all k-points within the FBZ can be circumvented by the use of a commensurate grid of k-points. To determine this set of k-points, the Pack-Monkhorst (PM) shrinking factor was used to specify the number of equidistant k-points taken along each direction of b1, b2 and b3 in one reciprocal lattice PUC. (https://books.google.co.uk/books?id=nX_wG7WaDJsC&pg=PA38&lpg=PA38&dq=Pack-Monkhorst+shrinking+factor&source=bl&ots=vL_-nToT5e&sig=SOp4EsY7oG-ki9tlvsMSRTJ-eTY&hl=en&sa=X&ved=0ahUKEwiG-PHajZfZAhXmLsAKHefcBDAQ6AEINzAC#v=onepage&q=Pack-Monkhorst%20shrinking%20factor&f=false) The Cartesian coordinates of the k-points calculated are given by the equation

[EQUATION FROM https://journals.aps.org/prb/pdf/10.1103/PhysRevB.93.155109]

A major advantage is its computational efficiency by restricting the number of k-points calculated to a finite value. Moreover, the accuracy obtained from calculations with a PUC can be comparable to that of a supercell as long as the shrinking factor is appropriate.

Table 1 illustrates the effect of modifying the PM shrinking factor on the number of k-points calculated.

1x1x1 1
2x2x2 4
3x3x3 18
4x4x4 32
5x5x5 75
6x6x6 108
8x8x8 256
10x10x10 500
16x16x16 2048
20x20x20 4000
64x64x64 >99999

As the mesh of k-points increases, the number of k-points calculated increases as well. This is contrary to the prediction from the above equation, where we would expect kx * ky* kz number of points. This can be attributed to the mapping of equivalent k-points onto each other and thus the number of k-points calculated is reduced.

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[GRAPH OF ENERGY AGAINST K-POINTS]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication, and will be further discussed in section 3.

An initial plot of the density of states was obtained from a 1x1x1 grid yielding six resultant modes. Sharp and distinct peaks are observed in the plot, since only one k-point was sampled.

Notably, only four unique peaks are observed even though we should observe 6 modes of vibrations. The final two modes are degenerate at _____ and _____ wavenumbers. Compared to the non-degenerate acoustic and optical peaks (___ and ____ respectively), the degenerate acoustic modes are higher in energy whereas the degenerate optical modes are lower in energy correspondingly. It can therefore be deduced that the degenerate acoustic and optical modes are transverse in nature.

The k-point used in the DOS calculation could be identified by comparing with the dispersion curve. Since point M contains all of the frequency values in Figure ___, it can be determined that the point represented in the DOS curve is M, where kx = 0.5, ky = 0.5 and kz = 0.5.

=== Relationship between the Dispersion Curve and DOS ===
The DOS curve illustrates the number of energy states per unit energy, demonstrating a mode at 414 cm-1. This correlates well with Figure ____. By constructing a horizontal line at frequency = 414 cm-1, it can be observed that the branches intersect this line frequently. This implies that a significant proportion of k-points have vibrational modes of frequency 414 cm-1.

Loss of information in dispersion curve: how is that so?
Orthogonal to the dispersion curve: gives the number of energy states.
Information provided for the dispersion curve vs the density of states.
=== Computing the Free Energy Using the Harmonic Approximation ===
The figure below demonstrates the relationship between the PM shrinking factor used and the computed Helmholtz free energy of the system.

[[file:Syl815FEvsT.png|thumb|center|600px|Free Energy vs. PM shrinking factor.]]

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[[file:Syl815FEshinking.png|thumb|center|600px|DOS against wavenumber for various shrinking factors.]]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication.

From Figure ____, the free energy of MgO is observed to increase and converge to a value of -40.926 483 eV, and it is observed that this occurs for a grid size of 8x8x8. (Why does it do so? Explain here)

A 2x2x2 grid is sufficient for calculating the free energy of MgO to 1 meV. A 4x4x4 grid is necessary for a precision to 0.5 meV and 0.1 meV.

=== Thermal Expansion ===
The Helmholtz free energy of a crystal is given by the sum of the energies of independent vibrational waves. The energy level <math display="inline">n</math> of a quantum harmonic oscillator are given by:

<math display="block">E_n = \left( n+ \frac{1}{2} \right) h \nu </math>

where <math display="inline">h</math> is Planck's constant and <math display="inline">\nu</math> is the frequency of energy level <math display="inline">n</math>. For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators, the vibrational energy is given by:

<math display="block">E_{vib} = \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu</math>

For a canonical <math display="inline">NVT</math> ensemble, the partition function is given by

<math display="block">Z = \sum_n e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

where <math display="inline">\beta = \frac{1}{kT}</math> and <math display="inline">E_n</math> enumerates all vibrational energy states.

For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators,

<math display="block">Z_N = \prod_n^{3N} e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

The phonon entropy can then be expressed in terms of the partition function:
<math display="block">S=-k_B \ln Z_N</math> where <math display="inline">k_B</math> is the Boltzmann constant.

Given the relation
<math display="block">F=U+TS</math>
where <math display="inline">U</math> is the internal energy of the system— for a crystal this is its electric potential energy
<math display="block">U_E = \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}}</math>

where <math display="inline">i</math> and <math display="inline">j</math> are the indices of the ions, <math display="inline">r_{ij}</math> is the distance between <math display="inline">i</math> and <math display="inline">j</math> and <math display="inline">\varepsilon_0 = 8.8542 \times 10^{-12} F\cdot m^{-1} </math>

The Helmholtz free energy of a crystal is thus given by

<math display="block">F= \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}} + \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu + k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math>

This equation could be used to qualitatively rationalise the free energy dependence on temperature. The data obtained is plotted in Figure ____.
[FIGURE]
Particularly, there are two salient regimes of interest. At low temperatures, T < 100 K, the graph is flat. However, at high temperatures, the behaviour is approximately linear. These observations are in agreement with the above equation, which highlights the temperature dependence of entropy <math display="inline">S</math>. At low temperatures, the term <math> k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> is extremely small, and hence the free energy term is dominated by the internal energy of the crystal. At high temperatures, the term <math>-k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> dominates and therefore the free energy of the system appears to have a dependence in temperature.

==== Variation of Lattice Parameter with Temperature ====
[[File:Syl815CellPvsT.png|thumb|center|600px|This figure illustrates the variation in cell parameter of MgO with temperature]]

As the temperature increases, the lattice parameter increases. It can thus be observed that the cell volume <math display="inline">V </math> has a dependence on temperature <math display="inline">T </math>, and the thermal expansion coefficient <math display="inline">\alpha </math> is given by

<math display="block">\alpha = \frac{1}{3V} \left( \frac{\partial V}{\partial T}\right)_P = \frac{1}{3B} \left( \frac{\partial P}{\partial T}\right)_V </math>
Where <math display="inline">B </math> is the bulk modulus and <math display="inline">P </math> is the pressure.
At 300 K, <math display="inline">\alpha = 2.06\times 10^{-5} K^{-1}</math>, compared to a literature value of <math display="inline">\alpha = 3.06\times 10^{-5} K^{-1}</math>

=== MD Simulation ===
The cell volume per formula unit of MgO was plotted against temperatures between 20 K to 4000 K.
[[File:Syl815MD.png|thumb|center|600px|"Experimental Data MD"]]
Under MD, the cell volume generally increases linearly with temperature throughout. By considering the mean kinetic energy of the crystal
<math display="block">\left \langle E_k \right \rangle = \frac{1}{2} M \sum_{i=1}^N v_i^2 = \frac{3}{2} Nk_BT_{MD} </math>
Where <math display="inline">\left \langle E_k \right \rangle </math> is the average kinetic energy of the atoms, <math display="inline">M</math> is the mass of the crystal lattice, and <math display="inline">v_i</math> represents the velocity of the atom <math display="inline">i</math>. It can be observed that the cell energy is linearly dependent on temperature. In a constant pressure system, this would result in volume expansion as temperature increases.

It can be observed that at high temperatures when <math display="inline">T\ge 2000 K </math>, more noise is present in the data due to the large cell volume and the large kinetic energy of the atoms.

[[File:Syl_MDvsQHA.png|thumb|center|600px|"This figure compares the data obtained for the thermal expansion of MgO under QHA and under MD."]]

At extremely low temperatures of <math display="inline">T\le 200 K </math>, QHA predicts a larger cell volume than MD. This can be attributed to the significant quantum effects at such low temperatures. Since MD only accounts for the kinetic energy of the atoms and neglects zero point vibrations, it predicts a smaller cell volume with the atoms closer together.

The data obtained for MD and QHA demonstrate strong agreement for temperatures between 200 to 1000 K. At these temperatures, the thermal energy of the system is sufficiently large such that the motion of the particles can be described classically.

Links:
https://pubs.acs.org/doi/pdf/10.1021/ct0500904
https://journals.aps.org/prb/pdf/10.1103/PhysRevB.43.5024

Conclusion

Between MD and Thermal expansion model
Compare quantitatively with values
Why you think it shouldn’t be linear

== References ==

File:Syl815FEshinking.png

2018-02-14T11:51:11Z

Syl815:

File:Syl815FEvsT.png

2018-02-14T11:49:39Z

Syl815:

Rep:MgO:syl815

2018-02-14T11:49:12Z

Syl815: /* Variation of Lattice Parameter with Temperature */

== Abstract ==
== Introduction ==
Thermal properties of materials are described by their vibrational free energies, which can be described in terms of the relative motion of atoms or the motion of their centre-of-mass. <ref name='phonons'>G. Srivastava, The physics of phonons, A. Hilger, Bristol, 1990.</ref> These concepts give rise to different approaches in calculating vibrational free energy, and both methodologies will be explored in greater detail.

=== Quasi-harmonic Approximation (QHA) ===
Fundamentally, QHA invokes the description of a crystalline solid as a primitive unit cell. This is essential due to the impracticality of calculating all the vibrational degrees of freedom in a crystal—for a crystal of size <math>N</math>, there are <math>3N</math> degrees of vibrational freedom, and in an infinitely large crystal lattice, 3N --> infinity.
Nonetheless, the translational periodicity of the crystal lattice, where <math display="inline">f(r + T) = f(r)</math>, simplifies the dynamics of all atoms in the lattice into that of a unit cell. For such a simplification to be appropriate, the following assumptions are made.

==== The Adiabatic Approximation ====
The adiabatic approximation separates the motion of the ion cores from that of the electrons since former are much more massive than the latter. Hence, the ion cores can be assumed to be in their equilibrium positions and that their motion is dependent on the potential field generated from the average motion of electrons.<ref name='phonons' />

==== The Harmonic Approximation ====
The total potential energy of a crystal can be expressed as the sum of all interatomic potentials. A two-body system typically has an anharmonic potential energy surface <math display="inline">U\left( r\right) </math>, where <math display="inline">r</math> is the interatomic separation. By considering a small displacement <math display="inline">x</math>,
<math display="block">x = r-r_0 \left( 1 \right)</math>

Where <math display="inline">r_0</math> is the equilibrium distance between the first and second atoms and is a minimum on the potential energy surface.

<math display="inline">U\left( r\right) </math> can be expanded in a Taylor series about <math display="inline">r_0</math>.
<math display="block">U\left( r\right) = U(r_0) + \frac{\partial U}{\partial x} x + \frac{\partial^2 U}{\partial x^2} x^2 + ... </math>
Since <math display="inline">U(r_0) </math> is unimportant in dynamics, <math display = "inline">\frac{\partial U}{\partial x}</math> is a force term and must be 0 for an equilibrium configuration, and all higher order terms <math display = "inline">x_n</math> , where <math display = "inline">n \ge 3</math> are assumed to be close to 0. As such, only the quadratic term is considered in the harmonic approximation. The solutions are the normal modes of vibrations for a system of independent quantum oscillators.

A phonon is a quantum of vibrational energy, hw, associated with a wave vector k.

Hence, for a crystal, its potential energy is given in the following equation.
[First principles, pdf, page 2]
Where l and k are the labels of the unit cells and atoms in each unit cell respectively. (http://dx.doi.org/10.1016/j.scriptamat.2015.07.021)

===== Limitations of Harmonic Approximation =====
The harmonic approximation predicts symmetric atomic vibrations about r0 at all temperatures, and is therefore incongruent with observed phenomena such as thermal expansion and heat conductivity. (https://core.ac.uk/download/pdf/42335965.pdf) The QHA causes renormalisation of the phonon frequencies and atomic force constants as is appropriate for the thermal equation of state. (https://ac.els-cdn.com/S0081194708606566/1-s2.0-S0081194708606566-main.pdf?_tid=f9617aac-104c-11e8-87bc-00000aab0f27&acdnat=1518478522_8ee2607e5ac970b2f00844b9ef59e8dc)

=== MD Simulation ===

=== System considered ===
MgO has a face-centred cubic lattice; more specifically, an alkali halide type structure with a Fm3m space group. (Citation)
MgO: insulator
Lattice parameters, symmetry, space group, etc.

=== Aims of the Exercise ===

== Methodology ==
Unless otherwise stated, all calculations were performed on a primitive unit cell of MgO with lattice parameters a = 2.9783 Angstrom, a = b = c; alpha = 60 degrees, and alpha = beta = gamma with GULP version 1.4.43 and crystals visualised with DLV interface.

A phonon dispersion curve was computed by sampling 100 points within the first Brillouin zone. The phonon density of states (DOS) was calculated with various shrinking factors, and the graphs subsequently plotted with matplotlib. The free energy of MgO was calculated with different shrinking factors at 300 K, and a suitable shrinking factor selected for the subsequent investigation of the thermal expansion of MgO. For every run, the Gibbs free energy was optimised, and calculations were performed from 0 to 2960 K in temperature steps of 20 K.

All MD simulations were performed on an isothermal-isobaric ensemble of MgO supercell of 32 formula units, with the following cell parameters:
A = 8.4239 Angstrom
A = b = c
Alpha = 90 degrees
Alpha = beta = gamma

MD was performed over a temperature range of 20 K to 4000 K, with temperature steps of 20 K. All calculations were performed with a time step of 1 fs. From 20 K to 1680 K, the system was allowed to first equilibrate for 1 ps; this was increased to 5 ps from 1700 K to 4000 K. Following which, MD production was allowed to run for 5 ps for all temperatures.

All data was analysed with Python on Jupyter notebook, and all graphs plotted with matplotlib.

== Results and Discussion ==
The lattice energy of MgO was calculated to be -41.0753 eV per primitive unit cell.
=== Phonon Modes of MgO ===
Figure 1 illustrates the phonon dispersion curve computed at 100 points for the primitive unit cell.

[FIGURE HERE]

Computation of the

A salient feature is the presence of 6 branches in the dispersion diagram. Assuming that the Born-von Karman boundary condition is satisfied, the edge effects of cells on dynamics can be ignored and u(N+1) = u1, where u is the displacement and N is the number of unit cells. This also implies the translational symmetry in k-space, such that all information of phonon dispersion can be derived by sampling in the first Brillouin zone.

By considering a linear diatomic chain satisfying the periodic boundary condition, the solutions to the vibrational frequency can be expressed in the form
[EQUATION, pg. 27 of book]

Highlighting two possible solutions for each k-value in a linear chain. Moreover, when m1 /= m2, a gap is observed at k = pi/2a, which is observed in Figure 1. (https://hal.archives-ouvertes.fr/file/index/docid/247340/filename/ajp-jp1v7p509.pdf)

By extending the logic to a 3D crystal lattice, the number of branches observed is given by 3x, where x is the number of atoms per unit cell. This is in agreement with the observation in Figure 1.

By appraising the solutions for k = 0 (long wavelength limit),

w1 = 2lambda(1/m + 1/M)
w2 = 0

w1 corresponds to a high energy mode where the atoms in the unit cell are moving out-of-phase, where frequency values are within the visible electromagnetic spectrum. The atoms are able to interact with an electric field of appropriate frequency due to the presence of both a positive and negative charge within the unit cell. It is hence naturally termed the optical mode. (Introduction to lattice dynamics)

On the other hand, w2 corresponds to a low energy mode with the atoms moving in phase and the wave pattern is similar to sound waves—hence the term acoustic mode. For any crystal with N atoms in the unit cell, there are only 3 acoustic—2 transverse and 1 longitudinal—and 3N-3 optical branches. The transverse modes are perpendicular to k, while the longitudinal mode is parallel.

For a cubic crystal, the highly symmetric nature indicates the possibility for some vibrations to be degenerate. Typically, transverse modes are lower in energy due to the weaker interaction between atoms in the unit cell.

This is encapsulated in Bloch's theorem and thus the wavefunction psi can be expressed in the form as given in equation 1.

[EQUATION HERE]

While typically applied to electrons in crystals, Bloch's theorem is broadly applicable in describing periodic wave phenomena, such as in phononic crystals. The branches relate to the modes of vibration, as given by equation for a one-dimensional system.

[EQUATION HERE]

Since the system is a three-dimensional one, vibrations can occur independently in perpendicular planes, therefore giving rise to additional normal coordinates and vibrations.

A dispersion diagram enables clear visualisation of the nature of the band gap, as illustrated in Figure ___. Figure __ clearly illustrates acoustic and optical phonons, depicting the in-phase and out-of-phase movement of the ions respectively.

Nonetheless, it is limited since energy values are unknown for points which are not sampled. Calculation of the density of states has been proven to be more useful.

=== Computing Density of States (DOS) ===
The impracticality of sampling all k-points within the FBZ can be circumvented by the use of a commensurate grid of k-points. To determine this set of k-points, the Pack-Monkhorst (PM) shrinking factor was used to specify the number of equidistant k-points taken along each direction of b1, b2 and b3 in one reciprocal lattice PUC. (https://books.google.co.uk/books?id=nX_wG7WaDJsC&pg=PA38&lpg=PA38&dq=Pack-Monkhorst+shrinking+factor&source=bl&ots=vL_-nToT5e&sig=SOp4EsY7oG-ki9tlvsMSRTJ-eTY&hl=en&sa=X&ved=0ahUKEwiG-PHajZfZAhXmLsAKHefcBDAQ6AEINzAC#v=onepage&q=Pack-Monkhorst%20shrinking%20factor&f=false) The Cartesian coordinates of the k-points calculated are given by the equation

[EQUATION FROM https://journals.aps.org/prb/pdf/10.1103/PhysRevB.93.155109]

A major advantage is its computational efficiency by restricting the number of k-points calculated to a finite value. Moreover, the accuracy obtained from calculations with a PUC can be comparable to that of a supercell as long as the shrinking factor is appropriate.

Table 1 illustrates the effect of modifying the PM shrinking factor on the number of k-points calculated.

1x1x1 1
2x2x2 4
3x3x3 18
4x4x4 32
5x5x5 75
6x6x6 108
8x8x8 256
10x10x10 500
16x16x16 2048
20x20x20 4000
64x64x64 >99999

As the mesh of k-points increases, the number of k-points calculated increases as well. This is contrary to the prediction from the above equation, where we would expect kx * ky* kz number of points. This can be attributed to the mapping of equivalent k-points onto each other and thus the number of k-points calculated is reduced.

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[GRAPH OF ENERGY AGAINST K-POINTS]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication, and will be further discussed in section 3.

An initial plot of the density of states was obtained from a 1x1x1 grid yielding six resultant modes. Sharp and distinct peaks are observed in the plot, since only one k-point was sampled.

Notably, only four unique peaks are observed even though we should observe 6 modes of vibrations. The final two modes are degenerate at _____ and _____ wavenumbers. Compared to the non-degenerate acoustic and optical peaks (___ and ____ respectively), the degenerate acoustic modes are higher in energy whereas the degenerate optical modes are lower in energy correspondingly. It can therefore be deduced that the degenerate acoustic and optical modes are transverse in nature.

The k-point used in the DOS calculation could be identified by comparing with the dispersion curve. Since point M contains all of the frequency values in Figure ___, it can be determined that the point represented in the DOS curve is M, where kx = 0.5, ky = 0.5 and kz = 0.5.

=== Relationship between the Dispersion Curve and DOS ===
The DOS curve illustrates the number of energy states per unit energy, demonstrating a mode at 414 cm-1. This correlates well with Figure ____. By constructing a horizontal line at frequency = 414 cm-1, it can be observed that the branches intersect this line frequently. This implies that a significant proportion of k-points have vibrational modes of frequency 414 cm-1.

Loss of information in dispersion curve: how is that so?
Orthogonal to the dispersion curve: gives the number of energy states.
Information provided for the dispersion curve vs the density of states.
=== Computing the Free Energy Using the Harmonic Approximation ===
Figure ___ demonstrates the relationship between the PM shrinking factor used and the computed Helmholtz free energy of the system.
[FIGURE HERE]

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[GRAPH OF ENERGY AGAINST K-POINTS]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication.

From Figure ____, the free energy of MgO is observed to increase and converge to a value of -40.926 483 eV, and it is observed that this occurs for a grid size of 8x8x8. (Why does it do so? Explain here)

A 2x2x2 grid is sufficient for calculating the free energy of MgO to 1 meV. A 4x4x4 grid is necessary for a precision to 0.5 meV and 0.1 meV.

=== Thermal Expansion ===
The Helmholtz free energy of a crystal is given by the sum of the energies of independent vibrational waves. The energy level <math display="inline">n</math> of a quantum harmonic oscillator are given by:

<math display="block">E_n = \left( n+ \frac{1}{2} \right) h \nu </math>

where <math display="inline">h</math> is Planck's constant and <math display="inline">\nu</math> is the frequency of energy level <math display="inline">n</math>. For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators, the vibrational energy is given by:

<math display="block">E_{vib} = \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu</math>

For a canonical <math display="inline">NVT</math> ensemble, the partition function is given by

<math display="block">Z = \sum_n e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

where <math display="inline">\beta = \frac{1}{kT}</math> and <math display="inline">E_n</math> enumerates all vibrational energy states.

For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators,

<math display="block">Z_N = \prod_n^{3N} e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

The phonon entropy can then be expressed in terms of the partition function:
<math display="block">S=-k_B \ln Z_N</math> where <math display="inline">k_B</math> is the Boltzmann constant.

Given the relation
<math display="block">F=U+TS</math>
where <math display="inline">U</math> is the internal energy of the system— for a crystal this is its electric potential energy
<math display="block">U_E = \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}}</math>

where <math display="inline">i</math> and <math display="inline">j</math> are the indices of the ions, <math display="inline">r_{ij}</math> is the distance between <math display="inline">i</math> and <math display="inline">j</math> and <math display="inline">\varepsilon_0 = 8.8542 \times 10^{-12} F\cdot m^{-1} </math>

The Helmholtz free energy of a crystal is thus given by

<math display="block">F= \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}} + \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu + k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math>

This equation could be used to qualitatively rationalise the free energy dependence on temperature. The data obtained is plotted in Figure ____.
[FIGURE]
Particularly, there are two salient regimes of interest. At low temperatures, T < 100 K, the graph is flat. However, at high temperatures, the behaviour is approximately linear. These observations are in agreement with the above equation, which highlights the temperature dependence of entropy <math display="inline">S</math>. At low temperatures, the term <math> k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> is extremely small, and hence the free energy term is dominated by the internal energy of the crystal. At high temperatures, the term <math>-k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> dominates and therefore the free energy of the system appears to have a dependence in temperature.

==== Variation of Lattice Parameter with Temperature ====
[[File:Syl815CellPvsT.png|thumb|center|600px|This figure illustrates the variation in cell parameter of MgO with temperature]]

As the temperature increases, the lattice parameter increases. It can thus be observed that the cell volume <math display="inline">V </math> has a dependence on temperature <math display="inline">T </math>, and the thermal expansion coefficient <math display="inline">\alpha </math> is given by

<math display="block">\alpha = \frac{1}{3V} \left( \frac{\partial V}{\partial T}\right)_P = \frac{1}{3B} \left( \frac{\partial P}{\partial T}\right)_V </math>
Where <math display="inline">B </math> is the bulk modulus and <math display="inline">P </math> is the pressure.
At 300 K, <math display="inline">\alpha = 2.06\times 10^{-5} K^{-1}</math>, compared to a literature value of <math display="inline">\alpha = 3.06\times 10^{-5} K^{-1}</math>

=== MD Simulation ===
The cell volume per formula unit of MgO was plotted against temperatures between 20 K to 4000 K.
[[File:Syl815MD.png|thumb|center|600px|"Experimental Data MD"]]
Under MD, the cell volume generally increases linearly with temperature throughout. By considering the mean kinetic energy of the crystal
<math display="block">\left \langle E_k \right \rangle = \frac{1}{2} M \sum_{i=1}^N v_i^2 = \frac{3}{2} Nk_BT_{MD} </math>
Where <math display="inline">\left \langle E_k \right \rangle </math> is the average kinetic energy of the atoms, <math display="inline">M</math> is the mass of the crystal lattice, and <math display="inline">v_i</math> represents the velocity of the atom <math display="inline">i</math>. It can be observed that the cell energy is linearly dependent on temperature. In a constant pressure system, this would result in volume expansion as temperature increases.

It can be observed that at high temperatures when <math display="inline">T\ge 2000 K </math>, more noise is present in the data due to the large cell volume and the large kinetic energy of the atoms.

[[File:Syl_MDvsQHA.png|thumb|center|600px|"This figure compares the data obtained for the thermal expansion of MgO under QHA and under MD."]]

At extremely low temperatures of <math display="inline">T\le 200 K </math>, QHA predicts a larger cell volume than MD. This can be attributed to the significant quantum effects at such low temperatures. Since MD only accounts for the kinetic energy of the atoms and neglects zero point vibrations, it predicts a smaller cell volume with the atoms closer together.

The data obtained for MD and QHA demonstrate strong agreement for temperatures between 200 to 1000 K. At these temperatures, the thermal energy of the system is sufficiently large such that the motion of the particles can be described classically.

Links:
https://pubs.acs.org/doi/pdf/10.1021/ct0500904
https://journals.aps.org/prb/pdf/10.1103/PhysRevB.43.5024

Conclusion

Between MD and Thermal expansion model
Compare quantitatively with values
Why you think it shouldn’t be linear

== References ==

Rep:MgO:syl815

2018-02-14T11:48:29Z

Syl815: /* Variation of Lattice Parameter with Temperature */

== Abstract ==
== Introduction ==
Thermal properties of materials are described by their vibrational free energies, which can be described in terms of the relative motion of atoms or the motion of their centre-of-mass. <ref name='phonons'>G. Srivastava, The physics of phonons, A. Hilger, Bristol, 1990.</ref> These concepts give rise to different approaches in calculating vibrational free energy, and both methodologies will be explored in greater detail.

=== Quasi-harmonic Approximation (QHA) ===
Fundamentally, QHA invokes the description of a crystalline solid as a primitive unit cell. This is essential due to the impracticality of calculating all the vibrational degrees of freedom in a crystal—for a crystal of size <math>N</math>, there are <math>3N</math> degrees of vibrational freedom, and in an infinitely large crystal lattice, 3N --> infinity.
Nonetheless, the translational periodicity of the crystal lattice, where <math display="inline">f(r + T) = f(r)</math>, simplifies the dynamics of all atoms in the lattice into that of a unit cell. For such a simplification to be appropriate, the following assumptions are made.

==== The Adiabatic Approximation ====
The adiabatic approximation separates the motion of the ion cores from that of the electrons since former are much more massive than the latter. Hence, the ion cores can be assumed to be in their equilibrium positions and that their motion is dependent on the potential field generated from the average motion of electrons.<ref name='phonons' />

==== The Harmonic Approximation ====
The total potential energy of a crystal can be expressed as the sum of all interatomic potentials. A two-body system typically has an anharmonic potential energy surface <math display="inline">U\left( r\right) </math>, where <math display="inline">r</math> is the interatomic separation. By considering a small displacement <math display="inline">x</math>,
<math display="block">x = r-r_0 \left( 1 \right)</math>

Where <math display="inline">r_0</math> is the equilibrium distance between the first and second atoms and is a minimum on the potential energy surface.

<math display="inline">U\left( r\right) </math> can be expanded in a Taylor series about <math display="inline">r_0</math>.
<math display="block">U\left( r\right) = U(r_0) + \frac{\partial U}{\partial x} x + \frac{\partial^2 U}{\partial x^2} x^2 + ... </math>
Since <math display="inline">U(r_0) </math> is unimportant in dynamics, <math display = "inline">\frac{\partial U}{\partial x}</math> is a force term and must be 0 for an equilibrium configuration, and all higher order terms <math display = "inline">x_n</math> , where <math display = "inline">n \ge 3</math> are assumed to be close to 0. As such, only the quadratic term is considered in the harmonic approximation. The solutions are the normal modes of vibrations for a system of independent quantum oscillators.

A phonon is a quantum of vibrational energy, hw, associated with a wave vector k.

Hence, for a crystal, its potential energy is given in the following equation.
[First principles, pdf, page 2]
Where l and k are the labels of the unit cells and atoms in each unit cell respectively. (http://dx.doi.org/10.1016/j.scriptamat.2015.07.021)

===== Limitations of Harmonic Approximation =====
The harmonic approximation predicts symmetric atomic vibrations about r0 at all temperatures, and is therefore incongruent with observed phenomena such as thermal expansion and heat conductivity. (https://core.ac.uk/download/pdf/42335965.pdf) The QHA causes renormalisation of the phonon frequencies and atomic force constants as is appropriate for the thermal equation of state. (https://ac.els-cdn.com/S0081194708606566/1-s2.0-S0081194708606566-main.pdf?_tid=f9617aac-104c-11e8-87bc-00000aab0f27&acdnat=1518478522_8ee2607e5ac970b2f00844b9ef59e8dc)

=== MD Simulation ===

=== System considered ===
MgO has a face-centred cubic lattice; more specifically, an alkali halide type structure with a Fm3m space group. (Citation)
MgO: insulator
Lattice parameters, symmetry, space group, etc.

=== Aims of the Exercise ===

== Methodology ==
Unless otherwise stated, all calculations were performed on a primitive unit cell of MgO with lattice parameters a = 2.9783 Angstrom, a = b = c; alpha = 60 degrees, and alpha = beta = gamma with GULP version 1.4.43 and crystals visualised with DLV interface.

A phonon dispersion curve was computed by sampling 100 points within the first Brillouin zone. The phonon density of states (DOS) was calculated with various shrinking factors, and the graphs subsequently plotted with matplotlib. The free energy of MgO was calculated with different shrinking factors at 300 K, and a suitable shrinking factor selected for the subsequent investigation of the thermal expansion of MgO. For every run, the Gibbs free energy was optimised, and calculations were performed from 0 to 2960 K in temperature steps of 20 K.

All MD simulations were performed on an isothermal-isobaric ensemble of MgO supercell of 32 formula units, with the following cell parameters:
A = 8.4239 Angstrom
A = b = c
Alpha = 90 degrees
Alpha = beta = gamma

MD was performed over a temperature range of 20 K to 4000 K, with temperature steps of 20 K. All calculations were performed with a time step of 1 fs. From 20 K to 1680 K, the system was allowed to first equilibrate for 1 ps; this was increased to 5 ps from 1700 K to 4000 K. Following which, MD production was allowed to run for 5 ps for all temperatures.

All data was analysed with Python on Jupyter notebook, and all graphs plotted with matplotlib.

== Results and Discussion ==
The lattice energy of MgO was calculated to be -41.0753 eV per primitive unit cell.
=== Phonon Modes of MgO ===
Figure 1 illustrates the phonon dispersion curve computed at 100 points for the primitive unit cell.

[FIGURE HERE]

Computation of the

A salient feature is the presence of 6 branches in the dispersion diagram. Assuming that the Born-von Karman boundary condition is satisfied, the edge effects of cells on dynamics can be ignored and u(N+1) = u1, where u is the displacement and N is the number of unit cells. This also implies the translational symmetry in k-space, such that all information of phonon dispersion can be derived by sampling in the first Brillouin zone.

By considering a linear diatomic chain satisfying the periodic boundary condition, the solutions to the vibrational frequency can be expressed in the form
[EQUATION, pg. 27 of book]

Highlighting two possible solutions for each k-value in a linear chain. Moreover, when m1 /= m2, a gap is observed at k = pi/2a, which is observed in Figure 1. (https://hal.archives-ouvertes.fr/file/index/docid/247340/filename/ajp-jp1v7p509.pdf)

By extending the logic to a 3D crystal lattice, the number of branches observed is given by 3x, where x is the number of atoms per unit cell. This is in agreement with the observation in Figure 1.

By appraising the solutions for k = 0 (long wavelength limit),

w1 = 2lambda(1/m + 1/M)
w2 = 0

w1 corresponds to a high energy mode where the atoms in the unit cell are moving out-of-phase, where frequency values are within the visible electromagnetic spectrum. The atoms are able to interact with an electric field of appropriate frequency due to the presence of both a positive and negative charge within the unit cell. It is hence naturally termed the optical mode. (Introduction to lattice dynamics)

On the other hand, w2 corresponds to a low energy mode with the atoms moving in phase and the wave pattern is similar to sound waves—hence the term acoustic mode. For any crystal with N atoms in the unit cell, there are only 3 acoustic—2 transverse and 1 longitudinal—and 3N-3 optical branches. The transverse modes are perpendicular to k, while the longitudinal mode is parallel.

For a cubic crystal, the highly symmetric nature indicates the possibility for some vibrations to be degenerate. Typically, transverse modes are lower in energy due to the weaker interaction between atoms in the unit cell.

This is encapsulated in Bloch's theorem and thus the wavefunction psi can be expressed in the form as given in equation 1.

[EQUATION HERE]

While typically applied to electrons in crystals, Bloch's theorem is broadly applicable in describing periodic wave phenomena, such as in phononic crystals. The branches relate to the modes of vibration, as given by equation for a one-dimensional system.

[EQUATION HERE]

Since the system is a three-dimensional one, vibrations can occur independently in perpendicular planes, therefore giving rise to additional normal coordinates and vibrations.

A dispersion diagram enables clear visualisation of the nature of the band gap, as illustrated in Figure ___. Figure __ clearly illustrates acoustic and optical phonons, depicting the in-phase and out-of-phase movement of the ions respectively.

Nonetheless, it is limited since energy values are unknown for points which are not sampled. Calculation of the density of states has been proven to be more useful.

=== Computing Density of States (DOS) ===
The impracticality of sampling all k-points within the FBZ can be circumvented by the use of a commensurate grid of k-points. To determine this set of k-points, the Pack-Monkhorst (PM) shrinking factor was used to specify the number of equidistant k-points taken along each direction of b1, b2 and b3 in one reciprocal lattice PUC. (https://books.google.co.uk/books?id=nX_wG7WaDJsC&pg=PA38&lpg=PA38&dq=Pack-Monkhorst+shrinking+factor&source=bl&ots=vL_-nToT5e&sig=SOp4EsY7oG-ki9tlvsMSRTJ-eTY&hl=en&sa=X&ved=0ahUKEwiG-PHajZfZAhXmLsAKHefcBDAQ6AEINzAC#v=onepage&q=Pack-Monkhorst%20shrinking%20factor&f=false) The Cartesian coordinates of the k-points calculated are given by the equation

[EQUATION FROM https://journals.aps.org/prb/pdf/10.1103/PhysRevB.93.155109]

A major advantage is its computational efficiency by restricting the number of k-points calculated to a finite value. Moreover, the accuracy obtained from calculations with a PUC can be comparable to that of a supercell as long as the shrinking factor is appropriate.

Table 1 illustrates the effect of modifying the PM shrinking factor on the number of k-points calculated.

1x1x1 1
2x2x2 4
3x3x3 18
4x4x4 32
5x5x5 75
6x6x6 108
8x8x8 256
10x10x10 500
16x16x16 2048
20x20x20 4000
64x64x64 >99999

As the mesh of k-points increases, the number of k-points calculated increases as well. This is contrary to the prediction from the above equation, where we would expect kx * ky* kz number of points. This can be attributed to the mapping of equivalent k-points onto each other and thus the number of k-points calculated is reduced.

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[GRAPH OF ENERGY AGAINST K-POINTS]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication, and will be further discussed in section 3.

An initial plot of the density of states was obtained from a 1x1x1 grid yielding six resultant modes. Sharp and distinct peaks are observed in the plot, since only one k-point was sampled.

Notably, only four unique peaks are observed even though we should observe 6 modes of vibrations. The final two modes are degenerate at _____ and _____ wavenumbers. Compared to the non-degenerate acoustic and optical peaks (___ and ____ respectively), the degenerate acoustic modes are higher in energy whereas the degenerate optical modes are lower in energy correspondingly. It can therefore be deduced that the degenerate acoustic and optical modes are transverse in nature.

The k-point used in the DOS calculation could be identified by comparing with the dispersion curve. Since point M contains all of the frequency values in Figure ___, it can be determined that the point represented in the DOS curve is M, where kx = 0.5, ky = 0.5 and kz = 0.5.

=== Relationship between the Dispersion Curve and DOS ===
The DOS curve illustrates the number of energy states per unit energy, demonstrating a mode at 414 cm-1. This correlates well with Figure ____. By constructing a horizontal line at frequency = 414 cm-1, it can be observed that the branches intersect this line frequently. This implies that a significant proportion of k-points have vibrational modes of frequency 414 cm-1.

Loss of information in dispersion curve: how is that so?
Orthogonal to the dispersion curve: gives the number of energy states.
Information provided for the dispersion curve vs the density of states.
=== Computing the Free Energy Using the Harmonic Approximation ===
Figure ___ demonstrates the relationship between the PM shrinking factor used and the computed Helmholtz free energy of the system.
[FIGURE HERE]

==== Determining Optimal Grid Size for MgO ====
The suitability of the grid size used can be evaluated by investigating the variation in total energy of the unit cell when the number of k-points is increased. Figure ____ demonstrates the results obtained.

[GRAPH OF ENERGY AGAINST K-POINTS]

As the number of k-points increased, the number of peaks on the density of states plot increases accordingly. An 8x8x8 grid demonstrates most of the features present in the 16, 20 and 64 plots, and could possibly be the minimum reasonable approximation for the density of states.

Nonetheless, for a more conclusive appraisal, the convergence of energy values is a useful indication.

From Figure ____, the free energy of MgO is observed to increase and converge to a value of -40.926 483 eV, and it is observed that this occurs for a grid size of 8x8x8. (Why does it do so? Explain here)

A 2x2x2 grid is sufficient for calculating the free energy of MgO to 1 meV. A 4x4x4 grid is necessary for a precision to 0.5 meV and 0.1 meV.

=== Thermal Expansion ===
The Helmholtz free energy of a crystal is given by the sum of the energies of independent vibrational waves. The energy level <math display="inline">n</math> of a quantum harmonic oscillator are given by:

<math display="block">E_n = \left( n+ \frac{1}{2} \right) h \nu </math>

where <math display="inline">h</math> is Planck's constant and <math display="inline">\nu</math> is the frequency of energy level <math display="inline">n</math>. For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators, the vibrational energy is given by:

<math display="block">E_{vib} = \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu</math>

For a canonical <math display="inline">NVT</math> ensemble, the partition function is given by

<math display="block">Z = \sum_n e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

where <math display="inline">\beta = \frac{1}{kT}</math> and <math display="inline">E_n</math> enumerates all vibrational energy states.

For <math display="inline">N</math> atoms and <math display="inline">3N</math> independent harmonic oscillators,

<math display="block">Z_N = \prod_n^{3N} e^{\left( - \beta \varepsilon_n \right)} = \frac{e^{- \frac{\beta \varepsilon_n}{2}} }{1-e^{- \beta \varepsilon_n}}</math>

The phonon entropy can then be expressed in terms of the partition function:
<math display="block">S=-k_B \ln Z_N</math> where <math display="inline">k_B</math> is the Boltzmann constant.

Given the relation
<math display="block">F=U+TS</math>
where <math display="inline">U</math> is the internal energy of the system— for a crystal this is its electric potential energy
<math display="block">U_E = \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}}</math>

where <math display="inline">i</math> and <math display="inline">j</math> are the indices of the ions, <math display="inline">r_{ij}</math> is the distance between <math display="inline">i</math> and <math display="inline">j</math> and <math display="inline">\varepsilon_0 = 8.8542 \times 10^{-12} F\cdot m^{-1} </math>

The Helmholtz free energy of a crystal is thus given by

<math display="block">F= \frac{1}{2} \sum_{i=1}^{3N} q_i \sum_{j=1, j\ne i}^{3N} \frac{q_j}{4\pi \varepsilon_0r_{ij}} + \sum_{n=1}^{3N} \left( n+ \frac{1}{2} \right) h \nu + k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math>

This equation could be used to qualitatively rationalise the free energy dependence on temperature. The data obtained is plotted in Figure ____.
[FIGURE]
Particularly, there are two salient regimes of interest. At low temperatures, T < 100 K, the graph is flat. However, at high temperatures, the behaviour is approximately linear. These observations are in agreement with the above equation, which highlights the temperature dependence of entropy <math display="inline">S</math>. At low temperatures, the term <math> k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> is extremely small, and hence the free energy term is dominated by the internal energy of the crystal. At high temperatures, the term <math>-k_BT\sum_n^{3N}\ln \left( 1-e^{-\beta \varepsilon_n} \right)</math> dominates and therefore the free energy of the system appears to have a dependence in temperature.

==== Variation of Lattice Parameter with Temperature ====
[[File:Syl815CellPvsT.png|thumb|center|600px|This figure illustrates the variation in cell parameter of MgO with temperature]]
As the temperature increases, the lattice parameter increases. It can thus be observed that the cell volume <math display="inline">V </math> has a dependence on temperature <math display="inline">T </math>, and the thermal expansion coefficient <math display="inline">\alpha </math> is given by

<math display="block">\alpha = \frac{1}{V} \left( \frac{\partial V}{\partial T}\right)_P </math>

At the low temperature limit, <math display="inline">T\le 80</math> °C, <math display="inline">\alpha = 1.2215 \times 10^{-5}</math>.

Compute the coefficient of thermal expansion for MgO
How does this compare to that measured? Find a measurement in the literature or on the web - at what temperature was the measurement made?
What are the main approximations in your calculation?

=== MD Simulation ===
The cell volume per formula unit of MgO was plotted against temperatures between 20 K to 4000 K.
[[File:Syl815MD.png|thumb|center|600px|"Experimental Data MD"]]
Under MD, the cell volume generally increases linearly with temperature throughout. By considering the mean kinetic energy of the crystal
<math display="block">\left \langle E_k \right \rangle = \frac{1}{2} M \sum_{i=1}^N v_i^2 = \frac{3}{2} Nk_BT_{MD} </math>
Where <math display="inline">\left \langle E_k \right \rangle </math> is the average kinetic energy of the atoms, <math display="inline">M</math> is the mass of the crystal lattice, and <math display="inline">v_i</math> represents the velocity of the atom <math display="inline">i</math>. It can be observed that the cell energy is linearly dependent on temperature. In a constant pressure system, this would result in volume expansion as temperature increases.

It can be observed that at high temperatures when <math display="inline">T\ge 2000 K </math>, more noise is present in the data due to the large cell volume and the large kinetic energy of the atoms.

[[File:Syl_MDvsQHA.png|thumb|center|600px|"This figure compares the data obtained for the thermal expansion of MgO under QHA and under MD."]]

At extremely low temperatures of <math display="inline">T\le 200 K </math>, QHA predicts a larger cell volume than MD. This can be attributed to the significant quantum effects at such low temperatures. Since MD only accounts for the kinetic energy of the atoms and neglects zero point vibrations, it predicts a smaller cell volume with the atoms closer together.

The data obtained for MD and QHA demonstrate strong agreement for temperatures between 200 to 1000 K. At these temperatures, the thermal energy of the system is sufficiently large such that the motion of the particles can be described classically.

Links:
https://pubs.acs.org/doi/pdf/10.1021/ct0500904
https://journals.aps.org/prb/pdf/10.1103/PhysRevB.43.5024

Conclusion

Between MD and Thermal expansion model
Compare quantitatively with values
Why you think it shouldn’t be linear

== References ==

File:Syl MDvsQHA.png

2018-02-14T11:42:18Z

Syl815: