Rep:Mod:SG94

2016-01-29T00:15:18Z

Sg2613: /* Molecular Dynamics */

2016-01-28T14:59:46Z

Sg2613: /* Computing the Free Energy */

== Introduction ==
MgO - Aims/crystal structures - unit cells

DLV

GULP

== Internal Energy of MgO crystal ==
[[File:SG Primitive MgO.png|right|thumb|250px|Fig1. Primitive cell of MgO]]
GULP was used to define the MgO crystal structure using a primitive unit cell (Figure 1). The cell vectors of this cell are (‎Å) and the LogFile can be found [https://wiki.ch.ic.ac.uk/wiki/images/6/6c/SG_PMgO.txt here]:

{|
|0.00000||2.10597||2.10597
|-
|2.10597||0.00000||2.10597
|-
|2.10597||2.10597||0.00000
|}
The Log File shows that the primitive cell has sides of length 2.9783 ‎Å and internal angles of 60 degrees. The total lattice energy is calculated by taking short-range repulsions and electrostatic potentials into account and summing over the lattice and represents the energy required to separate the ions of the lattice to infinite separation. The energy of the MgO lattice is -41.07531759 eV and corresponds to the energy of the lattice at 0 K (with no anharmonic contributions).

== Phonon Modes of MgO ==
The lattice vibration modes (also known as phonon modes) can be represented by a dispersion curve or density of states (DOS).

=== Phonon Dispersion Curves ===
[[File:SG_Dispersion.png|right|thumb|300px|Fig2. Phonon Dispersion of MgO along 50 points]]

In an infinite lattice, there are an infinite number of vibrations which can be simplified by considering the vibrations in the primitive unit cell of the crystal. Every possible vibration of the crystal has a '''k'''-label (symmetry points: <math>W, L, G, W, X, K</math>) which is related to its direction and wavelength. By summing over '''k''', the total number of vibrations can be calculated. This is achieved by using GULP to calculate the phonon dispersion and by computing the phonons at 50 points along the path W-L-G-W-X-K in '''k'''-space.

The phonon dispersion curve (Figure 2) shows how the frequency of vibration changes with a change in k and the computed curve agrees with that found in literature (DOI:10.1016/j.cocom.2014.10.005).

As there is more than one atom in the primitive unit cell of MgO (N = 8*1/8 + 1 = 2), two types of phonons are displayed in the dispersion relation: acoustic and optical, corresponding to in-phase and out-of-phase movements respectively. As the calculations are taking three directions into account (x, y and z), there are a total of six phonon modes. This is explained by G. E. Peckham, where the three acoustic modes consist of one longitudinal and two transverse modes (which pass through the origin and represents the crystal behaving as a continuum). The three optical modes (3N-3) have non-zero frequencies and zero gradient at k=0. They are optical modes as they can interact with IR radiation.

By summing all of the possible k-points of the lattice, the free energy of the lattice can be computed.

=== Phonon Density of States (DOS) ===
The<nowiki> </nowiki>density of states summarises the dispersion curves by calculating the average number of vibrational modes at each frequency. Shrinking factors<nowiki> </nowiki>define the grid of '''k'''-points that will be used to compute the DOS curve and their effects on the curve were investigated.

A 1x1x1 grid computes the DOS from a single point and is shown in Table 1. It can be<nowiki> </nowiki>seen that there are four distinct peaks which corresponds to the '''k'''-point<nowiki> </nowiki>L. There are only four peaks present in the DOS rather than 6 as bands become degenerate at frequencies at 288.49 and 351.76 cm-1. This explains why the intensities of these two peaks are double that of the other two at higher frequencies (676.23 and 818.82 cm-1).

{|align="center"
|[[File:SG_DOS(1).png|thumb|250px|Fig3. Density of States at grid size 1x1x1]]
|[[File:SG_DOS(14).png|thumb|250px|Fig4. Density of States at grid size 14x14x14]]
|[[File:SG_DOS(15).png|thumb|250px|Fig5. Density of States at grid size 15x15x15]]
|[[File:SG_DOS(30).png|thumb|250px|Fig6. Density of States at grid size 30x30x30]]
|}

By increasing the grid size, more points in '''k'''-space<nowiki> </nowiki>are taken to compute a more representative DOS and the graph starts to form a curve at 16x16x16. Several grid sizes were investigated to find a<nowiki> </nowiki>compromise between accuracy of results and CPU time (below).<nowiki> </nowiki>The minimum grid size for a good approximation is the 15x15x15 grid. Below this size, there are cuts on the x-axis (not a curve) as the number of sampled '''k'''-points were insufficient; and above this grid<nowiki> </nowiki>size, the DOS becomes smoother and more representative with an increase<nowiki> </nowiki>in the number of sampled possible vibrations. It can be seen that there<nowiki> </nowiki>is very little difference in the DOS of the two largest grid sizes. This occurs when a sufficient number of points in '''k'''-space are taken at a given size; and a further increase in would make little/no difference. The CPU time for the 40x40x40 grid size was noticeably longer than the 30x30x30 with very little difference in the DOS and so it was decided that the latter was the optimum size used for future calculations.

The phonon dispersion curve states the possible frequencies of phonons at each '''k'''-label and by differentiating this curve, the DOS produced gives information about the entire set of phonon frequencies. The frequencies which more '''k'''-labels correspond to increases the density of states at that frequency.

== Computing the Free Energy ==
{| class="wikitable" style="float:right;margin-left: 10px;"
|+ style="caption-side:bottom;"|''Table 2. The Effects of Grid Size on the Helmholtz Free Energy of MgO''
|-
! scope="col"| Grid size
! scope="col"| Helmholtz Free Energy / eV
! scope="col"| Accuracy
|-
! scope="row"| 1x1x1
| align="center"| - 40.9303|| align="center"| 0.1 eV
|-
!scope="row"| 2x2x2
| align="center"| - 40.9266|| align="center"| 0.1 eV
|-
! scope="row"| 3x3x3
| align="center"| - 40.9264||align="center"| 1 meV
|-
! scope="row"| 4x4x4
| align="center"|- 40.9265||align="center"| 0.5 meV
|-
! scope="row"| 10x10x10
| align="center"|- 40.9265||align="center"| 0.1 meV
|-
! scope="row"| 20x20x20
| align="center"|- 40.9265||align="center"| 0.1 meV
|-
! scope="row"| 30x30x30
| align="center"|- 40.9265||align="center"| n/a
|}

The quasi-harmonic approximation was used to compute the free energy of MgO by summing all of the possible vibrational modes of the crystal at 300 K. The free energies computed at various sizes using GULP are displayed in Table 2. An increase in grid size lowers the free energy of the lattice by increasing the degree of accuracy. This is the case due to the increase in number of microstates sampled (<math>S=klnW</math>), increasing the entropic term and in turn, lowering the Gibbs free energy (<math> G=H-TS</math>).

Table 2 shows that 30x30x30 grid size is an appropriate choice for accurate results with an appropriate CPU time and found the free energy converge to - 40.9265 eV. These grid sizes are applicable to other similar crystal structures such as calcium oxide (CaO), where the lattice constant (4.800 Å: http://link.springer.com/chapter/10.1007%2F10681719_224#page-1) is similar to that of MgO (4.211 Å: http://link.springer.com/chapter/10.1007%2F10681719_206#page-1).

Running a calculation on structures with a significant difference in lattice constants would require a different sized grid for the same level of accuracy. For zeolite structures such as faujasite with a larger lattice constant (between 24.2-25.1 Å): http://www.rigaku.com/downloads/journal/Vol12.2.1995/kaduk.pdf), a smaller grid size (fewer sampled '''k'''-points) is required to achieve sufficiently accurate results.

==Thermal Expansion of MgO==
When a substance is heated, the energy is converted to kinetic energy in its molecules, resulting in the expansion of the material. The linear thermal expansion coefficient (CTE) indicates the degree of expansion in one dimension (i.e. the lattice constant) of a material upon heating (Lit: http://www.owlnet.rice.edu/~msci301/ThermalExpansion.pdf) and the equation used to calculate this value is given by:

:<math>\alpha_L=\frac{1}{L_0}\,\frac{dL}{dT}</math>

Where <math>L_0</math> is the lattice constant at 0 K and <math>dL/dT</math> is the rate of change of the lattice constant per unit change in temperature.

The volumetric CTE at constant pressure is given by:

:<math>\alpha_V = \frac{1}{V}\,\left(\frac{\partial V}{\partial T}\right)_p</math>

This equation is important for a gas as the volume will vary with a change in pressure but this is not important for MgO as it is a solid.

The crystal structure of MgO means that it has identical properties in all three directions, i.e. it is isotropic. This means that the volumetric CTE should be three times that of the linear CTE: <math>\alpha_V = 3\alpha_L</math>

The CTEs of MgO were predicted using the the quasi-harmonic approximation and molecular dynamics.

===Quasi-Harmonic Approximation===
{|style="float:right; margin-left: 10px;"
|[[File:SG_Gibbs_vs_T.png|thumb|400px|Graph 1. Effect of temperature on the Gibbs free energies of MgO]]
|-
|[[File:SG_Lattice_const_vs_T.png|thumb|400px|Graph 2. Effect of temperature on the lattice constants of MgO]]
|}

This method optimises the MgO structure by minimising the Gibbs free energy through computing the internal energy and phonons at a sequence of geometries. By varying the temperature from 0 K to 1000 K in 100 K intervals, the effects of temperature on the Gibbs free energy and lattice constants are discovered (Graphs 1 and 2).

The negative correlation between Gibbs Energy and temperature can be explained through the definition of Gibbs where temperature is negatively proportional to the free energy; while the positive correlation between the lattice constant (bond length) and temperature proves the thermal expansion of MgO upon an increase in temperature.

:<math>G=H-TS</math>

Where G is the Gibbs free energy, H is the enthalpy (J mol-1), S is the entropy (J K-1 mol-1) and T is the temperature in Kelvin (K).

The volumetric expansion coefficient from 0 - 1000 K was calculated to be 2.2485x10-5 K-1 while the linear expansion coefficient for the same temperature range is 7.4400x10-6 K-1. The volumetric CTE is approximately three times that of the linear CTE and so the isotropic property is confirmed. From 0 - 300 K, the thermal CTE is 9.47x10-6 K-1 which is in general agreement with the literature value of 10.5x10-6 K-1. (DOI: 10.1007/10681719_206. http://link.springer.com/chapter/10.1007/10681719_206)

This method explains the thermal expansion based on the harmonic approximation while including anharmonic contributions to an extent (http://www.sciencedirect.com/science/article/pii/S0010465503005472). This allows the harmonic approximation to hold for every value of the lattice constant (i.e. phonon frequencies are volume-dependent), and explains why the bond length increases with temperature (as volume increases). The quasi-harmonic model in reciprocal space reduces the size of the force constant matrix by approximating that the vibration of the atoms in the primitive unit cell has the same magnitude and direction and only differs in phase. (B. T. M. Willis and A. W. Pryor, Thermal Vibrations in Crystallography �Cambridge University Press, Cambridge, 1975. ). This simplifies the model and reduces the effort required to calculate the free energy and other thermal properties of the MgO crystal. (http://www.ces.clemson.edu/~gli/papers/jap06.pdf).

The linear CTE, however, suggests that the relationship between temperature and the lattice constant is a linear one, which is not reflected in Graph 2. This can be explained through an over-simplification of the QH model

===Molecular Dynamics===
[[File:SG_Cell_vol_vs_T.png|left|thumb|400px|Graph 3. Comparing the Quasi-Harmonic and Molecular Dynamics methods and their effects of temperature on cell volume]]

The molecular dynamics method studies the physical movements of atoms and molecules where atomic and molecular interactions are allowed. The system is studied and evolves over a period of time according to Newton's second law: <math>F=ma</math>. The timestep at which the vibrations to be sampled was taken to be 1 femtosecond (fs) and the number of equilbration steps and production steps were set to 500. This situation obeys the ergodic hypothesis where when the number of sampled configurations are sufficiently large, the average calculated is representative for the time averaged system. Originally due to L. Boltzmann. See part 2 of Vorlesungen über Gastheorie. Leipzig: J. A. Barth. 1898. OCLC 01712811. ('Ergoden' on p.89 in the 1923 reprint.) It was used to prove equipartition of energy in the kinetic theory of gases. Increasing the size of the crystal used in the MD calculations rather than a single primitive cell allows more flexibility when the system is vibrating. This allows the movement of atoms that are not in phase and so the structure used contains 32 units (64 atoms).

In both cases, the anharmonic effects become significant at high temperatures and Graph 3 shows that the cell volumes obtained from the QH model is in good agreement with the MD data. The values obtained from the QH model are slightly higher than those from the MD simulations. The strong linear fit between the cell volume and the temperature using this method suggests that the calculations are based on a less harmonic approximation.
http://www.ces.clemson.edu/~gli/papers/jap06.pdf

The thermal CTE calculated using the molecular dynamics method is 2.9993x10-5 K-1 between 100-1000 K, similar to the value computed from the QH model of 2.4984x10-5 K-1, and the value from 0 - 300 K (3.0280x10-5 K-1) agrees well with the literature value of 1.05x10-5 K-1.