Rep:MgOvl915

2018-03-23T15:58:55Z

Vl915:

== Abstract ==
In this investigation the General Utility Lattice Program (GULP) was used to model the thermal expansion of an MgO crystal by simulating its vibrations at different temperatures. This task was approached in two different ways using methods based on the quasi harmonic approximation and molecular dynamics, respectively. The optimal k grid size for the simulations was found out to be 32x32x32 and temperature-dependent expressions for the thermal volumetric expansion coefficients αV(T) were obtained for both approaches from third order polynomial regressions over a temperature range from 300 K to 2000 K. Comparison of the computed thermal expansion coefficients with experimentally determined expansion coefficients shows that throughout a temperature range from 300 K to 2000 K the quasi-harmonic approximation approach better models the thermal expansion coefficients than the molecular dynamics apporach.

==Introduction==

=== The importance of knowing the thermal expansion coefficient of a material===
Thermal expansion is the tendency of matter to cange in shape, area and volume as a consequence of a change in temperature <ref>P.Tipler, G.Mosca, ''Physics for Scientists and Engineers-Volume 1, ''Worth Publishers, New York, 2008</ref> and is directly related to the atomic spacing within the material. Thermal expansion is quantified by thermal expansion coefficients. The ability to estimate how a certain material will behave as temperature is altered is crucial for successful engineering of bridges, buildings and railways for instance and is also important to consider when constructing electronic devices where a mismatch of the thermal expansion coefficients of the different materials within the device can cause electronic malfunction. Hence, the ability to quantify thermal expansion is very important to guarantee the safety and reliabilty of constructions and devices humanity uses in everyday life. Failure to account for thermal exansion can have detrimental consequences.

There are various methods to determine thermal expansion coefficients experimentally, including X-ray diffraction and high-resolution microscopy,<ref>J.D. James, J.A. Spittle, S.G.R Brown, ''Meas.Sci.Technol.,'' 2001, '''12''', 1-15</ref> meanwhile these expensive experiments can be circumvented by modeling the thermal expansion of a material with computer simulations.

In this investigation the thermal expansion coefficients of a MgO crystal are obtained by modeling its lattice vibrations at different temperatures using two different apporaches, the quasi-harmonic approximation and the MDS method.

===The partition function of a vibrating crystal ===
A partition function of a system (shown in eq. 1) is the sum of Boltzmann factors for every state of that system. They allow to generate a link between properties at the molecular level and properties at the macroscopic level. Hence, once the partition function of a system is known, numerous thermodynamic properties can be deduced from it.
[[File:PARTITION FUNCTION EQVL915.PNG|400px|thumb|center|Equation 1]]

Meanwhile, it is generally diffcult to obtain a partition function of a system as all its states are not always known. An exception are vibrating crystals and it is very desirable to obtain a partition function for crystal vibrations as it allows to compute numerous important properties of a crystal, such as its free energy or its thermal expansion coefficients at different temperatures for example, which will be shown in this investigation.<ref>N. Harrison, ''Vibrations in crystals'' lecture notes</ref>

===Modeling crystal vibrations with the harmonic / quasi harmonic approximation and the Brioullin Zone===
[[File:1Dmodelvl915.PNG|400px|thumb|center|Fig.1: 1D model of a crystal structure]]
To get a better understanding of vibartions in crystals it is important to take a closer look at the model shown in figure 1 which represents a 1D crystal structure as a linear chain of atoms, connected by springs which represent the interatomic potential. At equilibrium the distance between the atoms is equal to the unit cell length 'a'. Meanwhile during a vibration of the crystal the atoms will be displaced from their equilibrium position which will cause a rise in the potential energy of the crystal. The potential of the crystal in figure 1, as a function of the displacement of the atoms from the equilibrium position, can be approximated by a Taylor expansion as shown in equation 2: <ref name=":0">M. Dove, ''Introduction to Lattice Dynamics, ''Cambridge University Press, Cambridge, 1993</ref>

[[File:Taylorvl915.PNG|400px|thumb|center|
Equation 2
]]

In equation 2, σ(a) is the interatomic interaction potential (which can be modeled as a Lennard Jones - or Buckingham potential for instance) evaluated at the equilibrium interatomic spacing 'a', N is the number of atoms in the chain and un-un+1 is the displacement from equilibrium of the nth atom. If we omit all terms beyond s=2 in this Taylor expansion, the potential of the crystal will be effectively a set of harmonic oscillators and this approximation is called the harmonic approximation. The harmonic approximation is mathematically very convenient as the harmonic equations of motions have exact solutions. The solutions of the harmonic equations of motion are sinusodial waves. Thus the motion of the whole 1D model (fig.1) can be described by a set of travelling waves with wavevectors k and corresponding angular frequencies ωk. Using this result and introducing it into equation 3 along with the expression for the harmonic potential, while applying the Born-von Kármán periodic boundary condition, results in an expression for ωk as a function of k (equation 4). <ref name=":0" /> [[File:Eaution of motion vl915.PNG|400px|thumb|center|Equation 3]]
[[File:Dispersion equationvl915.PNG|400px|thumb|center|Equation 4]]

