=Introduction=
The experiment aims to investigate the coefficient of thermal expansion α of a crystal of MgO, using quasi-harmonic and molecular dynamics approximations under RedHat Linux environment. This intrinsic property, α is calculated as 2.6544 x10-5 K-1 and 3.2024 x10-5 K-1 between 100 K and 1000 K for each approximation. The interface is DLV, a package for the visualisation of materials, structures and properties. The calculations are completed by a classical simulation program GULP in DLVisualize.

Energies, vibrations and lattice constants of MgO crystal are computed and are used to plot graphs to decide the thermal expansion coefficient. Vibrations are important in studying thermal properties, phase transitions, electrical properties,etc. The graph of vibrational frequency against wave vector k, generated by the computer, is called Phonon dispersion. The wave vector k is a vector which points in the direction of propagation of the wave. It defines the vector in reciprocal space and is calculated as following:

<math>k(cm^{-1})=\frac{2\pi}{\lambda}</math>(1), where λ is wavelength.

The phonon dispersion diagram contains acoustical and optical branches, sharing a similar concept as band structures[1]. For the three-dimensional crystal lattice that has two atoms Mg and O in its primitive cell, the dispersion relations exhibit a total of six branches: two types of phonons in x, y, z directions. The boundaries at −π/''a'' and π/''a'' are called Brillouin zone.

Vibrational density of states is also related. it is an average over all k-points yielding the number of vibrational modes at each frequency (the density of modes). The grid size chosen in the calculation has a strong impact on the shape: generally small grids give discrete DOS and large grids give continuous DOS.

The quasi-harmonic approximation(QHA) is a more adapted hypothesis than the harmonic approximation. The system is equivalent to a collection of independent harmonic oscillators and treats vibrations as if they did not interact. It is not appropriate at high temperatures because of phonon-phonon interactions[2]. It has limitation that it does not predict bond dissociation. QHA is used to compute vibrational Helmholtz energy of MgO, with equation as following:

<math>F=E_0+\frac{1}{2}\sum_{k,j}\hslash\omega_{k,j}+k_BT\sum_{k,j}ln[1-exp(-\hslash\omega_{k,j}/k_BT)]</math> (2)

where E0 is the internal energy; the sum of ħω for each k, j is the zero-point energy; and the last term is the temperature-dependent free energy. Zero-point energy is the lowest possible energy that a quantum mechanical system may have.

The molecular dynamics approximation(MD) can capture the full anharmonicity at high temperatures. It simulates the vibrations as random motions of atoms inside a supercell. The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic evolution of the system. The system using classical Newton's Law:

<math>F=ma</math> (3)

In the calculation, energy is minimized as a function of the atomic position and of the lattice parameters.

= An Initial Calculation for MgO on its Structure and Internal Energy =
[[File:Primitive cell(em2015).PNG|left|200x200px|thumb|Fig.1 The structure of a primitive unit cell of MgO. ]]
[[File:Conventional cell(em2015).PNG|thumb|200x200px|Fig.2 conventional cell of MgO crystal]]
MgO crystal is a face-centered cube as shown on the right hand side. The conventional cell has eight atoms in it.

In the experiment, its crystal structure is described using a primitive unit cell, with cell parameters :a=b=c=2.9783Å ; α=β=γ=60.0000'''°'''. The primitive cell include one Magnesium ion (each Magnesium ion at the corner of the cell is shared by seven other cells and there are eight of this Magnesium ions) and one Oxygen ion(each Oxygen ion at the center of the cell is shared by one cell only and there is one of this Oxygen ions). The ratio of atoms in conventional cell and that in primitive cell is 4:1, which indicates that the ratio of volume between conventional cell and primitive cell is also 4:1.

The cell vectors of this cell are (in Angstrom):

0.000000 2.105970 2.105970

2.105970 0.000000 2.105970

2.105970 2.105970 0.000000.

The inter-atomic potential is 6.72119980 eV, and the total energy due to the inter-atomic forces should be -41.07531759 eV per primitive unit cell; it is the "binding energy" of the crystal - ie: the energy required to pull all of the ions of the crystal to infinite separation.

The '''Include Shells''' option should ''NOT'' be selected during calculation settings. This is because each ion is considered as a non-polarisable hard sphere here.

= Phonon Dispersion, DOS and Grid size =
Here the free energy is computed for increasing sizes of grid and its convergence to the infinite grid value monitored.

The special points along the conventional path in k-space are W-L-G-W-X-K. They are used because it is difficult to display the variation of the frequency over the full 3D '''k'''-space.[[File:Dispersion x 1x1x1(em2015).png|left|thumb|630x630px|Fig.3&4 The left hand side is the phonon dispersion diagram and the right hand side is the DOS diagram with 1x1x1 grid size.]]
{| class="wikitable"
!path in k-space
!W
!L
!G
!X
!K
|-
|kx
|1/2
|1/2
|0
|1/2
|3/8
|-
|ky
|1/4
|1/2
|0
|0
|3/8
|-
|kz
|3/4
|1/2
|0
|1/2
|3/4
|}
As mentioned in the introduction, the phonon dispersion diagram of MgO should have six branches. The vertical axis is frequency is vibration and the horizontal axis represents the path in k-space. Γ is the origin in the reciprocal space, just like 0 in the direct space.

The 1x1x1 grid size gives 4 very sharp peaks which indicates two similar higher density of states at around 280 cm-1 and 350 cm-1, and two similar lower density of states at around 670 cm-1 and 810 cm-1. The higher peaks and the lower peaks have a approximate ratio of 2:1 in height.

Furthermore, the DOS diagram can be correlated to the dispersion diagram. The point L in the dispersion diagram intersects six branches, each at 810 cm-1, 670 cm-1, 350 cm-1,350 cm-1,280 cm-1 and 280 cm-1. The two 350 cm-1 frequencies are degenerate at point L and so are the two 280 cm-1 frequencies. This, corresponds to the DOS diagram, is the 2:1 ratio in height.

The L point is 1/2 1/2 1/2.

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

Thus the K point is 1/2 1/2 1/2. If check with the log. file:

'''K point 1 = 0.500000 0.500000 0.500000. '''

This means that the speculation of the relationship between phonon dispersion diagram and DOS diagram is correct.

<nowiki> </nowiki>In order to use more '''k''' points, the the density of states calculation is reran, increasing the shrinking factors from 1x1x1, to 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 24x24x24, 36x36x36, 64x64x64. And the density of states diagrams are shown below.
{| class="wikitable"
![[File:DOS2x2x2(em2015).PNG|none|thumb|220x220px|Fig.5 Grid size: 2x2x2]]
![[File:DOS3x3x3(em2015).PNG|none|thumb|223x223px|Fig.6 Grid size: 3x3x3]]
![[File:DOS4x4x4(em2015).PNG|none|thumb|237x237px|Fig.7 Grid size: 4x4x4]]
![[File:DOS8x8x8(em2015).PNG|none|thumb|222x222px|Fig.8 Grid size: 8x8x8]]
|-
|[[File:DOS16x16x16(em2015).PNG|none|thumb|222x222px|'''Fig.9 Grid size: 16x16x16''']]
|[[File:DOS24x24x24(em2015).PNG|none|thumb|221x221px|'''Fig.10 Grid size: 24x24x24''']]
|[[File:DOS32x32x32(em2015).PNG|none|thumb|232x232px|'''Fig.11 Grid size: 36x36x36''']]
|[[File:DOS64x64x64(em2015).PNG|none|thumb|224x224px|'''Fig.12 Grid size: 64x64x64''']]
|}

When increasing the grid size, there are more peaks. When using a 8x8x8 grid size, the peaks connect to each other, giving continuous segments of lines. More and more of the possible vibrations are sampled and more features appear. When higher grid size is applied, the shape becomes more and more smooth and no apparent changes in the shape of the diagram after grid size 24x24x24. This means that the 24x24x24 grid size is the minimum/optimum.

The optimal grid size for MgO is appropriate for a calculation on a similar oxide like calcium oxide. The distance between atoms a is similar for Ca-O and Mg-O.

According to the equation of reciprocal space,

<math>a^{\ast} = \frac{2\Pi}{a}</math>(4), a* is similar as a is similar.

For Zeolite, which its primitive cell contains more atoms resulting in more branches, and larger lattice size (15.56 Å compared to 2.9783 Å of MgO) resulting in smaller reciprocal space. Thus a smaller grid size is enough to catch all the features from its phonon dispersion diagram. This is just opposite to a metal like lithium. It has only one atom in its primitive cell so only three branches in the phonon dispersion diagram. Lithium also has smaller lattice size, giving a larger a*. To put all factors together, lithium needs larger optimal grid size to include all features.

= Free Energy in the Harmonic Approximation =
The free energies of the different grid sizes in last section are reported here.
{| class="wikitable"
|+Table. 1.
!Grid Size
!Free Energy/eV
!Difference in energy/eV
!Accuracy
!Comments
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|<nowiki>-</nowiki>
|
| rowspan="9" |How does the free energy vary with grid size ?
As grid size increases, the free energy increases as well.

|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|<nowiki>-0.003692</nowiki>
|
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|0.000177
|accurate to 1 meV
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|<nowiki>-0.000018</nowiki>
|
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|<nowiki>-0.000028</nowiki>
|accurate to 0.5 meV
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|<nowiki>-0.000004</nowiki>
|accurate to 0.1 meV
|-
|24x24x24
|<nowiki>-40.926483</nowiki>
|<nowiki>-0.000001</nowiki>
|
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0.000000
|
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|0.000000
|}
The optimal grid size speculation is similar to that in the previous section.

= The Thermal Expansion Coefficient of MgO under Quasi-harmonic Approximation =
[[File:Energy vs t(em2015).PNG|thumb|592x592px|Fig.13 Graph of free energy against temperature between 0K and 1000K]]
[[File:Lattice vs t(em2015).PNG|thumb|591x591px|Fig.14 Graph of lattice constant against temperature between 0K and 1000K]]
Vary the temperature from 0K to 1000K in steps of 100K and compute the variation in the free energy and optimal lattice constant by optimizing the structure at each temperature. The numbers are extracted from the log. files and organized in the table 2.
{| class="wikitable"
|+Table. 2. The table of free energy and lattice constant at temperatures between 0K and 1000K with a variation of 100K
!Temperature/K
!Free Energy/eV
!Lattice Constant/Å
!Volume/Å3
|-
|0
|<nowiki>-40.90190628</nowiki>
|2.986563
|26.638824
|-
|100
|<nowiki>-40.90241965</nowiki>
|2.986563
|26.638824
|-
|200
|<nowiki>-40.90937738</nowiki>
|2.987605
|26.666716
|-
|300
|<nowiki>-40.92812471</nowiki>
|2.989390
|26.714542
|-
|400
|<nowiki>-40.95859417</nowiki>
|2.991629
|26.774613
|-
|500
|<nowiki>-40.99943595</nowiki>
|2.994135
|26.841955
|-
|600
|<nowiki>-41.04931542</nowiki>
|2.996820
|26.914231
|-
|700
|<nowiki>-41.10711924</nowiki>
|2.999643
|26.990362
|-
|800
|<nowiki>-41.17189187</nowiki>
|3.002588
|27.069936
|-
|900
|<nowiki>-41.24301814</nowiki>
|3.005635
|27.152431
|-
|1000
|<nowiki>-41.31984836</nowiki>
|3.008784
|27.237863
|}
The graph of free energy against temperature is shown in Fig.13, and the plot of lattice constant against temperature is shown in Fig.14. Fig.13 contains a downward curve and its gradient increases as temperature increases. The gradient becomes constant at higher temperatures, which indicated a straight line. Fig. 14 contains an upward curve and its gradient increases as temperature increases. The line becomes straight after 400K.
[[File:HightT(em2015).png|thumb|589x589px|Fig.15 Volume per formula unit at high temperatures]]
The coefficient of thermal expansion of a material is :

<math>\alpha=\frac{1}{Vo}(\frac{dV}{dT})P</math>= 2.6544 x10-5 K-1. (5)

The coefficient of volume thermal expansion is 4.0 × 10−5 K−1[2] in the temperature range from 300 K to 1250 K from A. S. Madhusudhan Rao and K. Narender's paper, which is larger than the calculated value.

For the harmonic approximation, the bond between atoms will eventually break if temperature is high enough, which indicates that the equilibrium distance between atoms is independent of temperature.The quasi-harmonic approximation expands upon the harmonic phonon model of lattice dynamics. In the quasi-harmonic approximation, bonds will never break but lengthen to infinite. In the diagrams below, the center of bonds are indicated in red in each approximations.
{| class="wikitable"
|[[File:Mp-5598-qha(em2015).png|thumb|200x200px|Fig.16 The Quasi-harmonic approximation[3] graph where the center of the bond changes. ]]
|[[File:800px-Morse-potential(em2015).png|thumb|217x217px|Fig.17 The harmonic approximation[4] graph where the center of the bond remains. ]]
|}
When a material is heated, the kinetic energy of that material is increased and it's atoms and molecules move about more. This means that each atom will take up more space due to it's movement so the material will expand.

Fig.15 is the graph of volume per formula unit under QHA at high temperatures. The blue points are the ones used to calculate the thermal expansion coefficient. The It is clear in the graph that the line deviates from the best fit line at high temperatures. Since the QHA assumes no bond-breaking between atoms even at high temperatures, the line should be linear. Thus the approximation is invalid at high temperatures.

= The Thermal Expansion Coefficient of MgO under Molecular Dynamics Approximation =
[[File:QH-MD(em2015).png|thumb|443x443px|Fig.18 Graph of volume per formula unit under Quasi-harmonic and Molecular Dynamics Approximation against temperature|left]][[File:QHA VS MD(em2015).png|thumb|497x497px|centre|Fig.19 Graph of Volume per formula unit against temperature of QHA and MD at high temperatures]]Fig.18 shows the volume per formula unit under two approximations between 100 K and 1000 K. The coefficient of thermal expansion under molecular dynamics approximation is 3.2024 x10-5 K-1, which is larger than that under quasi-harmonic approximation. However, this number is closer to the measured literature value. From the diagram, it is obvious that the shape of the line of quasi-harmonic is flat at first and becomes linear afterwards. This is because at very low temperatures(0 K to 100 K), the Helmholtz free energy is dominated by zero point energy. And after that, the temperature-dependent part of equation 1 dominates. The points from MD are scattered evenly near the best fit line from the very beginning.

The difference between Fig18 and Fig.19 is the range of temperature. Fig.19 expand the temperature range to 4000 K and the MD line still shows well-scattered points at high temperatures. Therefore, at low temperatures, there is only slight difference between QHA and MD approximation; but, at higher temperatures, the QHA fails and the MD approximation works. The cell volume in the QHA is larger than expected and that in the MD approximation is anticipated using its thermal expansion coefficient. Nevertheless, the melting temperature of MgO is 3125 K, and bonds break after that point. The MD approximation shows that even after melting point of MgO, there is still no bond-breaking process. This indicates that MD approximation is only appropriate below melting point of the material, not beyond.

In order to compute an accurate Free Energy, the best situation is to use a cell size as large as possible, like an infinite cell for MD calculations. This is because in small cells, all atoms move in the same direction and distance when given an energy. This is not the real life situation, where atoms move randomly and freely. A large cell averages out their motion and give a unifying result. However, this is not realistic because the calculation will last infinitely. Therefore, the alternative is to start calculation with small sizes of MgO and then increase the cell size to see how the free energy goes. The number should converge and at one cell size, the free energy won't change much. This size is a good choice to perform reliably MD for MgO.

= Conclusion =

give a general description of your calculations and your main findings

outline the differences between the methods in use and the results obtained

analyse critically these differences

