Rep:Mod:TAB9192

2017-01-26T23:05:33Z

Aob114: /* The Bulk Modulus of MgO */

File:Alwin Cmprs.png

2017-01-26T22:52:36Z

Aob114: Aob114 uploaded a new version of File:Alwin Cmprs.png

File:Alwin Cmprs.png

2017-01-26T22:51:28Z

Aob114: Aob114 uploaded a new version of File:Alwin Cmprs.png

File:Alwin Cmprs.png

2017-01-26T22:49:26Z

Aob114:

2017-01-26T21:31:57Z

Aob114: Aob114 uploaded a new version of File:Alwin ThExpMelt.png

File:Alwin ThExpMelt.png

2017-01-26T21:23:10Z

Aob114:

Rep:Mod:TAB9192

2017-01-26T17:36:42Z

Aob114: /* Results and Discussion */

= A Computational Study of the Thermodynamic Properties of MgO =

== Abstract ==
Thermal Expansion of MgO, MD and QH, expansion coefficient. Compressibility of MgO. Lattice Energy of MgO. Phonon DOS dependence on shrinking factor.

The volumetric coefficient of thermal expansion (<math> \alpha </math>) of MgO was computed through the use of both molecular dynamics (MD) simulations and phonon calculations under the quasi-harmonic approximation (QHA). Good agreement between both methods and existing literature was found for temperatures within a range of ???-???K. At high temperatures, the MD simulations showed anomalous behaviour in <math> \alpha </math> at temperatures between <math> 3000 - 3050 \, \mathrm{K} </math>, matching the experimental melting range at <math> 0 \, \mathrm{Pa} </math> closely. The volume at <math> 300 \mathrm{K} </math> was determined as a function of pressure, and used to calculate the bulk modulus at <math> 0 \, Pa </math> (<math> \mathrm{K}_{0} </math>) of ??? GPa with a pressure derivative <math> \mathrm{K}'_{0} </math> of ???. Again, both MD simulations and QHA-based calculations were performed

Zerr and Boehler (1994) for melting point
== Introduction ==
Structure, Importance in geochemistry, Advantages of theoretical over experimental characterisation, Comparison with literature.

Solid MgO adopts the well-known rock-salt structure, which consists of two superimposed but offset face-centred cubic lattices. A detailed theoretical understanding of its thermodynamic properties is paramount to a variety of applications. For example, MgO is a commonly used refractory material, in part due to its chemical stability at high temperatures. Also, the high-pressure melting of MgO in the earth's lower mantle plays a major role in the seismology of so-called ultra-low velocity zones(http://dx.doi.org/10.1016/j.physleta.2005.12.065). However, experimental investigations at such conditions rely on either shockwaves or laser-heated diamond anvil cells, both of which require highly sophisticated equipment and are thus very expensive. A convenient alternative is the use of computational methods. The present work therefore characterises the compression and thermal expansion of MgO through the use of both MD and QHA-based methods, comparing their accuracy and range of validity.

== Methodology ==
The software used to perform all calculations and simulations was the General Utility Lattice Program (GULP) (GULP: A computer program for the symmetry-adapted simulation of solids). The internal energy (excluding the zero-point phonon energies) of the perfect crystal at <math> 0 \, \mathrm{K} </math> was calculated
=== Lattice Energy of MgO ===
Both the MD simulations and the phonon calculations under the QHA require a pair potential that determines the inter-particle interactions. Throughout this study, the Coulomb-Buckingham potential was used, taking the general form below:

<math> U(r) = A\exp(-Br) - \frac{C}{r^6} + \frac{q_1 q_2}{4\pi \epsilon_{0}r} </math>

Where the constants <math> A </math>, <math> B </math> and <math> C </math> are chosen to match experimental data, and are shown in table 1.