In equation 4, J is the harmonic force constant and m the mass of the atoms (i.e the spheres) in figure 1. The graph of ωk against k is called a dispersion curve. If for a dispersion curve such as the ones shown further below for example (fig. 4a) a sum over all vibrational bands and all possible k points is taken, then all vibrational levels of the crystal will have been taken into account so that a partition function can be expressed from which various thermodynamic properties of the crystal can be obtained such as the Helmholtz free energy for instance for which the expression (in the harmonic approximation) is shown in equation 5: <ref name=":0" />
[[File:Helmholtz vl915.PNG|400px|thumb|center|Equation 5]]
In equation 5 the first term represents the internal energy contribution while the second and third terms stand for the zero-point energy of the vibrations and for their entropic contribution, respectively. J and k are indices for the vibrational bands and k points, respectively.

It should also be noted that the dispersion curve described by equation 4 is periodic with a repeat period of 2π/a. This means that all useful information is contained in waves with wavevectors in the range -π/a < k < π/a. This range can be interpreted as a unit cell length in reciprocal space (k - space) and is called the first Brioullin zone and -π/a and π/a are called the Brioullin zone boundaries. As consequence of this periodicity is that vibrational bands are generally only represented for k values of the first Brioullin zone. Hence a unit cell length a* in reciprocal space can be expressed by equation 6: <ref name=":0" /> [[File:Brioullinzonevl915.PNG|400px|thumb|center|Equation 6]]

The concepts discussed so far can be extended to 3D and to unit cells containing more than one atom. In 3D, Brioullin zone boundaries remain defined as lying half-way between two reciprocal lattice points but now the boundaries are represented by planes rather than points and a Brioullin zone becomes the minimum volume enclosed by a set of such boundary planes that encloses a reciprocal lattice point. In 3D the Brioullin zone can thus be interpreted as a unit cell of reciprocal space and just as it was the case for 1D, the dispersion curves are only shown for k values (which in 3D are broken down into a kx-,ky-, and a kz component) within the first Brioullin zone. The letters on the horizontal axis in such dispersion curves (like in fig. 4a) are simply high symmetry points in three dimensional k-space such as Γ for example which represents the point (0,0,0) where kx,ky, and kz are equal to zero.

The harmonic approximation works well as long as the displacements from equilibrium (un) are small realative to a, which can be assumed to be the case at low temepratures. For small displacements from the equilibrium spacing between the atoms the harmonic approximation will describe the interatomic potential reasonably well. Meanwhile, it should be noted that at higher temperatures the displacements from equilibrium will grow and more and more terms must be included in the Taylor expansion (equation 2) to obtain an accurate description of the interatomic potential which in reality is of course far from perfectly harmonic. It can thus be assumed that the harmonic approximation will fail at higher temperatures. One consequence of not taking into acoount anharmonic terms of the interatomic potential is that the harmonic approximation fails to predict thermal expansion (more on this later). Hence if thermal expansions are to be modeled the so-called quasi-harmonic approximation is often used, which takes into account the increased equilibirum spacing between the atoms as temperature increases but apart from that treats the dynamics of the crystal within the harmonic approximation. In the quasi-harmonic approximation the equilibrium spacing between atoms at a given temperature T is determined by minimising the Helmholtz free energy (equation 5) of the crystal for that temperature.

===The molecular dynamics simulation method===
The molecular dynamics simulation (MDS) method is an alternative to the harmomic / quasi-harmonic approximation to model crystal vibrations. The essence of this method is to solve Newton's equations of motion for a set of atoms. The computer randomly attributes velocities to the atoms in the crystal until the total internal kinetic energy of the set of atoms corresponds to the temperature being investigated. It is assumed that the atoms interact via a model interatomic anharmonic potential such as the Lennard-Jones or the Buckingham potential for example. The equations of motions are solved numerically using small discrete time steps which allows to determine the positon and the velocity of the atoms at a given time. This makes it possible to generate trajectories of each atom in time from which the computer determines time-averaged properties.

The molecular dynamics methods allows to simulate crystal vibrations at high temperatures (unlike for the harmonic / quasi-harmonic approximation) and can even be used to model phase transitions as the anharmonicity of the interatomic potential is fully taken into account. <ref name=":0" /> A drawback however is that the accuracy of the simulation is limited by the accuracy of this anharmonic potential. Furthermore, this method does not take into account the zero-point energy and can thus be expected to be inaccurate for simulations at low temeperatures near 0 K. Also, the MDS method is computationally more expensive than the quasi-harmonic approximation. Finite-sized samples and non-infinitessimal time steps impose severe limitations on the accuracy of the simulations but they are a necessary comrpomise between accuracy and available computational resources. Another commonly adapted simplification to reduce computational resources not only in MD simulations is the rigid-ion model which neglects the polarizability of ions in a crystal and treats them as non-deformable charged spheres. <ref name=":0" />