= References =
# Hoffmann, R. (1987), How Chemistry and Physics Meet in the Solid State. Angew. Chem. Int. Ed. Engl., 26: 846–878. doi:10.1002/anie.198708461
# <nowiki>http://indico.ictp.it/event/7921/session/323/contribution/1273/material/0/0.pdf</nowiki>
# A. S. Madhusudhan Rao and K. Narender, “Studies on Thermophysical Properties of CaO and MgO by γ-Ray Attenuation,” Journal of Thermodynamics, vol. 2014, Article ID 123478, 8 pages, 2014. doi:10.1155/2014/123478
# S. P. Ong, S. Cholia, A. Jain, M. Brafman, D. Gunter, G. Ceder, and K. A. Persson, Computational Materials Science, 2015, 97, 209–215. doi:10.1016/j.commatsci.2014.10.037
# https://chem.libretexts.org/Textbook_Maps/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_(McQuarrie_and_Simon)/05._The_Harmonic_Oscillator_and_the_Rigid_Rotator%3A_Two_Spectroscopic_Models/5.3%3A_The_Harmonic_Oscillator_is_an_Approximation

File:Dispersion x 1x1x1(em2015).png

2017-11-18T20:46:10Z

Em2015: Em2015 uploaded a new version of File:Dispersion x 1x1x1(em2015).png

Rep:MgO(em2015)

2017-11-18T20:38:23Z

Em2015: /* Phonon Dispersion, DOS and Grid size */

=Introduction=
The experiment aims to investigate the coefficient of thermal expansion α of a crystal of MgO, using quasi-harmonic and molecular dynamics approximations under RedHat Linux environment. This intrinsic property, α is calculated as 2.6544 x10-5 K-1 and 3.2024 x10-5 K-1 between 100 K and 1000 K for each approximation. The interface is DLV, a package for the visualisation of materials, structures and properties. The calculations are completed by a classical simulation program GULP in DLVisualize.

Energies, vibrations and lattice constants of MgO crystal are computed and are used to plot graphs to decide the thermal expansion coefficient. Vibrations are important in studying thermal properties, phase transitions, electrical properties,etc. The graph of vibrational frequency against wave vector k, generated by the computer, is called Phonon dispersion. The wave vector k is a vector which points in the direction of propagation of the wave. It defines the vector in reciprocal space and is calculated as following:

<math>k(cm^{-1})=\frac{2\pi}{\lambda}</math>(1), where λ is wavelength.

The phonon dispersion diagram contains acoustical and optical branches, sharing a similar concept as band structures[1]. For the three-dimensional crystal lattice that has two atoms Mg and O in its primitive cell, the dispersion relations exhibit a total of six branches: two types of phonons in x, y, z directions. The boundaries at −π/''a'' and π/''a'' are called Brillouin zone.

Vibrational density of states is also related. it is an average over all k-points yielding the number of vibrational modes at each frequency (the density of modes). The grid size chosen in the calculation has a strong impact on the shape: generally small grids give discrete DOS and large grids give continuous DOS.

The quasi-harmonic approximation(QHA) is a more adapted hypothesis than the harmonic approximation. The system is equivalent to a collection of independent harmonic oscillators and treats vibrations as if they did not interact. It is not appropriate at high temperatures because of phonon-phonon interactions[2]. It has limitation that it does not predict bond dissociation. QHA is used to compute vibrational Helmholtz energy of MgO, with equation as following:

<math>F=E_0+\frac{1}{2}\sum_{k,j}\hslash\omega_{k,j}+k_BT\sum_{k,j}ln[1-exp(-\hslash\omega_{k,j}/k_BT)]</math> (2)

where E0 is the internal energy; the sum of ħω for each k, j is the zero-point energy; and the last term is the temperature-dependent free energy. Zero-point energy is the lowest possible energy that a quantum mechanical system may have.

The molecular dynamics approximation(MD) can capture the full anharmonicity at high temperatures. It simulates the vibrations as random motions of atoms inside a supercell. The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic evolution of the system. The system using classical Newton's Law:

<math>F=ma</math> (3)

In the calculation, energy is minimized as a function of the atomic position and of the lattice parameters.

= An Initial Calculation for MgO on its Structure and Internal Energy =
[[File:Primitive cell(em2015).PNG|left|200x200px|thumb|Fig.1 The structure of a primitive unit cell of MgO. ]]
[[File:Conventional cell(em2015).PNG|thumb|200x200px|Fig.2 conventional cell of MgO crystal]]
MgO crystal is a face-centered cube as shown on the right hand side. The conventional cell has eight atoms in it.

In the experiment, its crystal structure is described using a primitive unit cell, with cell parameters :a=b=c=2.9783Å ; α=β=γ=60.0000'''°'''. The primitive cell include one Magnesium ion (each Magnesium ion at the corner of the cell is shared by seven other cells and there are eight of this Magnesium ions) and one Oxygen ion(each Oxygen ion at the center of the cell is shared by one cell only and there is one of this Oxygen ions). The ratio of atoms in conventional cell and that in primitive cell is 4:1, which indicates that the ratio of volume between conventional cell and primitive cell is also 4:1.

The cell vectors of this cell are (in Angstrom):

0.000000 2.105970 2.105970

2.105970 0.000000 2.105970

2.105970 2.105970 0.000000.

The inter-atomic potential is 6.72119980 eV, and the total energy due to the inter-atomic forces should be -41.07531759 eV per primitive unit cell; it is the "binding energy" of the crystal - ie: the energy required to pull all of the ions of the crystal to infinite separation.

The '''Include Shells''' option should ''NOT'' be selected during calculation settings. This is because each ion is considered as a non-polarisable hard sphere here.

= Phonon Dispersion, DOS and Grid size =
Here the free energy is computed for increasing sizes of grid and its convergence to the infinite grid value monitored.

The special points along the conventional path in k-space are W-L-G-W-X-K. They are used because it is difficult to display the variation of the frequency over the full 3D '''k'''-space.[[File:Dispersion x 1x1x1(em2015).png|left|thumb|560x560px|Fig.3&4 The left hand side is the phonon dispersion diagram and the right hand side is the DOS diagram with 1x1x1 grid size.]]
{| class="wikitable"
!path in k-space
!W
!L
!G
!X
!K
|-
|kx
|1/2
|1/2
|0
|1/2
|3/8
|-
|ky
|1/4
|1/2
|0
|0
|3/8
|-
|kz
|3/4
|1/2
|0
|1/2
|3/4
|}
As mentioned in the introduction, the phonon dispersion diagram of MgO should have six branches. The vertical axis is frequency is vibration and the horizontal axis represents the path in k-space. Γ is the origin in the reciprocal space, just like 0 in the direct space.

The 1x1x1 grid size gives 4 very sharp peaks which indicates two similar higher density of states at around 280 cm-1 and 350 cm-1, and two similar lower density of states at around 670 cm-1 and 810 cm-1. The higher peaks and the lower peaks have a approximate ratio of 2:1 in height.

Furthermore, the DOS diagram can be correlated to the dispersion diagram. The point L in the dispersion diagram intersects six branches, each at 810 cm-1, 670 cm-1, 350 cm-1,350 cm-1,280 cm-1 and 280 cm-1. The two 350 cm-1 frequencies are degenerate at point L and so are the two 280 cm-1 frequencies. This, corresponds to the DOS diagram, is the 2:1 ratio in height.

The L point is 1/2 1/2 1/2.

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

Thus the K point is 1/2 1/2 1/2. If check with the log. file:

'''K point 1 = 0.500000 0.500000 0.500000. '''

This means that the speculation of the relationship between phonon dispersion diagram and DOS diagram is correct.

<nowiki> </nowiki>In order to use more '''k''' points, the the density of states calculation is reran, increasing the shrinking factors from 1x1x1, to 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 24x24x24, 36x36x36, 64x64x64. And the density of states diagrams are shown below.
{| class="wikitable"
![[File:DOS2x2x2(em2015).PNG|none|thumb|220x220px|Fig.5 Grid size: 2x2x2]]
![[File:DOS3x3x3(em2015).PNG|none|thumb|223x223px|Fig.6 Grid size: 3x3x3]]
![[File:DOS4x4x4(em2015).PNG|none|thumb|237x237px|Fig.7 Grid size: 4x4x4]]
![[File:DOS8x8x8(em2015).PNG|none|thumb|222x222px|Fig.8 Grid size: 8x8x8]]
|-
|[[File:DOS16x16x16(em2015).PNG|none|thumb|222x222px|'''Fig.9 Grid size: 16x16x16''']]
|[[File:DOS24x24x24(em2015).PNG|none|thumb|221x221px|'''Fig.10 Grid size: 24x24x24''']]
|[[File:DOS32x32x32(em2015).PNG|none|thumb|232x232px|'''Fig.11 Grid size: 36x36x36''']]
|[[File:DOS64x64x64(em2015).PNG|none|thumb|224x224px|'''Fig.12 Grid size: 64x64x64''']]
|}

When increasing the grid size, there are more peaks. When using a 8x8x8 grid size, the peaks connect to each other, giving continuous segments of lines. More and more of the possible vibrations are sampled and more features appear. When higher grid size is applied, the shape becomes more and more smooth and no apparent changes in the shape of the diagram after grid size 24x24x24. This means that the 24x24x24 grid size is the minimum/optimum.

The optimal grid size for MgO is appropriate for a calculation on a similar oxide like calcium oxide. The distance between atoms a is similar for Ca-O and Mg-O.

According to the equation of reciprocal space,

<math>a^{\ast} = \frac{2\Pi}{a}</math>(4), a* is similar as a is similar.

For Zeolite, which its primitive cell contains more atoms resulting in more branches, and larger lattice size (15.56 Å compared to 2.9783 Å of MgO) resulting in smaller reciprocal space. Thus a smaller grid size is enough to catch all the features from its phonon dispersion diagram. This is just opposite to a metal like lithium. It has only one atom in its primitive cell so only three branches in the phonon dispersion diagram. Lithium also has smaller lattice size, giving a larger a*. To put all factors together, lithium needs larger optimal grid size to include all features.

= Free Energy in the Harmonic Approximation =
The free energies of the different grid sizes in last section are reported here.
{| class="wikitable"
|+Table. 1.
!Grid Size
!Free Energy/eV
!Difference in energy/eV
!Accuracy
!Comments
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|<nowiki>-</nowiki>
|
| rowspan="9" |How does the free energy vary with grid size ?
As grid size increases, the free energy increases as well.

|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|<nowiki>-0.003692</nowiki>
|
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|0.000177
|accurate to 1 meV
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|<nowiki>-0.000018</nowiki>
|
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|<nowiki>-0.000028</nowiki>
|accurate to 0.5 meV
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|<nowiki>-0.000004</nowiki>
|accurate to 0.1 meV
|-
|24x24x24
|<nowiki>-40.926483</nowiki>
|<nowiki>-0.000001</nowiki>
|
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0.000000
|
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|0.000000
|}
The optimal grid size speculation is similar to that in the previous section.

= The Thermal Expansion Coefficient of MgO under Quasi-harmonic Approximation =
[[File:Energy vs t(em2015).PNG|thumb|592x592px|Fig.13 Graph of free energy against temperature between 0K and 1000K]]
[[File:Lattice vs t(em2015).PNG|thumb|591x591px|Fig.14 Graph of lattice constant against temperature between 0K and 1000K]]
Vary the temperature from 0K to 1000K in steps of 100K and compute the variation in the free energy and optimal lattice constant by optimizing the structure at each temperature. The numbers are extracted from the log. files and organized in the table 2.
{| class="wikitable"
|+Table. 2. The table of free energy and lattice constant at temperatures between 0K and 1000K with a variation of 100K
!Temperature/K
!Free Energy/eV
!Lattice Constant/Å
!Volume/Å3
|-
|0
|<nowiki>-40.90190628</nowiki>
|2.986563
|26.638824
|-
|100
|<nowiki>-40.90241965</nowiki>
|2.986563
|26.638824
|-
|200
|<nowiki>-40.90937738</nowiki>
|2.987605
|26.666716
|-
|300
|<nowiki>-40.92812471</nowiki>
|2.989390
|26.714542
|-
|400
|<nowiki>-40.95859417</nowiki>
|2.991629
|26.774613
|-
|500
|<nowiki>-40.99943595</nowiki>
|2.994135
|26.841955
|-
|600
|<nowiki>-41.04931542</nowiki>
|2.996820
|26.914231
|-
|700
|<nowiki>-41.10711924</nowiki>
|2.999643
|26.990362
|-
|800
|<nowiki>-41.17189187</nowiki>
|3.002588
|27.069936
|-
|900
|<nowiki>-41.24301814</nowiki>
|3.005635
|27.152431
|-
|1000
|<nowiki>-41.31984836</nowiki>
|3.008784
|27.237863
|}
The graph of free energy against temperature is shown in Fig.13, and the plot of lattice constant against temperature is shown in Fig.14. Fig.13 contains a downward curve and its gradient increases as temperature increases. The gradient becomes constant at higher temperatures, which indicated a straight line. Fig. 14 contains an upward curve and its gradient increases as temperature increases. The line becomes straight after 400K.
[[File:HightT(em2015).png|thumb|589x589px|Fig.15 Volume per formula unit at high temperatures]]
The coefficient of thermal expansion of a material is :

<math>\alpha=\frac{1}{Vo}(\frac{dV}{dT})P</math>= 2.6544 x10-5 K-1. (5)

The coefficient of volume thermal expansion is 4.0 × 10−5 K−1[2] in the temperature range from 300 K to 1250 K from A. S. Madhusudhan Rao and K. Narender's paper, which is larger than the calculated value.

For the harmonic approximation, the bond between atoms will eventually break if temperature is high enough, which indicates that the equilibrium distance between atoms is independent of temperature.The quasi-harmonic approximation expands upon the harmonic phonon model of lattice dynamics. In the quasi-harmonic approximation, bonds will never break but lengthen to infinite. In the diagrams below, the center of bonds are indicated in red in each approximations.
{| class="wikitable"
|[[File:Mp-5598-qha(em2015).png|thumb|200x200px|Fig.16 The Quasi-harmonic approximation[3] graph where the center of the bond changes. ]]
|[[File:800px-Morse-potential(em2015).png|thumb|217x217px|Fig.17 The harmonic approximation[4] graph where the center of the bond remains. ]]
|}
When a material is heated, the kinetic energy of that material is increased and it's atoms and molecules move about more. This means that each atom will take up more space due to it's movement so the material will expand.

Fig.15 is the graph of volume per formula unit under QHA at high temperatures. The blue points are the ones used to calculate the thermal expansion coefficient. The It is clear in the graph that the line deviates from the best fit line at high temperatures. Since the QHA assumes no bond-breaking between atoms even at high temperatures, the line should be linear. Thus the approximation is invalid at high temperatures.

= The Thermal Expansion Coefficient of MgO under Molecular Dynamics Approximation =
[[File:QH-MD(em2015).png|thumb|443x443px|Fig.18 Graph of volume per formula unit under Quasi-harmonic and Molecular Dynamics Approximation against temperature|left]][[File:QHA VS MD(em2015).png|thumb|497x497px|centre|Fig.19 Graph of Volume per formula unit against temperature of QHA and MD at high temperatures]]Fig.18 shows the volume per formula unit under two approximations between 100 K and 1000 K. The coefficient of thermal expansion under molecular dynamics approximation is 3.2024 x10-5 K-1, which is larger than that under quasi-harmonic approximation. However, this number is closer to the measured literature value. From the diagram, it is obvious that the shape of the line of quasi-harmonic is flat at first and becomes linear afterwards. This is because at very low temperatures(0 K to 100 K), the Helmholtz free energy is dominated by zero point energy. And after that, the temperature-dependent part of equation 1 dominates. The points from MD are scattered evenly near the best fit line from the very beginning.

The difference between Fig18 and Fig.19 is the range of temperature. Fig.19 expand the temperature range to 4000 K and the MD line still shows well-scattered points at high temperatures. Therefore, at low temperatures, there is only slight difference between QHA and MD approximation; but, at higher temperatures, the QHA fails and the MD approximation works. The cell volume in the QHA is larger than expected and that in the MD approximation is anticipated using its thermal expansion coefficient. Nevertheless, the melting temperature of MgO is 3125 K, and bonds break after that point. The MD approximation shows that even after melting point of MgO, there is still no bond-breaking process. This indicates that MD approximation is only appropriate below melting point of the material, not beyond.

