View source for Rep:MgO zy3915

= '''Thermal expansion of MgO''' =

==== - Third Year Computational Laboratory  ====

==== Zhuohao You ====

==== In this experiment, two different model of a Face-centered cubic crystal of Magnesium oxide to simulate the motions of atoms and then calculate how the volume of MgO crystal changes under different temperature and the related constant, the thermal expansion coefficient, α, with comparison between the two models. ====

== 1-Introduction: ==
[[File:Bond Lengths and Energies.png|frameless|500px]]<ref><nowiki>https://chem.libretexts.org/Core/Physical_and_Theoretical_Chemistry/Chemical_Bonding/Fundamentals_of_Chemical_Bonding/Chemical_Bonds/Bond_Lengths_and_Energies</nowiki></ref>

As is known to all, atoms, even those in their solid states, are in constant vibrations. and this vibration is related to the thermal properties, for example, heat capacity and thermal expansion. Figure [1] shown above illustrates how the bond length varies with the potential energy of H<sub>2</sub> molecule, as temperature (energy) increases, the molecule vibrates between a larger boundary with the extension of the bond rise significantly faster than the compression of the bond. This asymmetry in vibration profile gives rise to the thermal expansion of MgO solid as the same theory applies to the solid, the average bond length increases with temperature. 
== ==
In this experiment, two method,QUASI-HARMONIC APPROXIMATION (LD) and MOLECULAR DYNAMICS (MD), was used to model the vibration of  magnesium oxide crystal under different temperature and therefore to calculate the volume change due to temperature change, which is the coefficient of thermal expansion, using the following formula:

[[File:Formula.png|thumb]|frameless|200px]]

The thermal expansion property is widely used in many fields, especially in material science, in which the volume change of substance is well concerned.  Although in general, the constant varies with temperatures, in this experiment we can still assume a linear relationship between temperature and volume in a certain range of temperature to find out the coefficient using the function linear-fitting in python. All the simulation was done by GULP calculation with DLVisualize in Linux, DLVisualize is a software specialized for quantum mechanical and empirical simulation and the powerful code code GULP (General Utility Lattice Program) was used here for both two model.

===== QUASI-HARMONIC APPROXIMATION (LD) =====
[[File:Shmzy666.png|frameless|350px]]<ref><nowiki>http://www.crystal.unito.it/mssc2016/Lectures/QHA.pdf</nowiki></ref>

As shown in figure [2], the harmonic approximation does not allow the shift in coordination for the atoms, which means there is no thermal expansion in this model. To fit the model that can be used to study volume change of solid, the adjusted model quasi-harmonic approximation was adopted to investigate how the volume of MgO crystal changes with temperature. The quasi-harmonic approximation takes into account that the dependence of phonon frequencies on volume and also allows to evaluate thermal expansion.