===The thermal expansion coefficent===
The origin of thermal expansion of a solid material can be explained by considering a typical anharmonic interatomic potential (fig.2):
[[File:Thermal expansionvl915graph.PNG|400px|thumb|center|Fig.2: Shows an anharmonic interatomic potential which and shows how the average interatomic separation increases with increased internal energy (i.e. increased temperature)]]
The temperature of a crystal is related to the kinetic energy of its atoms. Thus increasing the temeprature is equivalent to an increase in internal energy. It can be seen from figure 2 that for a larger internal energy (assuming a realistic anharmonic interatomic potential) the amplitude of the vibration and the average interatomic distance increase which results in thermal expansion of the crystal. If the potential in figure 2 was a harmonic potential (which is parabolic and therefore symmetric) the average interatomic distance would not change. Thus the harmonic approximation, as mentioned before, does not account for thermal expansion. Moreover, thermal expansion can also be explained from the viewpoint of free energy. It can be seen from figure 2 that a consequence of thermal expansion is a gradually decreasing energy difference between the allowed energy levels. For a harmonic potential these differences are all euqal and don't change. It can thus be concluded that thermal expansion allows for more energy states to be accessible at a certain temperature which makes it entropically favourable. Hence, a material expands by increasing the temperature as the energetic advantage of a small volume is overcome by the entropically favourable expanded volume. <ref name=":0" />

Thermal expansion can be quantified with thermal expansion coefficnets. An expression for the volumetric thermal expansion coefficent of a material is given in equation 7:
[[File:Thermal expansionvl9152.PNG|400px|thumb|center|Equation 7]]
It should be noted that the volume temperature-gradient (at constant pressure) of a material depends on the volume of the material and is hence an extensive property. The normalisation with the initial volume V makes of the expansion coefficent a volume-independent intensive property which is easier to tabulate. Meanwhile, it should be mentioned that thermal expansion coefficents are not temperature independent, as the volume of a material does not increase in a linear fashion as temperature is increased, and hence tabulated coefficents from the literature are always bound to a certain temperature. <ref>P. Atkins, J. de Paula, ''Physical Chemistry,'' Oxford University Press, Oxford, 2014</ref>
=== The unit cells of MgO===
Two different unit cells are commonly defined for an MgO crystal, a conventional unit cell (fig.3a) which is cubic with a face-centred lattice and contains 4 MgO units, and a rombohedral primitive unit cell which contains 1 MgO unit (fig.3b):

<center>
{| class="wikitable"
|[[File:Mgoconvvl915.gif|300px|thumb|center|Fig. 3a: Conventional unit cell of MgO]]
|[[File:mgo_primvl915.png|400px|thumb|center|Fig. 3b: Primitive unit cell of MgO]]
|} </center>

The conventional unit cell better reflects the symmetry of an MgO crystal and is thus more commonly used. Meanwhile software often outputs calculated properties such as lattice constants or untit cell volumes for the primitive unit cell. It should be thus noted that the density of MgO remains the same of course, completly independent of the unit cell used which is important for interconversions between different unit cells as the relation shown in equation 8 holds for any type of unit cell:
[[File:Density realtionvl915.PNG|550px|thumb|center|Equation 8]]
In equation 8, mi is the total mass contained in a certain unit cell and Vi is the volume of that unit cell.

Finally, it should be mentioned that for MD simulations supercells are used which contain many unit cells. This is necessary since, as mentioned before, the MDS method uses time-averaging to determine properties of the system such as the average interatomic distance for instance. Of course, the more atoms are studied, the more meaningful and accurate the averaging becomes which is why supercells, the larger the better, are used for these simulations.