In order to compute an accurate Free Energy, the best situation is to use a cell size as large as possible, like an infinite cell for MD calculations. This is because in small cells, all atoms move in the same direction and distance when given an energy. This is not the real life situation, where atoms move randomly and freely. A large cell averages out their motion and give a unifying result. However, this is not realistic because the calculation will last infinitely. Therefore, the alternative is to start calculation with small sizes of MgO and then increase the cell size to see how the free energy goes. The number should converge and at one cell size, the free energy won't change much. This size is a good choice to perform reliably MD for MgO.

= Conclusion =

give a general description of your calculations and your main findings

outline the differences between the methods in use and the results obtained

analyse critically these differences

= References =
# Hoffmann, R. (1987), How Chemistry and Physics Meet in the Solid State. Angew. Chem. Int. Ed. Engl., 26: 846–878. doi:10.1002/anie.198708461
# <nowiki>http://indico.ictp.it/event/7921/session/323/contribution/1273/material/0/0.pdf</nowiki>
# A. S. Madhusudhan Rao and K. Narender, “Studies on Thermophysical Properties of CaO and MgO by γ-Ray Attenuation,” Journal of Thermodynamics, vol. 2014, Article ID 123478, 8 pages, 2014. doi:10.1155/2014/123478
# S. P. Ong, S. Cholia, A. Jain, M. Brafman, D. Gunter, G. Ceder, and K. A. Persson, Computational Materials Science, 2015, 97, 209–215. doi:10.1016/j.commatsci.2014.10.037
# https://chem.libretexts.org/Textbook_Maps/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_(McQuarrie_and_Simon)/05._The_Harmonic_Oscillator_and_the_Rigid_Rotator%3A_Two_Spectroscopic_Models/5.3%3A_The_Harmonic_Oscillator_is_an_Approximation

Rep:MgO(em2015)

2017-11-18T20:31:31Z

Em2015: /* Lattice Vibrations - Computing the Phonons */

=Introduction=
The experiment aims to investigate the coefficient of thermal expansion α of a crystal of MgO, using quasi-harmonic and molecular dynamics approximations under RedHat Linux environment. This intrinsic property, α is calculated as 2.6544 x10-5 K-1 and 3.2024 x10-5 K-1 between 100 K and 1000 K for each approximation. The interface is DLV, a package for the visualisation of materials, structures and properties. The calculations are completed by a classical simulation program GULP in DLVisualize.

Energies, vibrations and lattice constants of MgO crystal are computed and are used to plot graphs to decide the thermal expansion coefficient. Vibrations are important in studying thermal properties, phase transitions, electrical properties,etc. The graph of vibrational frequency against wave vector k, generated by the computer, is called Phonon dispersion. The wave vector k is a vector which points in the direction of propagation of the wave. It defines the vector in reciprocal space and is calculated as following:

<math>k(cm^{-1})=\frac{2\pi}{\lambda}</math>(1), where λ is wavelength.

The phonon dispersion diagram contains acoustical and optical branches, sharing a similar concept as band structures[1]. For the three-dimensional crystal lattice that has two atoms Mg and O in its primitive cell, the dispersion relations exhibit a total of six branches: two types of phonons in x, y, z directions. The boundaries at −π/''a'' and π/''a'' are called Brillouin zone.

Vibrational density of states is also related. it is an average over all k-points yielding the number of vibrational modes at each frequency (the density of modes). The grid size chosen in the calculation has a strong impact on the shape: generally small grids give discrete DOS and large grids give continuous DOS.

The quasi-harmonic approximation(QHA) is a more adapted hypothesis than the harmonic approximation. The system is equivalent to a collection of independent harmonic oscillators and treats vibrations as if they did not interact. It is not appropriate at high temperatures because of phonon-phonon interactions[2]. It has limitation that it does not predict bond dissociation. QHA is used to compute vibrational Helmholtz energy of MgO, with equation as following:

<math>F=E_0+\frac{1}{2}\sum_{k,j}\hslash\omega_{k,j}+k_BT\sum_{k,j}ln[1-exp(-\hslash\omega_{k,j}/k_BT)]</math> (2)

where E0 is the internal energy; the sum of ħω for each k, j is the zero-point energy; and the last term is the temperature-dependent free energy. Zero-point energy is the lowest possible energy that a quantum mechanical system may have.

The molecular dynamics approximation(MD) can capture the full anharmonicity at high temperatures. It simulates the vibrations as random motions of atoms inside a supercell. The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic evolution of the system. The system using classical Newton's Law:

<math>F=ma</math> (3)

In the calculation, energy is minimized as a function of the atomic position and of the lattice parameters.

= An Initial Calculation for MgO on its Structure and Internal Energy =
[[File:Primitive cell(em2015).PNG|left|200x200px|thumb|Fig.1 The structure of a primitive unit cell of MgO. ]]
[[File:Conventional cell(em2015).PNG|thumb|200x200px|Fig.2 conventional cell of MgO crystal]]
MgO crystal is a face-centered cube as shown on the right hand side. The conventional cell has eight atoms in it.

In the experiment, its crystal structure is described using a primitive unit cell, with cell parameters :a=b=c=2.9783Å ; α=β=γ=60.0000'''°'''. The primitive cell include one Magnesium ion (each Magnesium ion at the corner of the cell is shared by seven other cells and there are eight of this Magnesium ions) and one Oxygen ion(each Oxygen ion at the center of the cell is shared by one cell only and there is one of this Oxygen ions). The ratio of atoms in conventional cell and that in primitive cell is 4:1, which indicates that the ratio of volume between conventional cell and primitive cell is also 4:1.

The cell vectors of this cell are (in Angstrom):

0.000000 2.105970 2.105970

2.105970 0.000000 2.105970

2.105970 2.105970 0.000000.

The inter-atomic potential is 6.72119980 eV, and the total energy due to the inter-atomic forces should be -41.07531759 eV per primitive unit cell; it is the "binding energy" of the crystal - ie: the energy required to pull all of the ions of the crystal to infinite separation.

The '''Include Shells''' option should ''NOT'' be selected during calculation settings. This is because each ion is considered as a non-polarisable hard sphere here.

= Phonon Dispersion, DOS and Grid size =
Here the free energy is computed for increasing sizes of grid and its convergence to the infinite grid value monitored.

The special points along the conventional path in '''k'''-space are W-L-G-W-X-K. They are used because it is difficult to display the variation of the frequency over the full 3D '''k'''-space.
{| class="wikitable"
!path in k-space
!W
!L
!G
!X
!W
!K
|-
|kx
|1/2
|1/2
|0
|1/2
|1/2
|3/8
|-
|ky
|1/4
|1/2
|0
|0
|1/4
|3/8
|-
|kz
|3/4
|1/2
|0
|1/2
|3/4
|
|}

[[File:Dispersion x 1x1x1(em2015).png|left|thumb|560x560px|Fig.3&4 The left hand side is the phonon dispersion diagram and the right hand side is the DOS diagram with 1x1x1 grid size.]]
MgO has two ions and in a 3D structure(kx, ky, kz in reciprocal space), overall of 6 branches, as shwon in the phonon dispersion diagram.

The 1x1x1 grid size gives 4 very sharp peaks which indicates two similar higher density of states at around 280 cm-1 and 350 cm-1, and two similar lower density of states at around 670 cm-1 and 810 cm-1. The higher peaks and the lower peaks have a approximate ratio of 2:1 in height.

Furthermore, the DOS diagram can be correlated to the dispersion diagram. The point L in the dispersion diagram intersects six branches, each at 810 cm-1, 670 cm-1, 350 cm-1,350 cm-1,280 cm-1 and 280 cm-1. The two 350 cm-1 frequencies are degenerate at point L and so are the two 280 cm-1 frequencies. This, corresponds to the DOS diagram, is the 2:1 ratio in height.

The L point is 1/2 1/2 1/2.

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

Thus the K point is 1/2 1/2 1/2. If check with the log. file:

'''K point 1 = 0.500000 0.500000 0.500000. '''

This means that the speculation of the relationship between phonon dispersion diagram and DOS diagram is correct.

<nowiki> </nowiki>In order to use more '''k''' points, the the density of states calculation is reran, increasing the shrinking factors from 1x1x1, to 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 24x24x24, 36x36x36, 64x64x64. And the density of states diagrams are shown below.
{| class="wikitable"
![[File:DOS2x2x2(em2015).PNG|none|thumb|220x220px|Fig.5 Grid size: 2x2x2]]
![[File:DOS3x3x3(em2015).PNG|none|thumb|223x223px|Fig.6 Grid size: 3x3x3]]
![[File:DOS4x4x4(em2015).PNG|none|thumb|237x237px|Fig.7 Grid size: 4x4x4]]
![[File:DOS8x8x8(em2015).PNG|none|thumb|222x222px|Fig.8 Grid size: 8x8x8]]
|-
|[[File:DOS16x16x16(em2015).PNG|none|thumb|222x222px|'''Fig.9 Grid size: 16x16x16''']]
|[[File:DOS24x24x24(em2015).PNG|none|thumb|221x221px|'''Fig.10 Grid size: 24x24x24''']]
|[[File:DOS32x32x32(em2015).PNG|none|thumb|232x232px|'''Fig.11 Grid size: 36x36x36''']]
|[[File:DOS64x64x64(em2015).PNG|none|thumb|224x224px|'''Fig.12 Grid size: 64x64x64''']]
|}

When increasing the grid size, there are more peaks. When using a 8x8x8 grid size, the peaks connect to each other, giving continuous segments of lines. More and more of the possible vibrations are sampled and more features appear. When higher grid size is applied, the shape becomes more and more smooth and no apparent changes in the shape of the diagram after grid size 24x24x24. This means that the 24x24x24 grid size is the minimum/optimum.

The optimal grid size for MgO is appropriate for a calculation on a similar oxide like calcium oxide. The distance between atoms a is similar for Ca-O and Mg-O.

According to the equation of reciprocal space,

<math>a^{\ast} = \frac{2\Pi}{a}</math>(4), a* is similar as a is similar.

For Zeolite, which its primitive cell contains more atoms resulting in more branches, and larger lattice size (15.56 Å compared to 2.9783 Å of MgO) resulting in smaller reciprocal space. Thus a smaller grid size is enough to catch all the features from its phonon dispersion diagram. This is just opposite to a metal like lithium. It has only one atom in its primitive cell so only three branches in the phonon dispersion diagram. Lithium also has smaller lattice size, giving a larger a*. To put all factors together, lithium needs larger optimal grid size to include all features.

= Calculating the Free Energy in the Harmonic Approximation =
The free energies of the different grid sizes in last section are reported here.
{| class="wikitable"
|+Table. 1.
!Grid Size
!Free Energy/eV
!Difference in energy/eV
!Accuracy
!Comments
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|<nowiki>-</nowiki>
|
| rowspan="9" |How does the free energy vary with grid size ?
As grid size increases, the free energy increases as well.

|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|<nowiki>-0.003692</nowiki>
|
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|0.000177
|accurate to 1 meV
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|<nowiki>-0.000018</nowiki>
|
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|<nowiki>-0.000028</nowiki>
|accurate to 0.5 meV
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|<nowiki>-0.000004</nowiki>
|accurate to 0.1 meV
|-
|24x24x24
|<nowiki>-40.926483</nowiki>
|<nowiki>-0.000001</nowiki>
|
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0.000000
|
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|0.000000
|}
The optimal grid size speculation is similar to that in the previous section.

= The Thermal Expansion Coefficient of MgO under Quasi-harmonic Approximation =
[[File:Energy vs t(em2015).PNG|thumb|592x592px|Fig.13 Graph of free energy against temperature between 0K and 1000K]]
[[File:Lattice vs t(em2015).PNG|thumb|591x591px|Fig.14 Graph of lattice constant against temperature between 0K and 1000K]]
Vary the temperature from 0K to 1000K in steps of 100K and compute the variation in the free energy and optimal lattice constant by optimizing the structure at each temperature. The numbers are extracted from the log. files and organized in the table 2.
{| class="wikitable"
|+Table. 2. The table of free energy and lattice constant at temperatures between 0K and 1000K with a variation of 100K
!Temperature/K
!Free Energy/eV
!Lattice Constant/Å
!Volume/Å3
|-
|0
|<nowiki>-40.90190628</nowiki>
|2.986563
|26.638824
|-
|100
|<nowiki>-40.90241965</nowiki>
|2.986563
|26.638824
|-
|200
|<nowiki>-40.90937738</nowiki>
|2.987605
|26.666716
|-
|300
|<nowiki>-40.92812471</nowiki>
|2.989390
|26.714542
|-
|400
|<nowiki>-40.95859417</nowiki>
|2.991629
|26.774613
|-
|500
|<nowiki>-40.99943595</nowiki>
|2.994135
|26.841955
|-
|600
|<nowiki>-41.04931542</nowiki>
|2.996820
|26.914231
|-
|700
|<nowiki>-41.10711924</nowiki>
|2.999643
|26.990362
|-
|800
|<nowiki>-41.17189187</nowiki>
|3.002588
|27.069936
|-
|900
|<nowiki>-41.24301814</nowiki>
|3.005635
|27.152431
|-
|1000
|<nowiki>-41.31984836</nowiki>
|3.008784
|27.237863
|}
The graph of free energy against temperature is shown in Fig.13, and the plot of lattice constant against temperature is shown in Fig.14. Fig.13 contains a downward curve and its gradient increases as temperature increases. The gradient becomes constant at higher temperatures, which indicated a straight line. Fig. 14 contains an upward curve and its gradient increases as temperature increases. The line becomes straight after 400K.
[[File:HightT(em2015).png|thumb|589x589px|Fig.15 Volume per formula unit at high temperatures]]
The coefficient of thermal expansion of a material is :

<math>\alpha=\frac{1}{Vo}(\frac{dV}{dT})P</math>= 2.6544 x10-5 K-1. (5)

The coefficient of volume thermal expansion is 4.0 × 10−5 K−1[2] in the temperature range from 300 K to 1250 K from A. S. Madhusudhan Rao and K. Narender's paper, which is larger than the calculated value.

For the harmonic approximation, the bond between atoms will eventually break if temperature is high enough, which indicates that the equilibrium distance between atoms is independent of temperature.The quasi-harmonic approximation expands upon the harmonic phonon model of lattice dynamics. In the quasi-harmonic approximation, bonds will never break but lengthen to infinite. In the diagrams below, the center of bonds are indicated in red in each approximations.
{| class="wikitable"
|[[File:Mp-5598-qha(em2015).png|thumb|200x200px|Fig.16 The Quasi-harmonic approximation[3] graph where the center of the bond changes. ]]
|[[File:800px-Morse-potential(em2015).png|thumb|217x217px|Fig.17 The harmonic approximation[4] graph where the center of the bond remains. ]]
|}
When a material is heated, the kinetic energy of that material is increased and it's atoms and molecules move about more. This means that each atom will take up more space due to it's movement so the material will expand.

Fig.15 is the graph of volume per formula unit under QHA at high temperatures. The blue points are the ones used to calculate the thermal expansion coefficient. The It is clear in the graph that the line deviates from the best fit line at high temperatures. Since the QHA assumes no bond-breaking between atoms even at high temperatures, the line should be linear. Thus the approximation is invalid at high temperatures.