{| class="wikitable" border="1"
|+ '''Table 1. ''' Coulomb-Buckingham parameters used.
|style="width: 110px;"|'''Pair'''
|style="width: 110px;"|'''A'''
|style="width: 110px;"|'''B'''
|style="width: 110px;"|'''C'''
|style="width: 110px;"|'''Cutoff / <math> \AA </math>'''
|-
|O-O
| <math> 0.130 \times 10^4 </math>
| <math> 0.300 </math>
| <math> 0.00 </math>
| <math> 12 </math>
|-
|Mg-O
| <math> 0.228 \times 10^4 </math>
| <math> 0.149 </math>
| <math> 27.9 </math>
| <math> 12 </math>
|-
|}
The cutoff indicates the value of <math> r </math> above which the short-range Pauli-repulsion and dispersive potential interactions are discarded (the first two terms of the Coulomb-Buckingham potential). Note that the Coulombic interactions were summed using the Ewald method, where interactions are calculated in Fourier-space, due to significantly faster convergence.

=== Phonon Calculations under the QHA ===
Phonons are the collective vibrational modes of a solid. For small distortions from some equilibrium interatomic spacing <math> r_{0} </math>, the potential can be approximated by a parabola, as can be easily verified by Taylor-expansion around <math> r_{0} </math>:

<math> V(r) \approx V(r_{0}) + \Bigg[ \frac{\mathrm{d}V}{\mathrm{d}r} \Bigg]_{r=r_{0}}(r-r_{0}) + \frac{1}{2}\Bigg[ \frac{\mathrm{d}^2 V}{\mathrm{d}r^2} \Bigg]_{r=r_{0}}(r-r_{0})^2 </math>

Since <math> r_{0} </math> is a potential energy minimum, <math> \Bigg[ \frac{\mathrm{d}V}{\mathrm{d}r} \Bigg]_{r=r_{0}} = 0 </math>, and so, defining <math> V(r_{0})=0</math> yields <math> V(r) \approx \frac{1}{2}\Bigg[ \frac{\mathrm{d}^2 V}{\mathrm{d}r^2} \Bigg]_{r=r_{0}}(r-r_{0})^2 = \frac{1}{2}k(r-r_{0})^2 </math>.

This is the origin of the 'harmonic' part of the QHA. However, this is clearly no a widely applicable model, as at any moment in time, the bond distance is equally likely to be shorter or longer than <math> r_{0} </math> for all excited states, and so no thermal expansion could take place.

Under the QHA, the Helmoltz free energy per unit cell consists of contributions from the lattice energy (described in the previous section), the zero-point energy of all phonons in the Brillouin zone, and the energy of occupied excited states:
<math> F(V,T) = U + \int \mathrm{d}\omega g(\omega ) \Bigg[ \frac{\hbar\omega}{2} + k_{\mathrm{B}}T\log (1- e^{-\frac{\hbar \omega}{k_{\mathrm{B}}T}}) \Bigg] </math>