==Experimental==
===General===
The computational modelling in this experiment was done using the General Utility Lattice Program (GULP) within in the DLVisualise graphical user interface. It was made sure that for all simulations 'ionic.lib' was selected and that the 'Include Shells' option was not selected on the 'GULP potential parameters' panel, where the latter ensures that a rigid-ion model is used. In these investigations the Buckingham potential was used to model the interatomic potential.
===Computing the phonon dispersion curve of MgO===
To obtain a phonon dispersion graph for MgO the phonons are calculated at 50 points along the W-L-G-W-X-K path in three dimensional k-space. This was done by loading a conventional MgO unit cell in the DlViusalise interface. Subsequently, 'Phonon Dispersion' was selected under 'General Opts' on the 'execute GULP' panel and it was made sure that 'Npoints'''<nowiki/>', '''which defines the number of k-points at which the phonons will be computed, was set to 50 and that 'Calculate Phonon Eigenvectors' was selected.
===Computing the density of dtates with different grid sizes===
A DOS for different grid sizes was obtained by loading a conventional MgO unit cell in the DLVisualise interface and selecting 'Phonon DOS' under 'General Opts' on the 'execute GULP' panel and by specifying the shrinking factors along the a,b and c directions (i.e. the grid size). A DOS was calculated for grid sizes of 1x1x1,2x2x2,3x3x3,4x4x4,5x5x5,8x8x8,16x16x16, 32x32x32, and 64x64x64 (fig.5).
===Calculating the free energy within the quasi-harmonic approximation at different temperatures===
The free energy of MgO can be obtained from a simulation output log file which is generated during the DOS calculations. In this experiment the free energy was calculated for the same grid sizes as as for the DOS in the previous section at a fixed temperature of 300 K. To do this the same procedure as for the DOS was followed but it was additionally made sure that the temperature was set to 300 K on the 'General Opts' panel (the pressure is set 0 GPa by default and was not changed).