= The Thermal Expansion Coefficient of MgO under Molecular Dynamics Approximation =
[[File:QH-MD(em2015).png|thumb|443x443px|Fig.18 Graph of volume per formula unit under Quasi-harmonic and Molecular Dynamics Approximation against temperature|left]][[File:QHA VS MD(em2015).png|thumb|497x497px|centre|Fig.19 Graph of Volume per formula unit against temperature of QHA and MD at high temperatures]]Fig.18 shows the volume per formula unit under two approximations between 100 K and 1000 K. The coefficient of thermal expansion under molecular dynamics approximation is 3.2024 x10-5 K-1, which is larger than that under quasi-harmonic approximation. However, this number is closer to the measured literature value. From the diagram, it is obvious that the shape of the line of quasi-harmonic is flat at first and becomes linear afterwards. This is because at very low temperatures(0 K to 100 K), the Helmholtz free energy is dominated by zero point energy. And after that, the temperature-dependent part of equation 1 dominates. The points from MD are scattered evenly near the best fit line from the very beginning.

The difference between Fig18 and Fig.19 is the range of temperature. Fig.19 expand the temperature range to 4000 K and the MD line still shows well-scattered points at high temperatures. Therefore, at low temperatures, there is only slight difference between QHA and MD approximation; but, at higher temperatures, the QHA fails and the MD approximation works. The cell volume in the QHA is larger than expected and that in the MD approximation is anticipated using its thermal expansion coefficient. Nevertheless, the melting temperature of MgO is 3125 K, and bonds break after that point. The MD approximation shows that even after melting point of MgO, there is still no bond-breaking process. This indicates that MD approximation is only appropriate below melting point of the material, not beyond.

In order to compute an accurate Free Energy, the best situation is to use a cell size as large as possible, like an infinite cell for MD calculations. This is because in small cells, all atoms move in the same direction and distance when given an energy. This is not the real life situation, where atoms move randomly and freely. A large cell averages out their motion and give a unifying result. However, this is not realistic because the calculation will last infinitely. Therefore, the alternative is to start calculation with small sizes of MgO and then increase the cell size to see how the free energy goes. The number should converge and at one cell size, the free energy won't change much. This size is a good choice to perform reliably MD for MgO.

= Conclusion =

give a general description of your calculations and your main findings

outline the differences between the methods in use and the results obtained

analyse critically these differences

= References =
# Hoffmann, R. (1987), How Chemistry and Physics Meet in the Solid State. Angew. Chem. Int. Ed. Engl., 26: 846–878. doi:10.1002/anie.198708461
# <nowiki>http://indico.ictp.it/event/7921/session/323/contribution/1273/material/0/0.pdf</nowiki>
# A. S. Madhusudhan Rao and K. Narender, “Studies on Thermophysical Properties of CaO and MgO by γ-Ray Attenuation,” Journal of Thermodynamics, vol. 2014, Article ID 123478, 8 pages, 2014. doi:10.1155/2014/123478
# S. P. Ong, S. Cholia, A. Jain, M. Brafman, D. Gunter, G. Ceder, and K. A. Persson, Computational Materials Science, 2015, 97, 209–215. doi:10.1016/j.commatsci.2014.10.037
# https://chem.libretexts.org/Textbook_Maps/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_(McQuarrie_and_Simon)/05._The_Harmonic_Oscillator_and_the_Rigid_Rotator%3A_Two_Spectroscopic_Models/5.3%3A_The_Harmonic_Oscillator_is_an_Approximation

Rep:MgO(em2015)

2017-11-18T20:30:01Z

2017-11-17T14:12:57Z

Em2015: /* Conclusion */

=Introduction=
The aim of the experiment is to investigate an intrinsic property, coefficient of thermal expansion α of a crystal of MgO, using quasi-harmonic and molecular dynamics approximations under RedHat Linux environment. And α is calculated as 2.6544 x10-5 K-1 between 100 K and 1000 K. The interface is DLV, a package for the visualisation of materials, structures and properties. The calculations are completed by a classical simulation program GULP in DLVisualize.

During the investigation of thermal expansion coefficient, the energies and vibrations of MgO crystal are also explored. Vibrations are important in studying thermal properties, phase transitions, electrical properties,etc. Phonon dispersion, which shares a similar concept as band structures[1], is a graph of the relationship between vibration frequency and energy. Vibrational density of states is also related. It shows the relative ratio of states at each vibration frequencies and it's derived from phonon dispersion. The grid size chosen in the calculation has a strong impact on the shape.

The quasi-harmonic approximation(QHA) is a more adapted hypothesis than the harmonic approximation. They have the same limitation that they do not predict bond dissociation. QHA is used to compute vibrational energy levels of MgO, with equation as following: <math>F=E_0+\frac{1}{2}\sum_{k,j}\hslash\omega_{k,j}+k_BT\sum_{k,j}ln[1-exp(-\hslash\omega_{k,j}/k_BT)]</math> (1)

where E0 is the internal energy; the sum of ħω for each k, j is the zero-point energy; and the last term is the temperature-dependent free energy. Zero-point energy is the lowest possible energy that a quantum mechanical system may have.

The molecular dynamics approximation simulates the vibrations as random motions of atoms inside a supercell, using classical Newton's Law:

<math>F=ma</math> (2)

The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic evolution of the system.

The quasi-harmonic approximation and molecular dynamics approximation are compared and contrasted at both low temperatures and high temperatures, to explore which methodology is more appropriate for thermal expansion modelling with MgO.

= An Initial Calculation for MgO on its Structure and Internal Energy =
[[File:Primitive cell(em2015).PNG|left|160x160px|thumb|Fig.1 The structure of a primitive unit cell of MgO. ]]
[[File:Conventional cell(em2015).PNG|thumb|179x179px|Fig.2 conventional cell of MgO crystal]]
MgO crystal is a face-centered cube as shown on the right hand side. In the experiment, its crystal structure is described using a primitive unit cell, with cell parameters :a=b=c=2.9783Å ; α=β=γ=60.0000'''°'''. The primitive cell include one magnesium ion (1/8 x 8)and one oxygen ion(1 x 1). Conventional cell has 8 atoms, which indicates that the ratio of volume between conventional cell and primitive cell is 4:1.

The cell vectors of this cell are (in Angstrom):

0.000000 2.105970 2.105970

2.105970 0.000000 2.105970

2.105970 2.105970 0.000000.

The inter-atomic potential is 6.72119980 eV, and the total energy due to the inter-atomic forces should be -41.07531759 eV per primitive unit cell; it is the "binding energy" of the crystal - ie: the energy required to pull all of the ions of the crystal to infinite separation.

The '''Include Shells''' option should ''NOT'' be selected during calculation settings. This is because each ion is considered as a non-polarisable hard sphere here.

= Lattice Vibrations - Computing the Phonons =
Here the free energy will be computed for ever increasing sizes of grid and its convergence to the infinite grid value monitored. The convergence of the density of states plot in the previous exercise is a useful guide for choosing the optimal grid size.

== The relationship between phonon dispersion and DOS ==
[[File:Dispersion x 1x1x1(em2015).png|left|thumb|560x560px|Fig.3&4 The left hand side is the phonon dispersion diagram and the right hand side is the DOS diagram with 1x1x1 grid size.]]
MgO has two ions and in a 3D structure(kx, ky, kz in reciprocal space), overall of 6 branches, as shwon in the phonon dispersion diagram.

The 1x1x1 grid size gives 4 very sharp peaks which indicates two similar higher density of states at around 280 cm-1 and 350 cm-1, and two similar lower density of states at around 670 cm-1 and 810 cm-1. The higher peaks and the lower peaks have a approximate ratio of 2:1 in height.

Furthermore, the DOS diagram can be correlated to the dispersion diagram. The point L in the dispersion diagram intersects six branches, each at 810 cm-1, 670 cm-1, 350 cm-1,350 cm-1,280 cm-1 and 280 cm-1. The two 350 cm-1 frequencies are degenerate at point L and so are the two 280 cm-1 frequencies. This, corresponds to the DOS diagram, is the 2:1 ratio in height.

The L point is 1/2 1/2 1/2.

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

Thus the K point is 1/2 1/2 1/2. If check with the log. file:

'''K point 1 = 0.500000 0.500000 0.500000. '''

This means that the speculation of the relationship between phonon dispersion diagram and DOS diagram is correct.

== Density of States and Grid Size ==
In order to use more '''k''' points, the the density of states calculation is reran, increasing the shrinking factors from 1x1x1, to 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 24x24x24, 36x36x36, 64x64x64. And the density of states diagrams are shown below.
{| class="wikitable"
![[File:DOS2x2x2(em2015).PNG|none|thumb|220x220px|Fig.5 Grid size: 2x2x2]]
![[File:DOS3x3x3(em2015).PNG|none|thumb|223x223px|Fig.6 Grid size: 3x3x3]]
![[File:DOS4x4x4(em2015).PNG|none|thumb|237x237px|Fig.7 Grid size: 4x4x4]]
![[File:DOS8x8x8(em2015).PNG|none|thumb|222x222px|Fig.8 Grid size: 8x8x8]]
|-
|[[File:DOS16x16x16(em2015).PNG|none|thumb|222x222px|'''Fig.9 Grid size: 16x16x16''']]
|[[File:DOS24x24x24(em2015).PNG|none|thumb|221x221px|'''Fig.10 Grid size: 24x24x24''']]
|[[File:DOS32x32x32(em2015).PNG|none|thumb|232x232px|'''Fig.11 Grid size: 36x36x36''']]
|[[File:DOS64x64x64(em2015).PNG|none|thumb|224x224px|'''Fig.12 Grid size: 64x64x64''']]
|}

When increasing the grid size, there are more peaks. When using a 8x8x8 grid size, the peaks connect to each other, giving continuous segments of lines. More and more of the possible vibrations are sampled and more features appear. When higher grid size is applied, the shape becomes more and more smooth and no apparent changes in the shape of the diagram after grid size 24x24x24. This means that the 24x24x24 grid size is the minimum/optimum.

The optimal grid size for MgO is appropriate for a calculation on a similar oxide like calcium oxide. The distance between atoms a is similar for Ca-O and Mg-O.

According to the equation of reciprocal space,

<math>a^{\ast} = \frac{2\Pi}{a}</math>, a* is similar as a is similar.

For Zeolite, which its primitive cell contains more atoms resulting in more branches, and larger lattice size (15.56 Å compared to 2.9783 Å of MgO) resulting in smaller reciprocal space. Thus a smaller grid size is enough to catch all the features from its phonon dispersion diagram. This is just opposite to a metal like lithium. It has only one atom in its primitive cell so only three branches in the phonon dispersion diagram. Lithium also has smaller lattice size, giving a larger a*. To put all factors together, lithium needs larger optimal grid size to include all features.

= Calculating the Free Energy in the Harmonic Approximation =
The free energies of the different grid sizes in last section are reported here.
{| class="wikitable"
|+Table. 1.
!Grid Size
!Free Energy/eV
!Difference in energy/eV
!Accuracy
!Comments
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|<nowiki>-</nowiki>
|
| rowspan="9" |How does the free energy vary with grid size ?
As grid size increases, the free energy increases as well.

|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|<nowiki>-0.003692</nowiki>
|
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|0.000177
|accurate to 1 meV
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|<nowiki>-0.000018</nowiki>
|
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|<nowiki>-0.000028</nowiki>
|accurate to 0.5 meV
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|<nowiki>-0.000004</nowiki>
|accurate to 0.1 meV
|-
|24x24x24
|<nowiki>-40.926483</nowiki>
|<nowiki>-0.000001</nowiki>
|
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0.000000
|
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|0.000000
|}
The optimal grid size speculation is similar to that in the previous section.

= The Thermal Expansion Coefficient of MgO under Quasi-harmonic Approximation =
[[File:Energy vs t(em2015).PNG|thumb|592x592px|Fig.13 Graph of free energy against temperature between 0K and 1000K]]
[[File:Lattice vs t(em2015).PNG|thumb|591x591px|Fig.14 Graph of lattice constant against temperature between 0K and 1000K]]
Vary the temperature from 0K to 1000K in steps of 100K and compute the variation in the free energy and optimal lattice constant by optimizing the structure at each temperature. The numbers are extracted from the log. files and organized in the table 2.
{| class="wikitable"
|+Table. 2. The table of free energy and lattice constant at temperatures between 0K and 1000K with a variation of 100K
!Temperature/K
!Free Energy/eV
!Lattice Constant/Å
!Volume/Å3
|-
|0
|<nowiki>-40.90190628</nowiki>
|2.986563
|26.638824
|-
|100
|<nowiki>-40.90241965</nowiki>
|2.986563
|26.638824
|-
|200
|<nowiki>-40.90937738</nowiki>
|2.987605
|26.666716
|-
|300
|<nowiki>-40.92812471</nowiki>
|2.989390
|26.714542
|-
|400
|<nowiki>-40.95859417</nowiki>
|2.991629
|26.774613
|-
|500
|<nowiki>-40.99943595</nowiki>
|2.994135
|26.841955
|-
|600
|<nowiki>-41.04931542</nowiki>
|2.996820
|26.914231
|-
|700
|<nowiki>-41.10711924</nowiki>
|2.999643
|26.990362
|-
|800
|<nowiki>-41.17189187</nowiki>
|3.002588
|27.069936
|-
|900
|<nowiki>-41.24301814</nowiki>
|3.005635
|27.152431
|-
|1000
|<nowiki>-41.31984836</nowiki>
|3.008784
|27.237863
|}
The graph of free energy against temperature is shown in Fig.13, and the plot of lattice constant against temperature is shown in Fig.14. Fig.13 contains a downward curve and its gradient increases as temperature increases. The gradient becomes constant at higher temperatures, which indicated a straight line. Fig. 14 contains an upward curve and its gradient increases as temperature increases. The line becomes straight after 400K.
[[File:HightT(em2015).png|thumb|589x589px|Fig.15 Volume per formula unit at high temperatures]]
The coefficient of thermal expansion of a material is :

<math>\alpha=\frac{1}{Vo}(\frac{dV}{dT})P</math>= 2.6544 x10-5 K-1. (3)

The coefficient of volume thermal expansion is 4.0 × 10−5 K−1[2] in the temperature range from 300 K to 1250 K from A. S. Madhusudhan Rao and K. Narender's paper, which is larger than the calculated value.

For the harmonic approximation, the bond between atoms will eventually break if temperature is high enough, which indicates that the equilibrium distance between atoms is independent of temperature.The quasi-harmonic approximation expands upon the harmonic phonon model of lattice dynamics. In the quasi-harmonic approximation, bonds will never break but lengthen to infinite. In the diagrams below, the center of bonds are indicated in red in each approximations.
{| class="wikitable"
|[[File:Mp-5598-qha(em2015).png|thumb|200x200px|Fig.16 The Quasi-harmonic approximation[3] graph where the center of the bond changes. ]]
|[[File:800px-Morse-potential(em2015).png|thumb|217x217px|Fig.17 The harmonic approximation[4] graph where the center of the bond remains. ]]
|}
When a material is heated, the kinetic energy of that material is increased and it's atoms and molecules move about more. This means that each atom will take up more space due to it's movement so the material will expand.

Fig.15 is the graph of volume per formula unit under QHA at high temperatures. The blue points are the ones used to calculate the thermal expansion coefficient. The It is clear in the graph that the line deviates from the best fit line at high temperatures. Since the QHA assumes no bond-breaking between atoms even at high temperatures, the line should be linear. Thus the approximation is invalid at high temperatures.

= The Thermal Expansion Coefficient of MgO under Molecular Dynamics Approximation =
[[File:QH-MD(em2015).png|thumb|443x443px|Fig.18 Graph of volume per formula unit under Quasi-harmonic and Molecular Dynamics Approximation against temperature|left]][[File:QHA VS MD(em2015).png|thumb|497x497px|centre|Fig.19 Graph of Volume per formula unit against temperature of QHA and MD at high temperatures]]Fig.18 shows the volume per formula unit under two approximations between 100 K and 1000 K. The coefficient of thermal expansion under molecular dynamics approximation is 3.2024 x10-5 K-1, which is larger than that under quasi-harmonic approximation. However, this number is closer to the measured literature value. From the diagram, it is obvious that the shape of the line of quasi-harmonic is flat at first and becomes linear afterwards. This is because at very low temperatures(0 K to 100 K), the Helmholtz free energy is dominated by zero point energy. And after that, the temperature-dependent part of equation 1 dominates. The points from MD are scattered evenly near the best fit line from the very beginning.