===== MOLECULAR DYNAMICS (MD) =====
This is a classical model that applies Newton's Second Law: '''F = ma '''to compute the random motions of atoms  in the crystal and then to obtain the average volume at equilibrium.

===== MgO crystal  =====
Three types of crystal cell was used in the experiment for different purposes：
# CONVENTIONAL CELL: This is a simple and clear cell due to it's cubic shape, it shows the bond relationship between atoms and the lattice constant of conventional cell is also important as it allows a easy calculation from lattice parameter a<sub>c</sub> to the volume of the cell, which is just a<sub>c</sub><sup>3</sup>, ''vice versa''.
# PRIMITIVE CELL: This is the simplest representation of the MgO crystal with a shape of rhombohedron. It was used in the quasi-harmonic approximation as the unit cell for repeating. Also the potential energy was calculated using this unit cell.
# SUPERCELL: This is just a combination of several conventional cells to built larger cubic cell that makes simulation of complex vibration modes possible in Molecular Dynamics method.

== 2-Experiment and Results ==

==== 2.1-An Initial Calculation on MgO: ====
In the first part a simple calculation of potential energy was carried using the primitive cell model to help familiarise the interface.  Interpretation of the output file ( LogFile) was showed to help extracting information, for instance lattice constants, form the data calculated by GULP. Using the atom ratio (i.e.  volume ratio) between primitive cell and conventional cell, which is 8:2, one can easily calculated the lattice parameter a<sub>c</sub>  with the equation V = a<sub>c</sub><sup>3</sup>,as shown above. And a 'binding energy ' of exact -41.07531759 eV (same as the script value) was given by GULP to confirm the calculation process was carried correctly.

==== 2.2-Calculating the phonon modes of MgO: ====
As introduced above, the quasi-harmonic approximation allows phonon frequency calculated with the consideration of thermal expansion. In this part the phonon dispersion curves was computed at 50 points of K-vectors as shown below:

[[File:Phonon disp1zyzy.gif|frameless|500px]]

It was 3-D calculation and because Mg atoms are non-identical with O atoms, 6 path was observed in the phonon diagram. Every single phonon mode on the curves could be animated by DLVisualize to give an intuitive ideal about how the crystal vibrates.

And then The phonon density of states was counted across the phonon dispersion curves:

[[File:Phonon dos1zyzy.gif|frameless|500px]]

The figure above shows the simplest situation when the grid size is 1*1*1. In this case only four phonon mode was covered with a ratio of 2:2:1:1. By checking the output Logfile, further datails was shown by the computed data:

<span lang="EN-GB">''Number
of k points for this configuration =    1''</span>

<span lang="EN-GB">''--------------------------------------------------------------------------------''</span>

<span lang="EN-GB">''  K point  
1 =   0.500000  0.500000 
0.500000  Weight =    1.000''</span>

<span lang="EN-GB">''--------------------------------------------------------------------------------''</span>

<span lang="EN-GB">''  Frequencies (cm-1) :''</span>

<span lang="EN-GB">''  288.49 
288.49  351.76  351.76 
676.23  818.82''</span>

It is clear to see that the K point for this configuration is the L point in the phonon dispersion curves, <span lang="EN-GB"> L: (1/2,1/2,1/2),    [π/''a'',π/''a'',π/''a''] , 3m</span><ref><span lang="EN-GB"><nowiki>http://lampx.tugraz.at/~hadley/ss1/bzones/fcc.php</nowiki></span>
</ref>

6 phonon modes was collected by a single K vector as the line went through 6 path for this crystal, and two pairs of phonon mode which were the same gave rise to the 2:2:1:1 ratio in the phonon density of states.

[[File:2^3zyzy.png|frameless|300px]][[File:Zyzy4^3.png|frameless|300px]][[File:Zyzy8^3.png|frameless|300px]][[File:Zyzy16^3.png|frameless|300px]]

Further calculation of 2*2*2; 4*4*4; 8*8*8; 16*16*16 was carried to give more and more smooth diagram with generally similar profile, more and more phonon modes was covered as more and more K-points was counted.

[[File:Zyzy24^3.png|frameless|450px]][[File:Zyzy32^3.png|frameless|450px]]

The shape of the DOS became quite stable when the grid size went above 24*24*24, as very slight different was observed between 24*24*24 and 32*32*32, which made 24*24*24 to be the minimum grid size for a reasonable approximation to the density of states.( a further 64*64*64 was carried taking about 10 mins of time but only gave the same profile) The calculation time required for the scale of 24*24*24 was less than 1 second, which was very efficient for this kind of simulation. DOS was  basically sampling along the  phonon dispersion path and it seemed that 24 sample along each axis was sufficient enough for this MgO crystal.

==== 2.3-Computing the Free Energy - The harmonic approximation: ====
{| class="wikitable"
!Free energy calculation for different K scale
|-
!
{|
  |Scale
  |Free energy (kJmol-1)
  |eV
  |Δ(meV)
 
 |-
  |1
  |<nowiki>-3949.148006</nowiki>
  |<nowiki>-40.930301</nowiki>
  |
 
 |-
  |2
  |<nowiki>-3948.79181</nowiki>
  |<nowiki>-40.926609</nowiki>
  |3.692
 
 |-
  |4
  |<nowiki>-3948.776419</nowiki>
  |<nowiki>-40.92645</nowiki>
  |0.159
 
 |-
  |8
  |<nowiki>-3948.779123</nowiki>
  |<nowiki>-40.926478</nowiki>
  |<nowiki>-0.028</nowiki>
 
 |-
  |16
  |<nowiki>-3948.77958</nowiki>
  |<nowiki>-40.926482</nowiki>
  |<nowiki>-0.004</nowiki>
 
 |-
  |24
  |<nowiki>-3948.779629</nowiki>
  |<nowiki>-40.926483</nowiki>
  |<nowiki>-0.001</nowiki>
 
 |-
  |32
  |<nowiki>-3948.779641</nowiki>
  |<nowiki>-40.926483</nowiki>
  |0
 |}
|}
The appropriate grid size was confirmed by the free energy calculated at 300K which are consistent with the smoothness of the DOS diagrams. As shown in the table above, the energy difference for each step down the table drop considerably fast and for the last step there was no difference to 6 decimal place. According to the table, a 4*4*4 grid size is good enough for an accuracy of 1meV or even 0.5 meV, and a 8*8*8 size is accurate to 0.1 meV. In the later calculation the size 24*24*24 was used as it gives an rather accurate result with reasonable time consumption.

For the crystal other than magnesium oxide, the symmetry of the lattice need to be take in to account when estimating the appropriate grid size:
# for a similar oxide (eg: CaO)  similar scale required as symmetry is the same, actual size needed might be slightly larger due to a larger size of Ca atoms.
# for a Zeolite (eg: Faujasite) large grid size need due to complex structure and therefore complex phonon modes requiring higher resolution of sampling. 
# for a metal (eg: lithium) small size will be expected to be sufficient as lithium atoms are small and only one type of atoms in the lattice therefore hight symmetric structure, very few of wave function combinations available.  

==== The Thermal Expansion of MgO ( LD & MD ): ====

[[File:MgOzyzy3.png|frameless|450px]]

In the last part, the thermal expansion of MgO was computed in two different modes and compared, and the data collected was plotted on the graph above. 

In general, both model gave a quite linear trend with quasi-harmonic approximation having a 'tail up' at very low temperature close to 0 K, this could be explained by the zero point energy from QM, which it a more realistically situation.  The gradient, which was the coefficient of thermal expansion studied in this experiment, was calculated to be 3.437*10<sup>-5</sup>K<sup>-1</sup> (LD) and

3.271*10<sup>-5</sup>K<sup>-1</sup> (MD) respectively ( uncertainty both about 0.1*10<sup>-5</sup>K<sup>-1</sup> ), the coefficient calculated was adjusted by a factor of intial volume at 300K(18.890024 A^3). And there 0-intercept was rather close to each other: 18.64 and 18.62 A^3 which means the two model met at 0 K and went further apart at higher temperature, if the zero-point energy did not exist.

[[File:Plozyzy.png|frameless|750px]]
* Plots of the free energy against temperature and the lattice constant against temperature
from the plots we can see that the lattice parameter give the same trend as the volume change while the free energy was increased with the temperature, a deep step went down the free energy from 0 to 100K as 0 K was considered to be no motion therefore very small vibrational energy with only the zero point energy.

comparing to the literature data 3.12-5.33*10<sup>-5</sup>K<sup>-1</sup> (300K-2000K<ref><nowiki>https://www.nist.gov/sites/default/files/documents/srd/jpcrd375.pdf</nowiki></ref>), our result was in the range but quite small considering the range of temperature we took, and that was because the linear fitting in python gave the same weight to each point (apart from the first few points in LD that tailed up), therefore the resulting valued was dragged down by those initial points. But, the quasi-harmonic model clearly illustrate the real situation better in terms of the curve graph with the zero-energy, and also a higher constant from linear fitting. The main approximations in the calculation was that the atoms in the solid act as anharmonic oscillators, which is the origin of thermal expansion: the bond stenches more significantly than it compress, leading to a increase in average bond length and therefore increased volume.    

However, this approximation can break at high temperature close to the melting point (about 3000K for MgO). the experiment was further carried at temperature above 2000K, and the linear trend seemed to be maintained by both model at up to 2600K, but a sharp increase was observed for LD at 2900K.      

== Conclusion ==
Overall, it seemed that the quasi-harmonic approximation is a better model compared to classical MD method as more factor (e.g. quantum effect was taken into account), and the data collected showed that the LD simulation tends to give more stable result (little fluctuation) because it was simply summing up states of vibrations, compare to a significant fluctuation for MD especially at high temperatures.