The origin of the 'quasi'-part of the QHA is that the volume of the crystal is allowed to vary with a change in conditions in such a way as to minimise the Helmholtz free energy. In general, the force constants of the phonons decrease with increasing volume, and so the vibrational frequencies <math> \omega </math> are lower (a relationship known as Badger's rule). On one hand, this causes the energy of occupied states (the final term) to become more negative and lowers the zero-point energies, but results in a less negative lattice energy for volumes larger than the equilibrium volume at <math> 0 \, \mathrm{K} </math>. There are therefore two opposing contributions to the Helmholtz free energy with a change in volume, and so the volume can be optimised in such a way as to minimise it, yielding the equilibrium structure at the given temperature and pressure.

Clearly, the only volume-dependent quantity in the expression for the Helmoltz free energy can be the phonon DOS, <math> g(\omega) </math>, which therefore contains all of the information required to calculate the Helmholtz free energy at a given temperature.
At each volume, this expression is the analytical result derived from the partition function of a harmonic oscillator, summed over all possible modes. Computationally, for any real system this integral must be approximated by a summation over a finite number of frequencies, as the phonon DOS cannot generally be determined analytically. The phonon DOS is itself summed over a selection of points in the Brillouin zone, which are determined by the shrinking factor. For a 3D lattice, this is a set of three numbers defining the 'resolution' of k-points in the Brilluoin zone at which the phonon DOS is calculated. However, a larger number of sampling points comes with greater computational expense, and so, prior to any studies of the thermodynamic properties of a real system, it is important to determine the convergence of the phonon DOS with the shrinking factor, such that a reasonable compromise between accuracy and speed can be reached. This can be assessed both visually from plots of the phonon DOS at different shrinking factors, or by monitoring associated properties such as the Helmoltz free energy.

=== Molecular Dynamics Simulations ===
Molecular dynamics (MD) simulations are a purely classical approach to modelling the behaviour of systems. The basic principle is that the time evolution of a system obeying Newtonian mechanics is simulated through the use of some iterative algorithm. In this case, the so-called leapfrog algorithm was used to integrate the equations of motion.

=== The Coefficient of Thermal Expansion ===
The thermal expansion of a material can be characterised by its coefficient of thermal expansion:
<math> \alpha = \frac{1}{V}\Bigg(\frac{\partial V}{\partial T} \Bigg)_{P} </math>

The division by the volume of the material is required for <math> \alpha </math> to be an intensive quantity. This generally varies with temperatue, and is thus conventionally computed by fitting the volume-temperature data with a polynomial, and evaluating its derivative in the requuired range (note hat this renders extrapolation meaningless).

=== The Bulk Modulus ===
The isothermal bulk modulus (the inverse of the compressibility) of a material can be obtained by fitting the Birch-Murnghan equation-of-state (EOS) (http://dx.doi.org/10.1016/S0038-1098(00)00412-9) to the volume-pressure data at constant temperature:

<math> P(V) = \frac{3}{2}K_{0}\Bigg[ \bigg(\frac{V_{0}}{V} \bigg)^{\frac{7}{3}} - \bigg(\frac{V_{0}}{V} \bigg)^{\frac{5}{3}} \Bigg] \Bigg[ 1 - \frac{3}{4}(4-K'_{0})\bigg[ \bigg(\frac{V_{0}}{V}\bigg)^{\frac{2}{3}} -1\bigg]\Bigg] </math>

Where <math> V_{0} </math> is the volume at <math> P = 0 \, \mathrm{Pa} </math>, <math> B_{0} = -V\Bigg(\frac{\partial P}{\partial V}\Bigg)</math> is the isothermal bulk modulus in the absence of an external pressure and <math> B'_{0} = \Bigg( \frac{\partial B}{\partial P}\Bigg)_{P=0} </math> its pressure derivative.

== Results and Discussion ==
=== The Total Lattice Energy ===
The total lattice energy per primitive unit cell was determined to be -41.07531759 eV.

=== The Phonon Dispersion ===
The phonon dispersion was calculated at 50 points along the W-L-G-X-W-K path in the Brillouin zone, the graph of which is shown in figure 1. Figure 2 contains an analogous (albeit the path in the first Brilluoin zone somewhat differs) graph derived from neutron scattering data, obtained by Sangster et al. (http://dx.doi.org/10.1088/0022-3719/3/5/017). As can be expected from the general rule that the total number of branches is given by <math> d\times m </math>, where <math> d </math> is the dimensionality of the system and <math> m </math> the number of atoms in the unit cell, there are six branches, since the primitive unit cell contains one Mg and one O ion. In addition, <math> d </math> of these branches should be acoustic, in other words, have a minimum at <math> k = \Gamma = 0 </math>. Clearly, this is also the case.

{|
|-
| [[File:Alwin_MgOPhDisp.png|500x500px|Figure 1: The phonon dispersion calculated at 50 points along the W-L-G-X-W-K path. Three acoustic and three optical branches are clearly visible.|none|thumb]]
| [[File:Alwin_MgOPhDispRef.png|500x500px|Figure 2: The phonon dispersion as obtained experimentally by Sangster et al. Clearly, the qualitative match is excellent, as the ordering of the bands is identical at aall symmetry points present in both figures.|none|thumb]]
|}

=== The Phonon DOS ===
Figures 3-8 illustrate the convergence of the phonon DOS with increasing mesh resolution (shrinking factor). As can be seen, the change in visual appearance of the DOS-curve as a result of a change in the mesh resolution decreases with increasing resolution, to the point where there is no significant visual difference between the DOS-curves at resolutions of 30x30x30 and 40x40x40, and so the convergence has satisfied at least the visual criterion.
{|
|-
| [[File:Alwin_MgOPhDOS1by1.png|thumb|530px|Figure 3: Phonon DOS with a mesh resolution of 1x1x1. Clearly, this corresponds to the symmetry point L in the phonon dispersion graph.]]
| [[File:Alwin_MgOPhDOS5by5.png|thumb|530px|Figure 4: Phonon DOS with a mesh resolution of 2x2x2.]]
|-
| [[File:Alwin_MgOPhDOS10by10.png|thumb|530px|Figure 5: Phonon DOS with a mesh resolution of 10x10x10.]]
| [[File:Alwin_MgOPhDOS20by20.png|thumb|530px|Figure 6: Phonon DOS with a mesh resolution of 20x20x20.]]
|-
| [[File:Alwin_MgOPhDOS30by30.png|thumb|530px|Figure 7: Phonon DOS with a mesh resolution of 30x30x30.]]
| [[File:Alwin_MgOPhDOS40by40.png|thumb|530px|Figure 8: Phonon DOS with a mesh resolution of 40x40x40.]]
|}

Comparison of the phonon DOS for evaluated at only a single point in the Brillouin zone (figure 3) with the phonon dispersion curve (figure 1) shows that this point must correspond to the symmetry point L, as this matches the observed frequencies.

In this context, the phonon DOS is defined as the number of states per unit increase in energy at a given energy i.e. <math> DOS(E_{0}) = \frac{\mathrm{d} n}{\mathrm{d} E} </math>. Using the chain rule, this can be rewritten as <math> DOS(E_{0}) = \int_{E(\vec{k'})=E_{0}}\mathrm{d}\vec{k'} \frac{\mathrm{d} n}{\mathrm{d} |\vec{k'}|} \frac{\mathrm{d}S(\vec{k'})}{|\nabla_{\vec{k'}}E(\vec{k'})|} </math>, where the integral along the surface of <math> E(\vec{k}) = E_{0} </math> arises from the need to consider the contributions of all k-values for which <math> E(\vec{k}) = E_{0} </math>. This expression can of course generally not be computed analytically, but the important thing to note is that the DOS varies inversely with <math> \nabla_{\vec{k}}E(\vec{k}) </math>, i.e. the slope of the phonon dispersion curve. In this specific case, the frequency <math> \omega </math> is used in place of the energy, however, the same relationship holds true, as the zero-point energy of a particular phonon mode is given by <math> E = \frac{\hbar \omega }{2} \implies E \propto \omega </math>

=== Free Energy in the Quasi-Harmonic Approximation ===

Following the above discussion of the visual convergence of the phonon DOS curve, it should be noted that for the purpose of calculating thermodynamic quantities, it is far more important to know the level of accuracy that can be expected at a given mesh resolution. Thus, the Helmholtz free energy per unit cell of the MgO lattice at <math>300 \, \mathrm{K} </math> was calculated at various mesh resolutions. The results are summarised in table 2 and illustrated in figures 9 and 10.

{| class="wikitable" border="1"
|+ '''Table 2. ''' Helmholtz free energy dependence on mesh resolution.
|style="width: 110px;"|'''Mesh resolution'''
|style="width: 110px;"|'''F/eV '''
|style="width: 110px;"|'''<math> \Delta F </math>/eV '''
|-
|1x1x1
| <math>-40.930301</math>
| N/A
|-
|2x2x2
| <math> -40.926609 </math>
| <math> 0.003692 </math>
|-
|3x3x3
| <math> -40.926432 </math>
| <math> -0.000177 </math>
|-
|4x4x4
| <math> -40.926450 </math>
| <math> -0.000018 </math>
|-
|5x5x5
| <math> -40.926463 </math>
| <math> -0.000013 </math>
|-
|6x6x6
| <math> -40.926471 </math>
| <math> -0.000008 </math>
|-
|7x7x7
| <math> -40.926475 </math>
| <math> -0.000004 </math>
|-
|8x8x8
| <math> -40.926478 </math>
| <math> -0.000003 </math>
|-
|9x9x9
| <math> -40.926479 </math>
| <math> -0.000001 </math>
|-
|10x10x10
| <math> -40.926480 </math>
| <math> -0.000001 </math>
|-
|15x15x15
| <math> -40.926482 </math>
| <math> -0.000002 </math>
|-
|20x20x20
| <math> -40.926483 </math>
| <math> -0.000001 </math>
|-
|25x25x25
| <math> -40.926483 </math>
| <math> 0.000000 </math>
|-
|30x30x30
| <math> -40.926483 </math>
| <math> 0.000000 </math>
|-
|40x40x40
| <math> -40.926483 </math>
| <math> 0.000000 </math>
|-
|}

{|
|-
| [[File:Alwin_FreeEnergy.png|thumb|530px|Figure 9: Variation of the Helmholtz free energy with mesh resolution (where a single number Y indicates a YxYxY grid).]]
| [[File:Alwin_FreeEnergyError.png|thumb|530px|Figure 10: Variation of the fractional change of the Helmholtz free energy, <math> |\frac{\Delta F}{F}| </math> , with mesh resolution.]]
|-
|}
The most distinctive feature of figures 9 and 10 is the rapid convergence of the fractional change in Helmholtz free energy to zero, and that of the Helmholtz free energy itself to a constant value of -40,926483 eV . To investigate the accuracy at each mesh resoltion it is necessary to compare not this change in Helmholtz free energy, which depends also on the choice of shrinking factors at which the free energy was evaluated, but instead the difference between the final converged free energy (taken to be -40,926483 eV) and the free energy at a particular mesh resolution. A plot of this is shown in figure 11.
{|
|-
| [[File:Alwin_FreeEnergyAcc.png|thumb|530px|Figure 11: Variation of the difference between the converged Helmholtz free energy and that for a given mesh resolution with mesh resolution.]]
|-
|}

Clearly then, an accuracy of 1meV and 0.5 meV can be achieved with a 2x2x2 grid, and one of 0.1 meV with a 3x3x3 grid. Whilst the discrepancy between the rapid convergence of the Helmholtz free energy, and the much more obvious changes in the visual appearance of the phonon DOS curve might seem counterintuitive, it is important to note that the energetic contribution (both zero-point and excited-state) from phonons is, particularly at 300 K, very small compared to the lattice energy. Therefore, the Helmholtz free energy can be determined with reasonable accuracy even when the phonon DOS curve does not strongly resemble its converged form.

However, as the calculations were not particularly computationally expensive, a mesh resolution of 40x40x40 was used throughout this investigation, to ensure that a lack of sample points in the Brillouin zone was not a significant contributor to the uncertainty in the final results. This grid size would certainly be appropriate for similar compounds such as CaO. Since the dimension of the first Brillouin zone (BZ) is proportional to <math> \frac{1}{a} </math>, where <math> a </math> is the characteristic length of the Wigner-Seitz cell, a material with a considerably larger Wigner-Seitz cell than MgO, such as a Zeolite (e.g. Faujasite), will have a smaller first BZ and thus require fewer sample points over which to evaluate the phonon DOS at similar levels of accuracy, i.e. a smaller mesh resolution. For a metal such as Lithium, another factor comes into play. Namely, the delocalised electrons are able to screen the repulsions between approaching positively charged nuclei, resulting in decreased phonon dispersion ('width' of the branches). As a result, the energies of the phonon modes at various k-points are more similar than in an insulator, and so an average over a fewer number of sample points describes the phonon DOS curve with equal accuracy. However, the Wigner-Seitz cell of Lithium is smaller than that of MgO, and so, conversely to the case of a Zeolite, this would tend to increase the mesh resolution required for the same level of accuracy. The question of the balance between these two competing factors is outside the scope of this investigation.

=== Thermal Expansion ===
Figures 13 and 14 show the behaviour of the lattice parameter (the edge-length of the primitive unit cell) and the Helmholtz free energy with increasing temperature, as calculated under the QHA.
{|
|-
| [[File:Alwin_ThExpFreeEn.png|thumb|530px|Figure 13: The Helmholtz free energy as a function of temperature]]
| [[File:Alwin_ThExpLatPar.png|thumb|530px|Figure 14: Dependence of the lattice parameter on temperature]]
|-
|}

== Conclusions ==

File:ThExpLatPar.png

2017-01-26T17:25:10Z

Aob114:

File:Alwin ThExpFreeEn.png

2017-01-26T16:54:24Z

Aob114:

Rep:Mod:TAB9192

2017-01-26T16:10:41Z