The difference between Fig18 and Fig.19 is the range of temperature. Fig.19 expand the temperature range to 4000 K and the MD line still shows well-scattered points at high temperatures. Therefore, at low temperatures, there is only slight difference between QHA and MD approximation; but, at higher temperatures, the QHA fails and the MD approximation works. The cell volume in the QHA is larger than expected and that in the MD approximation is anticipated using its thermal expansion coefficient. Nevertheless, the melting temperature of MgO is 3125 K, and bonds break after that point. The MD approximation shows that even after melting point of MgO, there is still no bond-breaking process. This indicates that MD approximation is only appropriate below melting point of the material, not beyond.

In order to compute an accurate Free Energy, the best situation is to use a cell size as large as possible, like an infinite cell for MD calculations. This is because in small cells, all atoms move in the same direction and distance when given an energy. This is not the real life situation, where atoms move randomly and freely. A large cell averages out their motion and give a unifying result. However, this is not realistic because the calculation will last infinitely. Therefore, the alternative is to start calculation with small sizes of MgO and then increase the cell size to see how the free energy goes. The number should converge and at one cell size, the free energy won't change much. This size is a good choice to perform reliably MD for MgO.

= Conclusion =

give a general description of your calculations and your main findings

outline the differences between the methods in use and the results obtained

analyse critically these differences

= References =
# Hoffmann, R. (1987), How Chemistry and Physics Meet in the Solid State. Angew. Chem. Int. Ed. Engl., 26: 846–878. doi:10.1002/anie.198708461
# A. S. Madhusudhan Rao and K. Narender, “Studies on Thermophysical Properties of CaO and MgO by γ-Ray Attenuation,” Journal of Thermodynamics, vol. 2014, Article ID 123478, 8 pages, 2014. doi:10.1155/2014/123478
# S. P. Ong, S. Cholia, A. Jain, M. Brafman, D. Gunter, G. Ceder, and K. A. Persson, Computational Materials Science, 2015, 97, 209–215. doi:10.1016/j.commatsci.2014.10.037
# https://chem.libretexts.org/Textbook_Maps/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_(McQuarrie_and_Simon)/05._The_Harmonic_Oscillator_and_the_Rigid_Rotator%3A_Two_Spectroscopic_Models/5.3%3A_The_Harmonic_Oscillator_is_an_Approximation

Rep:MgO(em2015)

2017-11-17T13:53:43Z

Em2015: /* Density of States and Grid Size */

=Introduction=
The aim of the experiment is to investigate an intrinsic property, coefficient of thermal expansion α of a crystal of MgO, using quasi-harmonic and molecular dynamics approximations under RedHat Linux environment. And α is calculated as 2.6544 x10-5 K-1 between 100 K and 1000 K. The interface is DLV, a package for the visualisation of materials, structures and properties. The calculations are completed by a classical simulation program GULP in DLVisualize.

During the investigation of thermal expansion coefficient, the energies and vibrations of MgO crystal are also explored. Vibrations are important in studying thermal properties, phase transitions, electrical properties,etc. Phonon dispersion, which shares a similar concept as band structures[1], is a graph of the relationship between vibration frequency and energy. Vibrational density of states is also related. It shows the relative ratio of states at each vibration frequencies and it's derived from phonon dispersion. The grid size chosen in the calculation has a strong impact on the shape.

The quasi-harmonic approximation(QHA) is a more adapted hypothesis than the harmonic approximation. They have the same limitation that they do not predict bond dissociation. QHA is used to compute vibrational energy levels of MgO, with equation as following: <math>F=E_0+\frac{1}{2}\sum_{k,j}\hslash\omega_{k,j}+k_BT\sum_{k,j}ln[1-exp(-\hslash\omega_{k,j}/k_BT)]</math> (1)

where E0 is the internal energy; the sum of ħω for each k, j is the zero-point energy; and the last term is the temperature-dependent free energy. Zero-point energy is the lowest possible energy that a quantum mechanical system may have.

The molecular dynamics approximation simulates the vibrations as random motions of atoms inside a supercell, using classical Newton's Law:

<math>F=ma</math> (2)

The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic evolution of the system.

The quasi-harmonic approximation and molecular dynamics approximation are compared and contrasted at both low temperatures and high temperatures, to explore which methodology is more appropriate for thermal expansion modelling with MgO.

= Calculating the internal energy of an MgO crystal =
[[File:Primitive cell(em2015).PNG|left|160x160px|thumb|Fig.1 The structure of a primitive unit cell of MgO. ]]
[[File:Conventional cell(em2015).PNG|thumb|179x179px|Fig.2 conventional cell of MgO crystal]]
MgO crystal is a face-centered cube as shown on the right hand side. In the experiment, its crystal structure is described using a primitive unit cell, with cell parameters :a=b=c=2.9783Å ; α=β=γ=60.0000'''°'''. The primitive cell include one magnesium ion (1/8 x 8)and one oxygen ion(1 x 1). Conventional cell has 8 atoms, which indicates that the ratio of volume between conventional cell and primitive cell is 4:1.

The cell vectors of this cell are (in Angstrom):

0.000000 2.105970 2.105970

2.105970 0.000000 2.105970

2.105970 2.105970 0.000000.

The inter-atomic potential is 6.72119980 eV, and the total energy due to the inter-atomic forces should be -41.07531759 eV per primitive unit cell; it is the "binding energy" of the crystal - ie: the energy required to pull all of the ions of the crystal to infinite separation.

The '''Include Shells''' option should ''NOT'' be selected during calculation settings. This is because each ion is considered as a non-polarisable hard sphere here.

= Lattice Vibrations - Computing the Phonons =
Here the free energy will be computed for ever increasing sizes of grid and its convergence to the infinite grid value monitored. The convergence of the density of states plot in the previous exercise is a useful guide for choosing the optimal grid size.

== The relationship between phonon dispersion and DOS ==
[[File:Dispersion x 1x1x1(em2015).png|left|thumb|560x560px|Fig.3&4 The left hand side is the phonon dispersion diagram and the right hand side is the DOS diagram with 1x1x1 grid size.]]
MgO has two ions and in a 3D structure(kx, ky, kz in reciprocal space), overall of 6 branches, as shwon in the phonon dispersion diagram.

The 1x1x1 grid size gives 4 very sharp peaks which indicates two similar higher density of states at around 280 cm-1 and 350 cm-1, and two similar lower density of states at around 670 cm-1 and 810 cm-1. The higher peaks and the lower peaks have a approximate ratio of 2:1 in height.

Furthermore, the DOS diagram can be correlated to the dispersion diagram. The point L in the dispersion diagram intersects six branches, each at 810 cm-1, 670 cm-1, 350 cm-1,350 cm-1,280 cm-1 and 280 cm-1. The two 350 cm-1 frequencies are degenerate at point L and so are the two 280 cm-1 frequencies. This, corresponds to the DOS diagram, is the 2:1 ratio in height.

The L point is 1/2 1/2 1/2.

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

Thus the K point is 1/2 1/2 1/2. If check with the log. file:

'''K point 1 = 0.500000 0.500000 0.500000. '''

This means that the speculation of the relationship between phonon dispersion diagram and DOS diagram is correct.

== Density of States and Grid Size ==
In order to use more '''k''' points, the the density of states calculation is reran, increasing the shrinking factors from 1x1x1, to 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 24x24x24, 36x36x36, 64x64x64. And the density of states diagrams are shown below.
{| class="wikitable"
![[File:DOS2x2x2(em2015).PNG|none|thumb|220x220px|Fig.5 Grid size: 2x2x2]]
![[File:DOS3x3x3(em2015).PNG|none|thumb|223x223px|Fig.6 Grid size: 3x3x3]]
![[File:DOS4x4x4(em2015).PNG|none|thumb|237x237px|Fig.7 Grid size: 4x4x4]]
![[File:DOS8x8x8(em2015).PNG|none|thumb|222x222px|Fig.8 Grid size: 8x8x8]]
|-
|[[File:DOS16x16x16(em2015).PNG|none|thumb|222x222px|'''Fig.9 Grid size: 16x16x16''']]
|[[File:DOS24x24x24(em2015).PNG|none|thumb|221x221px|'''Fig.10 Grid size: 24x24x24''']]
|[[File:DOS32x32x32(em2015).PNG|none|thumb|232x232px|'''Fig.11 Grid size: 36x36x36''']]
|[[File:DOS64x64x64(em2015).PNG|none|thumb|224x224px|'''Fig.12 Grid size: 64x64x64''']]
|}

When increasing the grid size, there are more peaks. When using a 8x8x8 grid size, the peaks connect to each other, giving continuous segments of lines. More and more of the possible vibrations are sampled and more features appear. When higher grid size is applied, the shape becomes more and more smooth and no apparent changes in the shape of the diagram after grid size 24x24x24. This means that the 24x24x24 grid size is the minimum/optimum.

The optimal grid size for MgO is appropriate for a calculation on a similar oxide like calcium oxide. The distance between atoms a is similar for Ca-O and Mg-O.

According to the equation of reciprocal space,

<math>a^{\ast} = \frac{2\Pi}{a}</math>, a* is similar as a is similar.

For Zeolite, which its primitive cell contains more atoms resulting in more branches, and larger lattice size (15.56 Å compared to 2.9783 Å of MgO) resulting in smaller reciprocal space. Thus a smaller grid size is enough to catch all the features from its phonon dispersion diagram. This is just opposite to a metal like lithium. It has only one atom in its primitive cell so only three branches in the phonon dispersion diagram. Lithium also has smaller lattice size, giving a larger a*. To put all factors together, lithium needs larger optimal grid size to include all features.

= Calculating the Free Energy in the Harmonic Approximation =
The free energies of the different grid sizes in last section are reported here.
{| class="wikitable"
|+Table. 1.
!Grid Size
!Free Energy/eV
!Difference in energy/eV
!Accuracy
!Comments
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|<nowiki>-</nowiki>
|
| rowspan="9" |How does the free energy vary with grid size ?
As grid size increases, the free energy increases as well.

|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|<nowiki>-0.003692</nowiki>
|
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|0.000177
|accurate to 1 meV
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|<nowiki>-0.000018</nowiki>
|
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|<nowiki>-0.000028</nowiki>
|accurate to 0.5 meV
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|<nowiki>-0.000004</nowiki>
|accurate to 0.1 meV
|-
|24x24x24
|<nowiki>-40.926483</nowiki>
|<nowiki>-0.000001</nowiki>
|
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0.000000
|
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|0.000000
|}
The optimal grid size speculation is similar to that in the previous section.

= The Thermal Expansion of MgO =
[[File:Energy vs t(em2015).PNG|thumb|592x592px|Fig.13 Graph of free energy against temperature between 0K and 1000K]]
[[File:Lattice vs t(em2015).PNG|thumb|591x591px|Fig.14 Graph of lattice constant against temperature between 0K and 1000K]]
Vary the temperature from 0K to 1000K in steps of 100K and compute the variation in the free energy and optimal lattice constant by optimizing the structure at each temperature. The numbers are extracted from the log. files and organized in the table 2.
{| class="wikitable"
|+Table. 2. The table of free energy and lattice constant at temperatures between 0K and 1000K with a variation of 100K
!Temperature/K
!Free Energy/eV
!Lattice Constant/Å
!Volume/Å3
|-
|0
|<nowiki>-40.90190628</nowiki>
|2.986563
|26.638824
|-
|100
|<nowiki>-40.90241965</nowiki>
|2.986563
|26.638824
|-
|200
|<nowiki>-40.90937738</nowiki>
|2.987605
|26.666716
|-
|300
|<nowiki>-40.92812471</nowiki>
|2.989390
|26.714542
|-
|400
|<nowiki>-40.95859417</nowiki>
|2.991629
|26.774613
|-
|500
|<nowiki>-40.99943595</nowiki>
|2.994135
|26.841955
|-
|600
|<nowiki>-41.04931542</nowiki>
|2.996820
|26.914231
|-
|700
|<nowiki>-41.10711924</nowiki>
|2.999643
|26.990362
|-
|800
|<nowiki>-41.17189187</nowiki>
|3.002588
|27.069936
|-
|900
|<nowiki>-41.24301814</nowiki>
|3.005635
|27.152431
|-
|1000
|<nowiki>-41.31984836</nowiki>
|3.008784
|27.237863
|}
The graph of free energy against temperature is shown in Fig.13, and the plot of lattice constant against temperature is shown in Fig.14. Fig.13 contains a downward curve and its gradient increases as temperature increases. The gradient becomes constant at higher temperatures, which indicated a straight line. Fig. 14 contains an upward curve and its gradient increases as temperature increases. The line becomes straight after 400K.
[[File:HightT(em2015).png|thumb|589x589px|Fig.15 Volume per formula unit at high temperatures]]
The coefficient of thermal expansion of a material is :

<math>\alpha=\frac{1}{Vo}(\frac{dV}{dT})P</math>= 2.6544 x10-5 K-1. (3)

The coefficient of volume thermal expansion is 4.0 × 10−5 K−1[2] in the temperature range from 300 K to 1250 K from A. S. Madhusudhan Rao and K. Narender's paper, which is larger than the calculated value.

For the harmonic approximation, the bond between atoms will eventually break if temperature is high enough, which indicates that the equilibrium distance between atoms is independent of temperature.The quasi-harmonic approximation expands upon the harmonic phonon model of lattice dynamics. In the quasi-harmonic approximation, bonds will never break but lengthen to infinite. In the diagrams below, the center of bonds are indicated in red in each approximations.
{| class="wikitable"
|[[File:Mp-5598-qha(em2015).png|thumb|200x200px|Fig.16 The Quasi-harmonic approximation[3] graph where the center of the bond changes. ]]
|[[File:800px-Morse-potential(em2015).png|thumb|217x217px|Fig.17 The harmonic approximation[4] graph where the center of the bond remains. ]]
|}
When a material is heated, the kinetic energy of that material is increased and it's atoms and molecules move about more. This means that each atom will take up more space due to it's movement so the material will expand.

Fig.15 is the graph of volume per formula unit under QHA at high temperatures. The blue points are the ones used to calculate the thermal expansion coefficient. The It is clear in the graph that the line deviates from the best fit line at high temperatures. Since the QHA assumes no bond-breaking between atoms even at high temperatures, the line should be linear. Thus the approximation is invalid at high temperatures.

= Molecular Dynamics =
[[File:QH-MD(em2015).png|thumb|443x443px|Fig.18 Graph of volume per formula unit under Quasi-harmonic and Molecular Dynamics Approximation against temperature|left]][[File:QHA VS MD(em2015).png|thumb|497x497px|centre|Fig.19 Graph of Volume per formula unit against temperature of QHA and MD at high temperatures]]Fig.18 shows the volume per formula unit under two approximations between 100 K and 1000 K. The coefficient of thermal expansion under molecular dynamics approximation is 3.2024 x10-5 K-1, which is larger than that under quasi-harmonic approximation. However, this number is closer to the measured literature value. From the diagram, it is obvious that the shape of the line of quasi-harmonic is flat at first and becomes linear afterwards. This is because at very low temperatures(0 K to 100 K), the Helmholtz free energy is dominated by zero point energy. And after that, the temperature-dependent part of equation 1 dominates. The points from MD are scattered evenly near the best fit line from the very beginning.

The difference between Fig18 and Fig.19 is the range of temperature. Fig.19 expand the temperature range to 4000 K and the MD line still shows well-scattered points at high temperatures. Therefore, at low temperatures, there is only slight difference between QHA and MD approximation; but, at higher temperatures, the QHA fails and the MD approximation works. The cell volume in the QHA is larger than expected and that in the MD approximation is anticipated using its thermal expansion coefficient. Nevertheless, the melting temperature of MgO is 3125 K, and bonds break after that point. The MD approximation shows that even after melting point of MgO, there is still no bond-breaking process. This indicates that MD approximation is only appropriate below melting point of the material, not beyond.