Moreover, the free energy of MgO was calculated at different temperatures with a grid size of 32x32x32. The temperature was increased in steps of 100 K, going from 0 K to 2900 K. The volumes of the primitive unit cell at different temperatures can be obtained from the same simulation output file as the free energies and was recorded for every temperature as it is needed for the calculation of thermal expansion coefficients (see Results and Discussion section).
===Calculating the free energy of MgO with the MDS method at different temperatures===
It was mentioned in the introduction that MD simulations on small samples such as a single MgO unit cell are meaningless which is why for these calculations a supercell of 32 MgO units was loaded into the the DLVisualise interface. Subsequently, 'Molecular Dynamics'<nowiki/>''''' '''''was selected on the 'Execute GULP'<nowiki/>'' ''panel'''''.''''' Under the''''' '<nowiki/>'''''MD Opts'<nowiki/>''''' '''''on the 'execute GULP' panel the 'ensemble'<nowiki/>''' '''was set to''''' ''''''NPT'. This fixes the number of particles N, the pressure P and the temperature T of the system, the volume is allowed to vary of course. The 'time step' was set to 1 femtosecond, the 'equilibration steps' and 'production steps' were set to 500, respectively, and the 'sampling steps' and the 'trajectory write steps' were both set to 500 as well. This calculation was repeated at different temperatures increasing the temperature in steps of 100 K, going from 0 K to 3000 K. For each temperature the total energy and the volume of the MgO supercell were recorded, which were again obtained from the simulation output log file.
==Results and Discussion==

=== Dependence of the density of states on the grid size ===
<center>
{| class="wikitable"
|[[File:Bandvl915.PNG|400px|thumb|center|Fig.4a: Dispersion graph for MgO.]]
|[[File:1x1x1picturevl915.PNG|400px|thumb|center|Fig.4b:DOS with 1x1x1 grid size.]]
|} </center>
To understand how the dispersion curves of MgO and the density of states (DOS) graphs are correlated it is important to take a closer look at figure 4: In the dispersion graph (fig.4a) six modes of oscillation per k-point (so six branches) would be expected as the primitive unit cell of MgO contains two atoms. This is the case left of the point L in figure 4a, meanwhile not on the right of the L point as two of the branches "merge" which means that at these k-points, two modes have the same energy. Looking at the density of states graph with shrinking factor 1x1x1 (fig.4b), four peaks can be seen at 286, 351, 676 and 793 cm-1, respectively. The two lower-frequency ones have twice the intensity of the latter higher-frequency two. When a branch in the dispersion graph intersects with the sampled k-point, a peak can be observed in the DOS graph at the corresponding frequency. Thus to locate the one sampled k-point that is sampled with a 1x1x1 grid, a value of k needs to be found on the vibrational bands graph which is crossed by four branches at the frequencies observed in the DOS graph. These criteria are met by the k-point labelled L in the vibrational bands graph. The stronger intensity of the two lower-frequency peaks in the DOS graph can be explained by the double-degeneracy of the two lower-frequency branches mentioned before. The modes that go to zero at Γ are called optical modes, the modes that do not come down to zero at Γ are called optical modes. In figure 4a, 3 acoustic modes and 3 optical modes can be observed .

As the grid size is increased more k-points are sampled resulting in much more peaks over a wider frequency-range in the DOS graphs (fig.5). At a grid size of 16x16x16 already, there are so many peaks that they strongly overlap so that a single broad peak is obtained. A relatively smooth curve is obtained at a grid size of 32x32x32 and it was found out that increasing the grid size beyond 32x32x32 does not result in significant changes in the DOS graphs. For instance the differences between the 32x32x32 and the 16x16x16 graph are much more pronounced than the differences between the 64x64x64 graph. Hence, from there on the 32x32x32 grid size was used for simulations, as it is only marginally less accurate but much less time-consuming than larger grid sizes.

<center>
{| class="wikitable"
![[File:2x2x2vl915.PNG|400px|thumb|left|Fig. 5a: DOS with a 2x2x2 grid size.]]
![[File:3x3x3vl915.PNG|400px|thumb|left|Fig. 5b: DOS with a 3x3x3 grid size.]]
|-
![[File:4x4x4vl915.PNG|400px|thumb|left|Fig. 5c: DOS with a 4x4x4 grid size.]]
![[File:5x5x5vl915.PNG|400px|thumb|left|Fig. 5d: DOS with a 5x5x5 grid size.]]
|-
![[File:8x8x8vl915.PNG|400px|thumb|left|Fig. 5e: DOS with a 8x8x8 grid size.]]
![[File:16x16x16vl915.PNG|400px|thumb|left|Fig. 5f: DOS with a 16x16x16 grid size.]]
|-
![[File:32x32x32vl915.PNG|400px|thumb|left|Fig. 5g: DOS with a 32x32x32 grid size.]]
![[File:64x64x64vl915.PNG|400px|thumb|left|Fig. 5h: DOS with a 64x64x64 grid size.]]
|} </center>
It can also be seen that the dispersion graph is related to the densitiy of states graph by the gradient of the branches. A steep gradient on a branch in the dispersion graph means a small DOS at the corresponding frequency, whereas a flat gradient on a branch results in a large DOS at the corresponding frequency. In figure 5i below the DOS for a 32x32x32 grid is shown rotated by 90° next to the dispersion of MgO so that this becomes easier to see.
[[File:Comparison DOS dispersionvl915.PNG|600px|thumb|center|Fig.5i: Comparison of the DOS graph and the dispersion curves.]]
It can be seen from figure 5i that the density of states is very large around 400 cm-1 which makes sense as many branches go through this frequency including one branch with a very shallow gradient in the middle which remains around 400 cm-1 throughout the whole k-path. In contrast, the density of states is considerably lower at the extremes looking at 200 cm-1 and 1000 cm-1 for instance. Much less branches go through these frequencies and if they do then with a quite steep gradient.
=== Dependence of the free energy on the grid size ===

[[File:Differencesvl915.png|600px|thumb|center|Fig.6: This graph shows the absolute difference between the free energy and the converged free energy at different grid sizes.]]

<center>
{| class="wikitable"
!Grid size
!Free energy / ev
!Differnece from converged
free energy / 105 eV
|-
|1x1x1
|<nowiki>-40.930301</nowiki>
|381.8
|-
|2x2x2
|<nowiki>-40.926609</nowiki>
|12.6
|-
|3x3x3
|<nowiki>-40.926432</nowiki>
|5.1
|-
|4x4x4
|<nowiki>-40.92645</nowiki>
|3.3
|-
|5x5x5
|<nowiki>-40.926463</nowiki>
|2.0
|-
|8x8x8
|<nowiki>-40.926478</nowiki>
|0.5
|-
|16x16x16
|<nowiki>-40.926482</nowiki>
|0.1
|-
|32x32x32
|<nowiki>-40.926483</nowiki>
|converged
|-
|64x64x64
|<nowiki>-40.926483</nowiki>
|converged
|}
Table 1: This table shows the grid size, the corresponding free energy and the absolute difference from the converged free energy.</center>

Equations for thermodynamic functions can be evaluated using a list of frequencies taken from the dispersion curves calculated over a fine grid of k-points within the first Brioullin zone. Since energies of vibrations depend on their frequency only, the density of states can be used for the evaluation of energies. <ref name=":0" /> It follows that the more accurate the density of states used the more accurate will be the energy value obtained. This implies that a finer grid of-k-points sampled, will improve the accuracy of the free energy obtained. However, just as with the DOS functions in the previous section, there is a certain threshold-grid size beyond which the energy values obtained do not change significantly. In fact in figure 6 and table 1 it can be seen that the energy values obtained converge, which is unsurprisingly at the grid-sizes for which the DOS graphs in the previous section were found out to converge as well. It can be seen from table 1 that at a grid size of 8x8x8 calculations are accurate to 0.5 eV already, at 16x16x16 to 0.1 eV. For a grid size of 32x32x32 the same energy value is obtained as for the much more time consuming 64x64x64 calculation. This underpins again that calculations with grid-sizes above a certain threshold result in unnecessarily long computation time without providing significalntly more accurate results.

Meanwhile, what that threshold grid size is depends on the sytsem studied. From equation 6 it can be seen that the larger the unit cell length a, the shorter will be the unit cell length in reciprocal space a* which is directly related to the size of the Brioullin zone. It can be assumed that less k-points need to be sampled to "explore" a smaller Brioullin zone and that the aforementioned threshold grid size will hence be lower, the smaller the Brioullin zone. Conversely, the larger the Brioullin zone of a crystal system, the larger will be the threshold grid size.

A CaO crystal for example, which just as MgO has a face-centred cubic unit cell, can be expected to have a larger unit cell length than MgO since Ca has an additional electron shell and will hence have a larger ionic radius than magnesium. Consequently, CaO will have a smaller Brioullin zone than MgO and will thus have a smaller threshold grid size.

Similarly for Zeolites such as Faujasite for example. Zeolites are porous aluminosilicate minerals and have very large unit cells. Faujasite for example has a unit cell length of 24.6 Å.<ref name=":1">G.Henderson, ''The Cambridge handbook of earth science data, C''ambridge University Press, Cambridge, 19968
</ref> For comparison the conventional MgO unit cell has a lattice constant of 4.2 Å.<ref name=":1" /> Hence, by the same token as for CaO, Faujasite can be assumed to have a smaller threshold grid size than both MgO and CaO.

Lithium is a metal and the bonding in metals ( metallic bonding) is commonly described as cations being held together by a shared "sea" of electrons. Since metallic bonding is non-directional the cations in metals (which can be treated as hard spheres) adopt a close packed structure for which the unit cell parameter can be expected to be smaller than for a metal oxide such as MgO. Thus Lithium can be expected to have the largest Brioullin zone and hence the largest threshold grid size of all solids discussed up until now.

=== The Helmholtz free energy and thermal expansion within the quasi-harmonic approximation ===
<center>
{| class="wikitable"
|[[File:Helmholtztempvl915.png|400px|thumb|center|Fig.7a: Shows how the Helmholtz free energy varies when increasing the temperature.]]
|[[File:clptempcl915.png|400px|thumb|center|Fig.7b:Shows how the conventional lattice parameter varies when increasing the temperature.]]
|} </center>

The two graphs in figure 7 show how the Helmholtz free energy of the MgO crytsal and the conveltional lattice parameter vary as temperature is increased in the framework of the quasi-harmonic approximation. It can immediately be seen from figure 7a that the Helmholtz free energy of the MgO crystal decreases with increasing temperature. And from ~1000 K onwards it even follows a linear trend. Considering the most common expression for the Helmholtz free energy, A = U-TS, this suggests that at least from a 1000 K onwards the entropy S and the internal energy U are constant and therefore temperature independent. This is surprising as if for instance the third term in equation 5 is taken into focus, which should be equivalent to the "-TS" term, two T depedencies are found suggesting that S is temperature dependent. Furthermore, the internal energy, E0 in equation 5, even though it is temperature independent, can be expected to change for every temperature as within the quasi-harmonic approximation the equilibrium interatomic spacing is adjusted at every temperature. Hence the linear trend of the Helmholtz free energy could be explained by assuming that the temperature dependencies of U and S become negligible at high temperatures. Conversely, at temperatures below a 1000 K the dependencies of U and S can be assumed to have a more significant effect so that a linear relation between the Helmholtz free energy and temperature is not established.

The overall decreasing trend of the free energy as temperature is increased can be rationalised considering that the temperature T is always positive (as the Kelvin scale is used) and that S, the absolute entropy, is always positive as well. As mentioned, the internal energy is not found to vary strongly with temperature, especially not at temperatures above a 1000 K, and hence increasing T leads to a decrease in the Helmholtz free energy.

The conventional lattice parameter was obtained by taking the cube-root of four times the primitive unit cell volume (which was obtained from the output log file) using equation 8 and the fact that a convetnional unit cell has the fourfold mass of a primitive one. As would be expected the conventional lattice parameter increases with increasing temperature as a consequence of thermal expansion. At low temepratures (0 K - 300 K) the lattice parameter is harldy changing at all which is not a surprising if one looks at how the Helmholtz energies evolve at these temperatures, they are hardly changing. Since the lattice paramter is obtained from the unit cell volume which in turn is obtained from the minimum Helmholtz free energy at that temperature, the lattice parameter will of course hardly change if the Helmholtz free energy hardly changes. In the intermediate regime at temepratures from 500 K to 2200 K the lattice paramter follows a roughly linear trend. Meanwhile, above 2200 K the lattice parameters starts to increase more steepily with temperature.

=== Comparison of the thermal expansion modeled with the quasi-harmonic approximation and the MDS method ===
[[File:Expansionvl915harmds2.png|400px|thumb|center|Fig.8: This figure shows how the convetnional unti cell volume is predicted to evolve as temperature increases comapring the quasi-harmonic approximation and with the MDS method.]]

Figure 8 above shows how the convetnional unit cell volume is predicted to eveolve as temperature is increased and compares the predictions of the of the quasi-harmonic model and the MDS metod. As mentioned in the introdcution the quasi-harmonic approximation, unlike the MDS method, takes into account the zero-point vibrational energy. Hence the MDS method tends to underestimate unit cell volumes at low temepratures. Indeed, it can be seen in figure 8 that at low temperatures (0 K - 400 K) the volumes predicted by the MDS method are lower than the ones predicted by the quasi-harmonic approximation. In the intermediate regime at temperatures between 500 K and 2000 K the two models appear to agree quite well while at 2000 K the two curves diverge from each other quite strongly. The unit cell volumes predicted by the quasi-harmonic approximation start to increase more steepily from 2000 K onwards than the MDS volumes. This divergence of the two models at high temperatures can again explained by the fact that anharmonic terms of the interatomic potential become important at high temperatures; the MDS method accounts for these terms while the quasi-harmonic approxiamtion doesn't which is why the two curves in figure 8 diverge at high temperatures.

The MDS method can in principle successfully model crystal vibrations at high temepratures and can even be used to model phase transitions in some cases. Meanwhile, from figure 8 it can be seen that the MDS method starts to break down around 1300 K already as the predicted volumes start to hover up and down. These unwanted fluctuations could be resolved by using a larger supercell or a smaller time step for the MD simulations for the following reason: Increasing the temperature results in more vigorous vibrations of the ions and hence there will be stronger fluctuations in the respective positions of the ions between each time step. The use of a larger supercell or of a smaller time step would thus improve the reliabilty of the time-averaged interatomic spacing determined by the software and prevent unwanted fluctuations of the predicted unit cell volumes as seen in figure 8.

As typically up to 10 steps along a given vibration is accepted to be sufficient using a time step of one femtosecond in combination with a supercell containing 32 MgO units seemed to be adequate to use in this investigation. Meanwhile, as the MDS method's results break down at a relatively low temeprature, far away from the melting temperature of MgO (3125 K), it should be tried in future investigations to carry out these simulations with different settings. For instance, a similar procedure to the optimal grid size investigation for the quasi harmonic approximation could be carried out. The time step could be gradually shortened (while keeping the super cell size constant) until the total energy of the simulations start to converge. The same procedure could also be carried out by varying the super- cell size and keeping the time step constant. This would possibly allow to determine settings for which the MDS method does break down at higher temperatures than around 1300 K which would yield more consistent results in the temperature range studied in this investigation. Meanwhile, as often the case for computational simulations, a compromise must be kept in mind between accuracy and available computational resources.

=== Determining the thermal expansion coefficents from the modeled thermal volumetric expansions ===
As can be seen from equation 7 the volumetric thermal expansion coefficient is related to the gradient of a volume - temperature plot. Hence an expression for the expansion coefficients at different temperatures can be obtained by taking the first derivative of a polynomial regerssion of the data points in figure 8 (without forgetting to normalise this expression of course). A set of experimentally determined thermal volumetric expansion coefficients for an MgO crystal was found in the literature<ref>O.L. Anderson, K. Zou, ''J. Phys. Chem. Ref. Data'', 1990, '''19'''(1), 69-81</ref> for a temeprature range between 300 K and 2000 K. Hence a polynomial regression was attempted over this temeprature range for the quasi-harmonic and the MDS method data points, respectively. A linear regression would only be possible over short temperature ranges and result in a temperature independent expansion coefficient. Hence a polyanomial regression was preferred as it could be carried out over the whole temperature range of interest (300 K - 2000 K) and because it would allow to see how the computed expansion coefficients evolve as temeprature is altered since the first derivative of a non-linear polynomial is of course not a constant. Furthermore, the higher the order of the polynomial regression, the larger the coefficient of determination (the R2 value) which means that the regression better describes how the individual data points are interconnected. Polynomial regressions of third order were assumed to be sufficient as they resulted in R2 values over 99 % (see table 3) for both, the MDS and the quasi-harmonic set of data points. The polynomial regressions are shown in figure 9 and the corresponding R2 values in table 3:
<center>
{| class="wikitable"
|[[File:Linear fit expansion coeffeicnet quasi harmonic2.png|400px|thumb|center|Fig. 9a: Polynomial regression (3rd order) of the quasi harmonic data points.]]
|[[File:Linear fit expansion coeffeicnet md2.png|400px|thumb|center|Fig. 9b: Polynomial regression (3rd order) of the MDS data points.]]
|} </center>

<center>
{| class="wikitable"
|
|R2
|-
|Quasi-harmonic (fig. 9a)
|0.999982
|-
|MDS (fig. 9b)
|0.991150
|} </center>
<center> Table 3: This table shows the coefficients of determination for the third order polynomial regressions on the MDS and quasi-harmonic data points, respectively. </center>

The first derivative of these third order polynomials (table 4 ) was determined and normalised by the volume predicted at 300 K for the MDS method and the quasi-harmonic approximation, respectively. This resulted in the two second order polynomials shown in table 4 which describe how the computed thermal volumetric expansion coefficients evolve as temperature is altered (in the temperature range from 300 K to 2000 K) according to the quasi harmonic approxicmation and the MDS method, respectively.

<center>
{| class="wikitable"
|
|Polynomial regressions
|αV(T)
|-
|Quasi-harmonic
|4.020 * 10-11 T3 + 3.638 * 10-7 T2 + 1.595* 10-3 T + 75.035
|1.596 * 10-12 T2 + 9.630 * 10-9 T + 2.111*10-5
|-
|MDS
|3.999 * 10-10 T3 - 1.306 * 10-6 T2 + 3.562 * 10-3 T + 74.282
|1.592 * 10-11 T2 - 3.468 * 10-8 T + 4.728 * 10-5
|} </center>

<center> Table 4: This table shows the polynomial regressions and αV(T) determined from the MDS and quasi-harmonic data points, respectively. </center>

The expressions obtained in table 4 were then used to plot the detemined values of αV(T) together with the experimental expansion coefficients from literature which can be seen in figure 10a. Figure 10b shows the error of the computed expansion coefficients (in %) from the experimental expansion coefficients in the temeprature range 300 K - 2000 K for the quasi harmonic approxiation and the MDS method, respectively.

<center>
{| class="wikitable"
|[[File:blavl915.png|500px|thumb|center|Fig. 10a: This graph shows how the thermal volumetric expansion coefficient evolves with temperature as measured experimentally and as predicted by the quasi-harmonic approx. and the MDS method, respectively.]]
|[[File:Main graph 2vl9152.png|500px|thumb|center|Fig.10b: This graph shows the error in % of the thermal expansion coefficients as predicted by the MDS method and quasi harmonic approx. from the experimental values in a temeprature range from 300 K to 2000 K]]
|} </center>

From figure 10 it can be seen that both models, the quasi-harmonic and the MRD method underestimate the thermal expansion coefficient in the temperature range of 300 K to 2000 K. Moreover, figure 10 shows that throughout the whole of this temperature range the quasi harmonic approximation is closer to the experimental values and it can therefore be concluded that the quasi-harmonic approximation better estimates the expansion coefficients than the MDS method in the temperature range of 300 K to 2000 K. From figure 10b it can be seen that the quasi-harmonic approximation is the furthest off from the experimental values in the range from 500 K to 1000 K but starts to better estimate the expansion coefficients below 500 K and beyond 1000 K. The MRD method expansion coefficients are the furthest off from experimental values in the temperature range from 1250 K to 1500 K and start to better estimate the expansion coefficients beyond and below this temperature range. The modeling with the MDS method could probably be improved significantly by increasing the supercell size or by shortening the time step and it should be a subject of further investigation to show this.

== Conclusion ==
In this investigation the lattice vibrations of an MgO crystal were modeled at different temperatures to predict the thermal expansion of the crystal using computer simulations. This task was approached in two different ways namely with the quasi-harmonic approximation and with the MDS method . It was furthermore found out that a k point grid size of 32x32x32 is a good compromise between accuracy and computation time to model the DOS of the MgO crystal from the crystal's dispersion. Increasing the grid size beyond 32x32x32 was shown not to result in signifacant changes in the appearence of the DOS and it was shown that the free energy values calculated at different grid sizes also unsurprisingly converge from a gridsize of 32x32x32 onwards. Moreover, it was shown that the thermal expansion predicted by both models agrees well in a temperature range from 500 K to 2000 K but below and beyond that range the quasi-harmoni model predicts larger volumes than the MDS model. Temperature dependent expressions for the thermal volumetric expansion coefficients of MgO were determined from third-order polynomial regressions of the volume-temperature curves of both models in a temperature range from 300 K to 2000 K. By comparison with experimentally determined expansion coefficients from literature it was shown that throughout this temperature range the quasi-harmonic approximation significantly better models thermal expansion than the MDS method.

== References ==

File:Blavl915.png

2018-03-23T15:57:24Z

Vl915:

File:Main graph 2vl9152.png

2018-03-23T15:56:05Z

Vl915: Vl915 uploaded a new version of File:Main graph 2vl9152.png

2018-03-20T19:31:44Z

Vl915:

2018-03-20T19:04:18Z

Vl915:

Bla

2018-03-20T19:03:52Z

Vl915:

<center>
{| class="wikitable"
|
|
|Quasi-harmonic
|5.321 * 10-13 x3 + 4.815 * 10-9 x2 + 2.111 * 10-5 x + 0.993
|-
|MDS
|5.308 * 10-12 x3 - 1.734 * 10-8 x2 + 4.728 * 10-5 x + 0.986
|} </center>