In order to compute an accurate Free Energy, the best situation is to use a cell size as large as possible, like an infinite cell for MD calculations. This is because in small cells, all atoms move in the same direction and distance when given an energy. This is not the real life situation, where atoms move randomly and freely. A large cell averages out their motion and give a unifying result. However, this is not realistic because the calculation will last infinitely. Therefore, the alternative is to start calculation with small sizes of MgO and then increase the cell size to see how the free energy goes. The number should converge and at one cell size, the free energy won't change much. This size is a good choice to perform reliably MD for MgO.

= Conclusion =

give a general description of your calculations and your main findings

outline the differences between the methods in use and the results obtained

analyse critically these differences

= References =
# Hoffmann, R. (1987), How Chemistry and Physics Meet in the Solid State. Angew. Chem. Int. Ed. Engl., 26: 846–878. doi:10.1002/anie.198708461
# A. S. Madhusudhan Rao and K. Narender, “Studies on Thermophysical Properties of CaO and MgO by γ-Ray Attenuation,” Journal of Thermodynamics, vol. 2014, Article ID 123478, 8 pages, 2014. doi:10.1155/2014/123478
# S. P. Ong, S. Cholia, A. Jain, M. Brafman, D. Gunter, G. Ceder, and K. A. Persson, Computational Materials Science, 2015, 97, 209–215. doi:10.1016/j.commatsci.2014.10.037
# https://chem.libretexts.org/Textbook_Maps/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_(McQuarrie_and_Simon)/05._The_Harmonic_Oscillator_and_the_Rigid_Rotator%3A_Two_Spectroscopic_Models/5.3%3A_The_Harmonic_Oscillator_is_an_Approximation

Rep:MgO(em2015)

2017-11-17T13:52:37Z

Em2015: /* Density of States and Grid Size */

=Introduction=
The aim of the experiment is to investigate an intrinsic property, coefficient of thermal expansion α of a crystal of MgO, using quasi-harmonic and molecular dynamics approximations under RedHat Linux environment. And α is calculated as 2.6544 x10-5 K-1 between 100 K and 1000 K. The interface is DLV, a package for the visualisation of materials, structures and properties. The calculations are completed by a classical simulation program GULP in DLVisualize.

During the investigation of thermal expansion coefficient, the energies and vibrations of MgO crystal are also explored. Vibrations are important in studying thermal properties, phase transitions, electrical properties,etc. Phonon dispersion, which shares a similar concept as band structures[1], is a graph of the relationship between vibration frequency and energy. Vibrational density of states is also related. It shows the relative ratio of states at each vibration frequencies and it's derived from phonon dispersion. The grid size chosen in the calculation has a strong impact on the shape.

The quasi-harmonic approximation(QHA) is a more adapted hypothesis than the harmonic approximation. They have the same limitation that they do not predict bond dissociation. QHA is used to compute vibrational energy levels of MgO, with equation as following: <math>F=E_0+\frac{1}{2}\sum_{k,j}\hslash\omega_{k,j}+k_BT\sum_{k,j}ln[1-exp(-\hslash\omega_{k,j}/k_BT)]</math> (1)

where E0 is the internal energy; the sum of ħω for each k, j is the zero-point energy; and the last term is the temperature-dependent free energy. Zero-point energy is the lowest possible energy that a quantum mechanical system may have.

The molecular dynamics approximation simulates the vibrations as random motions of atoms inside a supercell, using classical Newton's Law:

<math>F=ma</math> (2)

The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic evolution of the system.

The quasi-harmonic approximation and molecular dynamics approximation are compared and contrasted at both low temperatures and high temperatures, to explore which methodology is more appropriate for thermal expansion modelling with MgO.

= Calculating the internal energy of an MgO crystal =
[[File:Primitive cell(em2015).PNG|left|160x160px|thumb|Fig.1 The structure of a primitive unit cell of MgO. ]]
[[File:Conventional cell(em2015).PNG|thumb|179x179px|Fig.2 conventional cell of MgO crystal]]
MgO crystal is a face-centered cube as shown on the right hand side. In the experiment, its crystal structure is described using a primitive unit cell, with cell parameters :a=b=c=2.9783Å ; α=β=γ=60.0000'''°'''. The primitive cell include one magnesium ion (1/8 x 8)and one oxygen ion(1 x 1). Conventional cell has 8 atoms, which indicates that the ratio of volume between conventional cell and primitive cell is 4:1.

The cell vectors of this cell are (in Angstrom):

0.000000 2.105970 2.105970

2.105970 0.000000 2.105970

2.105970 2.105970 0.000000.

The inter-atomic potential is 6.72119980 eV, and the total energy due to the inter-atomic forces should be -41.07531759 eV per primitive unit cell; it is the "binding energy" of the crystal - ie: the energy required to pull all of the ions of the crystal to infinite separation.

The '''Include Shells''' option should ''NOT'' be selected during calculation settings. This is because each ion is considered as a non-polarisable hard sphere here.

= Lattice Vibrations - Computing the Phonons =
Here the free energy will be computed for ever increasing sizes of grid and its convergence to the infinite grid value monitored. The convergence of the density of states plot in the previous exercise is a useful guide for choosing the optimal grid size.

== The relationship between phonon dispersion and DOS ==
[[File:Dispersion x 1x1x1(em2015).png|left|thumb|560x560px|Fig.3&4 The left hand side is the phonon dispersion diagram and the right hand side is the DOS diagram with 1x1x1 grid size.]]
MgO has two ions and in a 3D structure(kx, ky, kz in reciprocal space), overall of 6 branches, as shwon in the phonon dispersion diagram.

The 1x1x1 grid size gives 4 very sharp peaks which indicates two similar higher density of states at around 280 cm-1 and 350 cm-1, and two similar lower density of states at around 670 cm-1 and 810 cm-1. The higher peaks and the lower peaks have a approximate ratio of 2:1 in height.

Furthermore, the DOS diagram can be correlated to the dispersion diagram. The point L in the dispersion diagram intersects six branches, each at 810 cm-1, 670 cm-1, 350 cm-1,350 cm-1,280 cm-1 and 280 cm-1. The two 350 cm-1 frequencies are degenerate at point L and so are the two 280 cm-1 frequencies. This, corresponds to the DOS diagram, is the 2:1 ratio in height.

The L point is 1/2 1/2 1/2.

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

Thus the K point is 1/2 1/2 1/2. If check with the log. file:

'''K point 1 = 0.500000 0.500000 0.500000. '''

This means that the speculation of the relationship between phonon dispersion diagram and DOS diagram is correct.

== Density of States and Grid Size ==
In order to use more '''k''' points, the the density of states calculation is reran, increasing the shrinking factors from 1x1x1, to 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 24x24x24, 36x36x36, 64x64x64. And the density of states diagrams are shown below.
{| class="wikitable"
![[File:DOS2x2x2(em2015).PNG|none|thumb|220x220px|Fig.5 Grid size: 2x2x2]]
![[File:DOS3x3x3(em2015).PNG|none|thumb|223x223px|Fig.6 Grid size: 3x3x3]]
![[File:DOS4x4x4(em2015).PNG|none|thumb|237x237px|Fig.7 Grid size: 4x4x4]]
![[File:DOS8x8x8(em2015).PNG|none|thumb|222x222px|Fig.8 Grid size: 8x8x8]]
|-
|[[File:DOS16x16x16(em2015).PNG|none|thumb|222x222px|'''Fig.9 Grid size: 16x16x16''']]
|[[File:DOS24x24x24(em2015).PNG|none|thumb|221x221px|'''Fig.10 Grid size: 24x24x24''']]
|[[File:DOS32x32x32(em2015).PNG|none|thumb|232x232px|'''Fig.11 Grid size: 36x36x36''']]
|[[File:DOS64x64x64(em2015).PNG|none|thumb|224x224px|'''Fig.12 Grid size: 64x64x64''']]
|}

When increasing the grid size, there are more peaks. When using a 8x8x8 grid size, the peaks connect to each other, giving continuous segments of lines. More and more of the possible vibrations are sampled and more features appear. When higher grid size is applied, the shape becomes more and more smooth and no apparent changes in the shape of the diagram after grid size 24x24x24. This means that the 24x24x24 grid size is the minimum/optimum.

The optimal grid size for MgO is appropriate for a calculation on a similar oxide like calcium oxide. The distance between atoms a is similar for Ca-O and Mg-O.

According to the equation of reciprocal space,

<math>a^{\ast} = \frac{2\Pi}{a}</math>, a* is similar as a is similar.

For Zeolite, which its primitive cell contains more atoms resulting in more branches, and larger lattice size (15.56 Å compared to 2.9783 Å of MgO) resulting in smaller reciprocal space. Thus a smaller grid size is enough to catch all the features from its phonon dispersion diagram. This is just opposite to a metal like lithium. It has only one atom in its primitive cell so only three branches in the phonon dispersion diagram. Lithium also has smaller lattice size, giving a larger a*. To put all factors together, lithium needs larger optimal grid size to include all features.

= Calculating the Free Energy in the Harmonic Approximation =
The free energies of the different grid sizes in last section are reported here.
{| class="wikitable"
|+Table. 1.
!Grid Size
!Free Energy/eV
!Difference in energy/eV
!Accuracy
!Comments
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|<nowiki>-</nowiki>
|
| rowspan="9" |How does the free energy vary with grid size ?
As grid size increases, the free energy increases as well.

|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|<nowiki>-0.003692</nowiki>
|
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|0.000177
|accurate to 1 meV
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|<nowiki>-0.000018</nowiki>
|
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|<nowiki>-0.000028</nowiki>
|accurate to 0.5 meV
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|<nowiki>-0.000004</nowiki>
|accurate to 0.1 meV
|-
|24x24x24
|<nowiki>-40.926483</nowiki>
|<nowiki>-0.000001</nowiki>
|
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0.000000
|
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|0.000000
|}
The optimal grid size speculation is similar to that in the previous section.

= The Thermal Expansion of MgO =
[[File:Energy vs t(em2015).PNG|thumb|592x592px|Fig.13 Graph of free energy against temperature between 0K and 1000K]]
[[File:Lattice vs t(em2015).PNG|thumb|591x591px|Fig.14 Graph of lattice constant against temperature between 0K and 1000K]]
Vary the temperature from 0K to 1000K in steps of 100K and compute the variation in the free energy and optimal lattice constant by optimizing the structure at each temperature. The numbers are extracted from the log. files and organized in the table 2.
{| class="wikitable"
|+Table. 2. The table of free energy and lattice constant at temperatures between 0K and 1000K with a variation of 100K
!Temperature/K
!Free Energy/eV
!Lattice Constant/Å
!Volume/Å3
|-
|0
|<nowiki>-40.90190628</nowiki>
|2.986563
|26.638824
|-
|100
|<nowiki>-40.90241965</nowiki>
|2.986563
|26.638824
|-
|200
|<nowiki>-40.90937738</nowiki>
|2.987605
|26.666716
|-
|300
|<nowiki>-40.92812471</nowiki>
|2.989390
|26.714542
|-
|400
|<nowiki>-40.95859417</nowiki>
|2.991629
|26.774613
|-
|500
|<nowiki>-40.99943595</nowiki>
|2.994135
|26.841955
|-
|600
|<nowiki>-41.04931542</nowiki>
|2.996820
|26.914231
|-
|700
|<nowiki>-41.10711924</nowiki>
|2.999643
|26.990362
|-
|800
|<nowiki>-41.17189187</nowiki>
|3.002588
|27.069936
|-
|900
|<nowiki>-41.24301814</nowiki>
|3.005635
|27.152431
|-
|1000
|<nowiki>-41.31984836</nowiki>
|3.008784
|27.237863
|}
The graph of free energy against temperature is shown in Fig.13, and the plot of lattice constant against temperature is shown in Fig.14. Fig.13 contains a downward curve and its gradient increases as temperature increases. The gradient becomes constant at higher temperatures, which indicated a straight line. Fig. 14 contains an upward curve and its gradient increases as temperature increases. The line becomes straight after 400K.
[[File:HightT(em2015).png|thumb|589x589px|Fig.15 Volume per formula unit at high temperatures]]
The coefficient of thermal expansion of a material is :

<math>\alpha=\frac{1}{Vo}(\frac{dV}{dT})P</math>= 2.6544 x10-5 K-1. (3)

The coefficient of volume thermal expansion is 4.0 × 10−5 K−1[2] in the temperature range from 300 K to 1250 K from A. S. Madhusudhan Rao and K. Narender's paper, which is larger than the calculated value.

For the harmonic approximation, the bond between atoms will eventually break if temperature is high enough, which indicates that the equilibrium distance between atoms is independent of temperature.The quasi-harmonic approximation expands upon the harmonic phonon model of lattice dynamics. In the quasi-harmonic approximation, bonds will never break but lengthen to infinite. In the diagrams below, the center of bonds are indicated in red in each approximations.
{| class="wikitable"
|[[File:Mp-5598-qha(em2015).png|thumb|200x200px|Fig.16 The Quasi-harmonic approximation[3] graph where the center of the bond changes. ]]
|[[File:800px-Morse-potential(em2015).png|thumb|217x217px|Fig.17 The harmonic approximation[4] graph where the center of the bond remains. ]]
|}
When a material is heated, the kinetic energy of that material is increased and it's atoms and molecules move about more. This means that each atom will take up more space due to it's movement so the material will expand.

Fig.15 is the graph of volume per formula unit under QHA at high temperatures. The blue points are the ones used to calculate the thermal expansion coefficient. The It is clear in the graph that the line deviates from the best fit line at high temperatures. Since the QHA assumes no bond-breaking between atoms even at high temperatures, the line should be linear. Thus the approximation is invalid at high temperatures.

= Molecular Dynamics =
[[File:QH-MD(em2015).png|thumb|443x443px|Fig.18 Graph of volume per formula unit under Quasi-harmonic and Molecular Dynamics Approximation against temperature|left]][[File:QHA VS MD(em2015).png|thumb|497x497px|centre|Fig.19 Graph of Volume per formula unit against temperature of QHA and MD at high temperatures]]Fig.18 shows the volume per formula unit under two approximations between 100 K and 1000 K. The coefficient of thermal expansion under molecular dynamics approximation is 3.2024 x10-5 K-1, which is larger than that under quasi-harmonic approximation. However, this number is closer to the measured literature value. From the diagram, it is obvious that the shape of the line of quasi-harmonic is flat at first and becomes linear afterwards. This is because at very low temperatures(0 K to 100 K), the Helmholtz free energy is dominated by zero point energy. And after that, the temperature-dependent part of equation 1 dominates. The points from MD are scattered evenly near the best fit line from the very beginning.

The difference between Fig18 and Fig.19 is the range of temperature. Fig.19 expand the temperature range to 4000 K and the MD line still shows well-scattered points at high temperatures. Therefore, at low temperatures, there is only slight difference between QHA and MD approximation; but, at higher temperatures, the QHA fails and the MD approximation works. The cell volume in the QHA is larger than expected and that in the MD approximation is anticipated using its thermal expansion coefficient. Nevertheless, the melting temperature of MgO is 3125 K, and bonds break after that point. The MD approximation shows that even after melting point of MgO, there is still no bond-breaking process. This indicates that MD approximation is only appropriate below melting point of the material, not beyond.

In order to compute an accurate Free Energy, the best situation is to use a cell size as large as possible, like an infinite cell for MD calculations. This is because in small cells, all atoms move in the same direction and distance when given an energy. This is not the real life situation, where atoms move randomly and freely. A large cell averages out their motion and give a unifying result. However, this is not realistic because the calculation will last infinitely. Therefore, the alternative is to start calculation with small sizes of MgO and then increase the cell size to see how the free energy goes. The number should converge and at one cell size, the free energy won't change much. This size is a good choice to perform reliably MD for MgO.

= Conclusion =

give a general description of your calculations and your main findings

outline the differences between the methods in use and the results obtained

analyse critically these differences

= References =
# Hoffmann, R. (1987), How Chemistry and Physics Meet in the Solid State. Angew. Chem. Int. Ed. Engl., 26: 846–878. doi:10.1002/anie.198708461
# A. S. Madhusudhan Rao and K. Narender, “Studies on Thermophysical Properties of CaO and MgO by γ-Ray Attenuation,” Journal of Thermodynamics, vol. 2014, Article ID 123478, 8 pages, 2014. doi:10.1155/2014/123478
# S. P. Ong, S. Cholia, A. Jain, M. Brafman, D. Gunter, G. Ceder, and K. A. Persson, Computational Materials Science, 2015, 97, 209–215. doi:10.1016/j.commatsci.2014.10.037
# https://chem.libretexts.org/Textbook_Maps/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_(McQuarrie_and_Simon)/05._The_Harmonic_Oscillator_and_the_Rigid_Rotator%3A_Two_Spectroscopic_Models/5.3%3A_The_Harmonic_Oscillator_is_an_Approximation

Rep:MgO(em2015)

2017-11-17T13:51:56Z

Em2015: /* Density of States and Grid Size */

=Introduction=
The aim of the experiment is to investigate an intrinsic property, coefficient of thermal expansion α of a crystal of MgO, using quasi-harmonic and molecular dynamics approximations under RedHat Linux environment. And α is calculated as 2.6544 x10-5 K-1 between 100 K and 1000 K. The interface is DLV, a package for the visualisation of materials, structures and properties. The calculations are completed by a classical simulation program GULP in DLVisualize.

During the investigation of thermal expansion coefficient, the energies and vibrations of MgO crystal are also explored. Vibrations are important in studying thermal properties, phase transitions, electrical properties,etc. Phonon dispersion, which shares a similar concept as band structures[1], is a graph of the relationship between vibration frequency and energy. Vibrational density of states is also related. It shows the relative ratio of states at each vibration frequencies and it's derived from phonon dispersion. The grid size chosen in the calculation has a strong impact on the shape.

The quasi-harmonic approximation(QHA) is a more adapted hypothesis than the harmonic approximation. They have the same limitation that they do not predict bond dissociation. QHA is used to compute vibrational energy levels of MgO, with equation as following: <math>F=E_0+\frac{1}{2}\sum_{k,j}\hslash\omega_{k,j}+k_BT\sum_{k,j}ln[1-exp(-\hslash\omega_{k,j}/k_BT)]</math> (1)

where E0 is the internal energy; the sum of ħω for each k, j is the zero-point energy; and the last term is the temperature-dependent free energy. Zero-point energy is the lowest possible energy that a quantum mechanical system may have.

The molecular dynamics approximation simulates the vibrations as random motions of atoms inside a supercell, using classical Newton's Law:

<math>F=ma</math> (2)

The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic evolution of the system.

The quasi-harmonic approximation and molecular dynamics approximation are compared and contrasted at both low temperatures and high temperatures, to explore which methodology is more appropriate for thermal expansion modelling with MgO.

= Calculating the internal energy of an MgO crystal =
[[File:Primitive cell(em2015).PNG|left|160x160px|thumb|Fig.1 The structure of a primitive unit cell of MgO. ]]
[[File:Conventional cell(em2015).PNG|thumb|179x179px|Fig.2 conventional cell of MgO crystal]]
MgO crystal is a face-centered cube as shown on the right hand side. In the experiment, its crystal structure is described using a primitive unit cell, with cell parameters :a=b=c=2.9783Å ; α=β=γ=60.0000'''°'''. The primitive cell include one magnesium ion (1/8 x 8)and one oxygen ion(1 x 1). Conventional cell has 8 atoms, which indicates that the ratio of volume between conventional cell and primitive cell is 4:1.

The cell vectors of this cell are (in Angstrom):

0.000000 2.105970 2.105970

2.105970 0.000000 2.105970

2.105970 2.105970 0.000000.

The inter-atomic potential is 6.72119980 eV, and the total energy due to the inter-atomic forces should be -41.07531759 eV per primitive unit cell; it is the "binding energy" of the crystal - ie: the energy required to pull all of the ions of the crystal to infinite separation.

The '''Include Shells''' option should ''NOT'' be selected during calculation settings. This is because each ion is considered as a non-polarisable hard sphere here.

= Lattice Vibrations - Computing the Phonons =
Here the free energy will be computed for ever increasing sizes of grid and its convergence to the infinite grid value monitored. The convergence of the density of states plot in the previous exercise is a useful guide for choosing the optimal grid size.

== The relationship between phonon dispersion and DOS ==
[[File:Dispersion x 1x1x1(em2015).png|left|thumb|560x560px|Fig.3&4 The left hand side is the phonon dispersion diagram and the right hand side is the DOS diagram with 1x1x1 grid size.]]
MgO has two ions and in a 3D structure(kx, ky, kz in reciprocal space), overall of 6 branches, as shwon in the phonon dispersion diagram.

The 1x1x1 grid size gives 4 very sharp peaks which indicates two similar higher density of states at around 280 cm-1 and 350 cm-1, and two similar lower density of states at around 670 cm-1 and 810 cm-1. The higher peaks and the lower peaks have a approximate ratio of 2:1 in height.

Furthermore, the DOS diagram can be correlated to the dispersion diagram. The point L in the dispersion diagram intersects six branches, each at 810 cm-1, 670 cm-1, 350 cm-1,350 cm-1,280 cm-1 and 280 cm-1. The two 350 cm-1 frequencies are degenerate at point L and so are the two 280 cm-1 frequencies. This, corresponds to the DOS diagram, is the 2:1 ratio in height.

The L point is 1/2 1/2 1/2.

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

Thus the K point is 1/2 1/2 1/2. If check with the log. file:

'''K point 1 = 0.500000 0.500000 0.500000. '''

This means that the speculation of the relationship between phonon dispersion diagram and DOS diagram is correct.

== Density of States and Grid Size ==
In order to use more '''k''' points, the the density of states calculation is reran, increasing the shrinking factors from 1x1x1, to 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 24x24x24, 36x36x36, 64x64x64. And the density of states diagrams are shown below.
{| class="wikitable"
![[File:DOS2x2x2(em2015).PNG|none|thumb|220x220px|Fig.5 Grid size: 2x2x2]]
![[File:DOS3x3x3(em2015).PNG|none|thumb|223x223px|Fig.6 Grid size: 3x3x3]]
![[File:DOS4x4x4(em2015).PNG|none|thumb|237x237px|Fig.7 Grid size: 4x4x4]]
![[File:DOS8x8x8(em2015).PNG|none|thumb|222x222px|Fig.8 Grid size: 8x8x8]]
|-
|[[File:DOS16x16x16(em2015).PNG|none|thumb|222x222px|'''Fig.9 Grid size: 16x16x16''']]
|[[File:DOS24x24x24(em2015).PNG|none|thumb|221x221px|'''Fig.10 Grid size: 24x24x24''']]
|[[File:DOS32x32x32(em2015).PNG|none|thumb|232x232px|'''Fig.11 Grid size: 36x36x36''']]
|[[File:DOS64x64x64(em2015).PNG|none|thumb|224x224px|'''Fig.12 Grid size: 64x64x64''']]
|}

When increasing the grid size, there are more peaks. When using a 8x8x8 grid size, the peaks connect to each other, giving continuous segments of lines. More and more of the possible vibrations are sampled and more features appear. When higher grid size is applied, the shape becomes more and more smooth and no apparent changes in the shape of the diagram after grid size 24x24x24. This means that the 24x24x24 grid size is the minimum/optimum.

The optimal grid size for MgO is appropriate for a calculation on a similar oxide like calcium oxide. The distance between atoms a is similar for Ca-O and Mg-O. According to the equation of reciprocal space,
<math>a^{\ast} = \frac{2\Pi}{a}</math>,
a* is similar as a is similar. For Zeolite, which its primitive cell contains more atoms resulting in more branches, and larger lattice size (15.56 Å compared to 2.9783 Å of MgO) resulting in smaller reciprocal space. Thus a smaller grid size is enough to catch all the features from its phonon dispersion diagram. This is just opposite to a metal like lithium. It has only one atom in its primitive cell so only three branches in the phonon dispersion diagram. Lithium also has smaller lattice size, giving a larger a*. To put all factors together, lithium needs larger optimal grid size to include all features.

= Calculating the Free Energy in the Harmonic Approximation =
The free energies of the different grid sizes in last section are reported here.
{| class="wikitable"
|+Table. 1.
!Grid Size
!Free Energy/eV
!Difference in energy/eV
!Accuracy
!Comments
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|<nowiki>-</nowiki>
|
| rowspan="9" |How does the free energy vary with grid size ?
As grid size increases, the free energy increases as well.

|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|<nowiki>-0.003692</nowiki>
|
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|0.000177
|accurate to 1 meV
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|<nowiki>-0.000018</nowiki>
|
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|<nowiki>-0.000028</nowiki>
|accurate to 0.5 meV
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|<nowiki>-0.000004</nowiki>
|accurate to 0.1 meV
|-
|24x24x24
|<nowiki>-40.926483</nowiki>
|<nowiki>-0.000001</nowiki>
|
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0.000000
|
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|0.000000
|}
The optimal grid size speculation is similar to that in the previous section.

= The Thermal Expansion of MgO =
[[File:Energy vs t(em2015).PNG|thumb|592x592px|Fig.13 Graph of free energy against temperature between 0K and 1000K]]
[[File:Lattice vs t(em2015).PNG|thumb|591x591px|Fig.14 Graph of lattice constant against temperature between 0K and 1000K]]
Vary the temperature from 0K to 1000K in steps of 100K and compute the variation in the free energy and optimal lattice constant by optimizing the structure at each temperature. The numbers are extracted from the log. files and organized in the table 2.
{| class="wikitable"
|+Table. 2. The table of free energy and lattice constant at temperatures between 0K and 1000K with a variation of 100K
!Temperature/K
!Free Energy/eV
!Lattice Constant/Å
!Volume/Å3
|-
|0
|<nowiki>-40.90190628</nowiki>
|2.986563
|26.638824
|-
|100
|<nowiki>-40.90241965</nowiki>
|2.986563
|26.638824
|-
|200
|<nowiki>-40.90937738</nowiki>
|2.987605
|26.666716
|-
|300
|<nowiki>-40.92812471</nowiki>
|2.989390
|26.714542
|-
|400
|<nowiki>-40.95859417</nowiki>
|2.991629
|26.774613
|-
|500
|<nowiki>-40.99943595</nowiki>
|2.994135
|26.841955
|-
|600
|<nowiki>-41.04931542</nowiki>
|2.996820
|26.914231
|-
|700
|<nowiki>-41.10711924</nowiki>
|2.999643
|26.990362
|-
|800
|<nowiki>-41.17189187</nowiki>
|3.002588
|27.069936
|-
|900
|<nowiki>-41.24301814</nowiki>
|3.005635
|27.152431
|-
|1000
|<nowiki>-41.31984836</nowiki>
|3.008784
|27.237863
|}
The graph of free energy against temperature is shown in Fig.13, and the plot of lattice constant against temperature is shown in Fig.14. Fig.13 contains a downward curve and its gradient increases as temperature increases. The gradient becomes constant at higher temperatures, which indicated a straight line. Fig. 14 contains an upward curve and its gradient increases as temperature increases. The line becomes straight after 400K.
[[File:HightT(em2015).png|thumb|589x589px|Fig.15 Volume per formula unit at high temperatures]]
The coefficient of thermal expansion of a material is :

<math>\alpha=\frac{1}{Vo}(\frac{dV}{dT})P</math>= 2.6544 x10-5 K-1. (3)

The coefficient of volume thermal expansion is 4.0 × 10−5 K−1[2] in the temperature range from 300 K to 1250 K from A. S. Madhusudhan Rao and K. Narender's paper, which is larger than the calculated value.

For the harmonic approximation, the bond between atoms will eventually break if temperature is high enough, which indicates that the equilibrium distance between atoms is independent of temperature.The quasi-harmonic approximation expands upon the harmonic phonon model of lattice dynamics. In the quasi-harmonic approximation, bonds will never break but lengthen to infinite. In the diagrams below, the center of bonds are indicated in red in each approximations.
{| class="wikitable"
|[[File:Mp-5598-qha(em2015).png|thumb|200x200px|Fig.16 The Quasi-harmonic approximation[3] graph where the center of the bond changes. ]]
|[[File:800px-Morse-potential(em2015).png|thumb|217x217px|Fig.17 The harmonic approximation[4] graph where the center of the bond remains. ]]
|}
When a material is heated, the kinetic energy of that material is increased and it's atoms and molecules move about more. This means that each atom will take up more space due to it's movement so the material will expand.

Fig.15 is the graph of volume per formula unit under QHA at high temperatures. The blue points are the ones used to calculate the thermal expansion coefficient. The It is clear in the graph that the line deviates from the best fit line at high temperatures. Since the QHA assumes no bond-breaking between atoms even at high temperatures, the line should be linear. Thus the approximation is invalid at high temperatures.

= Molecular Dynamics =
[[File:QH-MD(em2015).png|thumb|443x443px|Fig.18 Graph of volume per formula unit under Quasi-harmonic and Molecular Dynamics Approximation against temperature|left]][[File:QHA VS MD(em2015).png|thumb|497x497px|centre|Fig.19 Graph of Volume per formula unit against temperature of QHA and MD at high temperatures]]Fig.18 shows the volume per formula unit under two approximations between 100 K and 1000 K. The coefficient of thermal expansion under molecular dynamics approximation is 3.2024 x10-5 K-1, which is larger than that under quasi-harmonic approximation. However, this number is closer to the measured literature value. From the diagram, it is obvious that the shape of the line of quasi-harmonic is flat at first and becomes linear afterwards. This is because at very low temperatures(0 K to 100 K), the Helmholtz free energy is dominated by zero point energy. And after that, the temperature-dependent part of equation 1 dominates. The points from MD are scattered evenly near the best fit line from the very beginning.

The difference between Fig18 and Fig.19 is the range of temperature. Fig.19 expand the temperature range to 4000 K and the MD line still shows well-scattered points at high temperatures. Therefore, at low temperatures, there is only slight difference between QHA and MD approximation; but, at higher temperatures, the QHA fails and the MD approximation works. The cell volume in the QHA is larger than expected and that in the MD approximation is anticipated using its thermal expansion coefficient. Nevertheless, the melting temperature of MgO is 3125 K, and bonds break after that point. The MD approximation shows that even after melting point of MgO, there is still no bond-breaking process. This indicates that MD approximation is only appropriate below melting point of the material, not beyond.

In order to compute an accurate Free Energy, the best situation is to use a cell size as large as possible, like an infinite cell for MD calculations. This is because in small cells, all atoms move in the same direction and distance when given an energy. This is not the real life situation, where atoms move randomly and freely. A large cell averages out their motion and give a unifying result. However, this is not realistic because the calculation will last infinitely. Therefore, the alternative is to start calculation with small sizes of MgO and then increase the cell size to see how the free energy goes. The number should converge and at one cell size, the free energy won't change much. This size is a good choice to perform reliably MD for MgO.

= Conclusion =

give a general description of your calculations and your main findings

outline the differences between the methods in use and the results obtained

analyse critically these differences

= References =
# Hoffmann, R. (1987), How Chemistry and Physics Meet in the Solid State. Angew. Chem. Int. Ed. Engl., 26: 846–878. doi:10.1002/anie.198708461
# A. S. Madhusudhan Rao and K. Narender, “Studies on Thermophysical Properties of CaO and MgO by γ-Ray Attenuation,” Journal of Thermodynamics, vol. 2014, Article ID 123478, 8 pages, 2014. doi:10.1155/2014/123478
# S. P. Ong, S. Cholia, A. Jain, M. Brafman, D. Gunter, G. Ceder, and K. A. Persson, Computational Materials Science, 2015, 97, 209–215. doi:10.1016/j.commatsci.2014.10.037
# https://chem.libretexts.org/Textbook_Maps/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_(McQuarrie_and_Simon)/05._The_Harmonic_Oscillator_and_the_Rigid_Rotator%3A_Two_Spectroscopic_Models/5.3%3A_The_Harmonic_Oscillator_is_an_Approximation

Rep:MgO(em2015)

2017-11-17T13:42:55Z

Em2015: /* Introduction */

=Introduction=
The aim of the experiment is to investigate an intrinsic property, coefficient of thermal expansion α of a crystal of MgO, using quasi-harmonic and molecular dynamics approximations under RedHat Linux environment. And α is calculated as 2.6544 x10-5 K-1 between 100 K and 1000 K. The interface is DLV, a package for the visualisation of materials, structures and properties. The calculations are completed by a classical simulation program GULP in DLVisualize.

During the investigation of thermal expansion coefficient, the energies and vibrations of MgO crystal are also explored. Vibrations are important in studying thermal properties, phase transitions, electrical properties,etc. Phonon dispersion, which shares a similar concept as band structures[1], is a graph of the relationship between vibration frequency and energy. Vibrational density of states is also related. It shows the relative ratio of states at each vibration frequencies and it's derived from phonon dispersion. The grid size chosen in the calculation has a strong impact on the shape.

The quasi-harmonic approximation(QHA) is a more adapted hypothesis than the harmonic approximation. They have the same limitation that they do not predict bond dissociation. QHA is used to compute vibrational energy levels of MgO, with equation as following: <math>F=E_0+\frac{1}{2}\sum_{k,j}\hslash\omega_{k,j}+k_BT\sum_{k,j}ln[1-exp(-\hslash\omega_{k,j}/k_BT)]</math> (1)

where E0 is the internal energy; the sum of ħω for each k, j is the zero-point energy; and the last term is the temperature-dependent free energy. Zero-point energy is the lowest possible energy that a quantum mechanical system may have.

The molecular dynamics approximation simulates the vibrations as random motions of atoms inside a supercell, using classical Newton's Law:

<math>F=ma</math> (2)

The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic evolution of the system.

The quasi-harmonic approximation and molecular dynamics approximation are compared and contrasted at both low temperatures and high temperatures, to explore which methodology is more appropriate for thermal expansion modelling with MgO.

= Calculating the internal energy of an MgO crystal =
[[File:Primitive cell(em2015).PNG|left|160x160px|thumb|Fig.1 The structure of a primitive unit cell of MgO. ]]
[[File:Conventional cell(em2015).PNG|thumb|179x179px|Fig.2 conventional cell of MgO crystal]]
MgO crystal is a face-centered cube as shown on the right hand side. In the experiment, its crystal structure is described using a primitive unit cell, with cell parameters :a=b=c=2.9783Å ; α=β=γ=60.0000'''°'''. The primitive cell include one magnesium ion (1/8 x 8)and one oxygen ion(1 x 1). Conventional cell has 8 atoms, which indicates that the ratio of volume between conventional cell and primitive cell is 4:1.

The cell vectors of this cell are (in Angstrom):

0.000000 2.105970 2.105970

2.105970 0.000000 2.105970

2.105970 2.105970 0.000000.

The inter-atomic potential is 6.72119980 eV, and the total energy due to the inter-atomic forces should be -41.07531759 eV per primitive unit cell; it is the "binding energy" of the crystal - ie: the energy required to pull all of the ions of the crystal to infinite separation.

The '''Include Shells''' option should ''NOT'' be selected during calculation settings. This is because each ion is considered as a non-polarisable hard sphere here.

= Lattice Vibrations - Computing the Phonons =
Here the free energy will be computed for ever increasing sizes of grid and its convergence to the infinite grid value monitored. The convergence of the density of states plot in the previous exercise is a useful guide for choosing the optimal grid size.

== The relationship between phonon dispersion and DOS ==
[[File:Dispersion x 1x1x1(em2015).png|left|thumb|560x560px|Fig.3&4 The left hand side is the phonon dispersion diagram and the right hand side is the DOS diagram with 1x1x1 grid size.]]
MgO has two ions and in a 3D structure(kx, ky, kz in reciprocal space), overall of 6 branches, as shwon in the phonon dispersion diagram.

The 1x1x1 grid size gives 4 very sharp peaks which indicates two similar higher density of states at around 280 cm-1 and 350 cm-1, and two similar lower density of states at around 670 cm-1 and 810 cm-1. The higher peaks and the lower peaks have a approximate ratio of 2:1 in height.

Furthermore, the DOS diagram can be correlated to the dispersion diagram. The point L in the dispersion diagram intersects six branches, each at 810 cm-1, 670 cm-1, 350 cm-1,350 cm-1,280 cm-1 and 280 cm-1. The two 350 cm-1 frequencies are degenerate at point L and so are the two 280 cm-1 frequencies. This, corresponds to the DOS diagram, is the 2:1 ratio in height.

The L point is 1/2 1/2 1/2.

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

Thus the K point is 1/2 1/2 1/2. If check with the log. file:

'''K point 1 = 0.500000 0.500000 0.500000. '''

This means that the speculation of the relationship between phonon dispersion diagram and DOS diagram is correct.

== Density of States and Grid Size ==
In order to use more '''k''' points, the the density of states calculation is reran, increasing the shrinking factors from 1x1x1, to 2x2x2, 3x3x3, 4x4x4, 8x8x8, 16x16x16, 24x24x24, 36x36x36, 64x64x64. And the density of states diagrams are shown below.
{| class="wikitable"
![[File:DOS2x2x2(em2015).PNG|none|thumb|220x220px|Fig.5 Grid size: 2x2x2]]
![[File:DOS3x3x3(em2015).PNG|none|thumb|223x223px|Fig.6 Grid size: 3x3x3]]
![[File:DOS4x4x4(em2015).PNG|none|thumb|237x237px|Fig.7 Grid size: 4x4x4]]
![[File:DOS8x8x8(em2015).PNG|none|thumb|222x222px|Fig.8 Grid size: 8x8x8]]
|-
|[[File:DOS16x16x16(em2015).PNG|none|thumb|222x222px|'''Fig.9 Grid size: 16x16x16''']]
|[[File:DOS24x24x24(em2015).PNG|none|thumb|221x221px|'''Fig.10 Grid size: 24x24x24''']]
|[[File:DOS32x32x32(em2015).PNG|none|thumb|232x232px|'''Fig.11 Grid size: 36x36x36''']]
|[[File:DOS64x64x64(em2015).PNG|none|thumb|224x224px|'''Fig.12 Grid size: 64x64x64''']]
|}

When increasing the grid size, there are more peaks. When using a 8x8x8 grid size, the peaks connect to each other, giving continuous segments of lines. More and more of the possible vibrations are sampled and more features appear. When higher grid size is applied, the shape becomes more and more smooth and no apparent changes in the shape of the diagram after grid size 24x24x24. This means that the 24x24x24 grid size is the minimum/optimum.

The optimal grid size for MgO is appropriate for a calculation on a similar oxide like calcium oxide. The distance between atoms a is similar for Ca-O and Mg-O. According to the equation of reciprocal space,
<math>a*{*}=\frac{2π}{a}</math>,
a* is similar as a is similar. For Zeolite, which its primitive cell contains more atoms resulting in more branches, and larger lattice size (15.56 Å compared to 2.9783 Å of MgO) resulting in smaller reciprocal space. Thus a smaller grid size is enough to catch all the features from its phonon dispersion diagram. This is just opposite to a metal like lithium. It has only one atom in its primitive cell so only three branches in the phonon dispersion diagram. Lithium also has smaller lattice size, giving a larger a*. To put all factors together, lithium needs larger optimal grid size to include all features.

= Calculating the Free Energy in the Harmonic Approximation =
The free energies of the different grid sizes in last section are reported here.
{| class="wikitable"
|+Table. 1.
!Grid Size
!Free Energy/eV
!Difference in energy/eV
!Accuracy
!Comments
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|<nowiki>-</nowiki>
|
| rowspan="9" |How does the free energy vary with grid size ?
As grid size increases, the free energy increases as well.

|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|<nowiki>-0.003692</nowiki>
|
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|0.000177
|accurate to 1 meV
|-
|4x4x4
|<nowiki>-40.926450</nowiki>
|<nowiki>-0.000018</nowiki>
|
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|<nowiki>-0.000028</nowiki>
|accurate to 0.5 meV
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|<nowiki>-0.000004</nowiki>
|accurate to 0.1 meV
|-
|24x24x24
|<nowiki>-40.926483</nowiki>
|<nowiki>-0.000001</nowiki>
|
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|0.000000
|
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|0.000000
|}
The optimal grid size speculation is similar to that in the previous section.

= The Thermal Expansion of MgO =
[[File:Energy vs t(em2015).PNG|thumb|592x592px|Fig.13 Graph of free energy against temperature between 0K and 1000K]]
[[File:Lattice vs t(em2015).PNG|thumb|591x591px|Fig.14 Graph of lattice constant against temperature between 0K and 1000K]]
Vary the temperature from 0K to 1000K in steps of 100K and compute the variation in the free energy and optimal lattice constant by optimizing the structure at each temperature. The numbers are extracted from the log. files and organized in the table 2.
{| class="wikitable"
|+Table. 2. The table of free energy and lattice constant at temperatures between 0K and 1000K with a variation of 100K
!Temperature/K
!Free Energy/eV
!Lattice Constant/Å
!Volume/Å3
|-
|0
|<nowiki>-40.90190628</nowiki>
|2.986563
|26.638824
|-
|100
|<nowiki>-40.90241965</nowiki>
|2.986563
|26.638824
|-
|200
|<nowiki>-40.90937738</nowiki>
|2.987605
|26.666716
|-
|300
|<nowiki>-40.92812471</nowiki>
|2.989390
|26.714542
|-
|400
|<nowiki>-40.95859417</nowiki>
|2.991629
|26.774613
|-
|500
|<nowiki>-40.99943595</nowiki>
|2.994135
|26.841955
|-
|600
|<nowiki>-41.04931542</nowiki>
|2.996820
|26.914231
|-
|700
|<nowiki>-41.10711924</nowiki>
|2.999643
|26.990362
|-
|800
|<nowiki>-41.17189187</nowiki>
|3.002588
|27.069936
|-
|900
|<nowiki>-41.24301814</nowiki>
|3.005635
|27.152431
|-
|1000
|<nowiki>-41.31984836</nowiki>
|3.008784
|27.237863
|}
The graph of free energy against temperature is shown in Fig.13, and the plot of lattice constant against temperature is shown in Fig.14. Fig.13 contains a downward curve and its gradient increases as temperature increases. The gradient becomes constant at higher temperatures, which indicated a straight line. Fig. 14 contains an upward curve and its gradient increases as temperature increases. The line becomes straight after 400K.
[[File:HightT(em2015).png|thumb|589x589px|Fig.15 Volume per formula unit at high temperatures]]
The coefficient of thermal expansion of a material is :

<math>\alpha=\frac{1}{Vo}(\frac{dV}{dT})P</math>= 2.6544 x10-5 K-1. (3)

The coefficient of volume thermal expansion is 4.0 × 10−5 K−1[2] in the temperature range from 300 K to 1250 K from A. S. Madhusudhan Rao and K. Narender's paper, which is larger than the calculated value.

For the harmonic approximation, the bond between atoms will eventually break if temperature is high enough, which indicates that the equilibrium distance between atoms is independent of temperature.The quasi-harmonic approximation expands upon the harmonic phonon model of lattice dynamics. In the quasi-harmonic approximation, bonds will never break but lengthen to infinite. In the diagrams below, the center of bonds are indicated in red in each approximations.
{| class="wikitable"
|[[File:Mp-5598-qha(em2015).png|thumb|200x200px|Fig.16 The Quasi-harmonic approximation[3] graph where the center of the bond changes. ]]
|[[File:800px-Morse-potential(em2015).png|thumb|217x217px|Fig.17 The harmonic approximation[4] graph where the center of the bond remains. ]]
|}
When a material is heated, the kinetic energy of that material is increased and it's atoms and molecules move about more. This means that each atom will take up more space due to it's movement so the material will expand.

Fig.15 is the graph of volume per formula unit under QHA at high temperatures. The blue points are the ones used to calculate the thermal expansion coefficient. The It is clear in the graph that the line deviates from the best fit line at high temperatures. Since the QHA assumes no bond-breaking between atoms even at high temperatures, the line should be linear. Thus the approximation is invalid at high temperatures.

= Molecular Dynamics =
[[File:QH-MD(em2015).png|thumb|443x443px|Fig.18 Graph of volume per formula unit under Quasi-harmonic and Molecular Dynamics Approximation against temperature|left]][[File:QHA VS MD(em2015).png|thumb|497x497px|centre|Fig.19 Graph of Volume per formula unit against temperature of QHA and MD at high temperatures]]Fig.18 shows the volume per formula unit under two approximations between 100 K and 1000 K. The coefficient of thermal expansion under molecular dynamics approximation is 3.2024 x10-5 K-1, which is larger than that under quasi-harmonic approximation. However, this number is closer to the measured literature value. From the diagram, it is obvious that the shape of the line of quasi-harmonic is flat at first and becomes linear afterwards. This is because at very low temperatures(0 K to 100 K), the Helmholtz free energy is dominated by zero point energy. And after that, the temperature-dependent part of equation 1 dominates. The points from MD are scattered evenly near the best fit line from the very beginning.

The difference between Fig18 and Fig.19 is the range of temperature. Fig.19 expand the temperature range to 4000 K and the MD line still shows well-scattered points at high temperatures. Therefore, at low temperatures, there is only slight difference between QHA and MD approximation; but, at higher temperatures, the QHA fails and the MD approximation works. The cell volume in the QHA is larger than expected and that in the MD approximation is anticipated using its thermal expansion coefficient. Nevertheless, the melting temperature of MgO is 3125 K, and bonds break after that point. The MD approximation shows that even after melting point of MgO, there is still no bond-breaking process. This indicates that MD approximation is only appropriate below melting point of the material, not beyond.

In order to compute an accurate Free Energy, the best situation is to use a cell size as large as possible, like an infinite cell for MD calculations. This is because in small cells, all atoms move in the same direction and distance when given an energy. This is not the real life situation, where atoms move randomly and freely. A large cell averages out their motion and give a unifying result. However, this is not realistic because the calculation will last infinitely. Therefore, the alternative is to start calculation with small sizes of MgO and then increase the cell size to see how the free energy goes. The number should converge and at one cell size, the free energy won't change much. This size is a good choice to perform reliably MD for MgO.

= Conclusion =

give a general description of your calculations and your main findings

outline the differences between the methods in use and the results obtained

analyse critically these differences

= References =
# Hoffmann, R. (1987), How Chemistry and Physics Meet in the Solid State. Angew. Chem. Int. Ed. Engl., 26: 846–878. doi:10.1002/anie.198708461
# A. S. Madhusudhan Rao and K. Narender, “Studies on Thermophysical Properties of CaO and MgO by γ-Ray Attenuation,” Journal of Thermodynamics, vol. 2014, Article ID 123478, 8 pages, 2014. doi:10.1155/2014/123478
# S. P. Ong, S. Cholia, A. Jain, M. Brafman, D. Gunter, G. Ceder, and K. A. Persson, Computational Materials Science, 2015, 97, 209–215. doi:10.1016/j.commatsci.2014.10.037
# https://chem.libretexts.org/Textbook_Maps/Physical_and_Theoretical_Chemistry_Textbook_Maps/Map%3A_Physical_Chemistry_(McQuarrie_and_Simon)/05._The_Harmonic_Oscillator_and_the_Rigid_Rotator%3A_Two_Spectroscopic_Models/5.3%3A_The_Harmonic_Oscillator_is_an_Approximation