2018-02-27T15:21:59Z

Ar3015:

<h1>MgO Report - Antonia Reeves</h1>
<h2>Introduction</h2>
In chemistry, it is much easier to treat crystals as rigid structures with set lattice points in unit cells. This can often hold as a good approximation at very low temperatures. However, this is obviously not the case. There are many different vibrations happening throughout the lattice structure leading to thermal expansion. These vibrations can be quantified by using a wave vector - this gives the direction of the propagation of the wave in reciprocal space. We use reciprocal space because FILL IN HERE.
<h3>Phonons and phonon dispersions</h3>
In quantum chemistry, we can treat vibrations as both particles and waves. When treating vibrations as particles - just as particles of light are photons - they are called phonons. We can plot our frequency CHECK THIS of the phonons against the wave vector. This effectively gives a plot of the energy against the 'the way that vibrations move' INSERT DIAGRAM OF POTENTIAL WAYS TO DEMONSTRATE which is the k vector. As the phonon dispersion plot investigates the different ways vibrations can take place in the same time, if there were for example only a few atoms, this would give rise to only a few different k points. However, as we are dealing with an infinite lattice, instead of giving quantised k points, a continuous line of k points can be generated, giving rise to a band. One band is called a branch. The number of branches is dependent on the number of different atoms in the lattice. In MgO lattice there are 2 atoms. Using the equation 3N (not 3N-6 as this takes into account rotational and translation degeneracy for single molecules that we don't need to consider for crystal vibrations) where N is the number of different atoms, we get 6 branches which can be seen in the phonon dispersion plot in Figure X. Using the dispersion plots are useful for seeing the frequency vs k vector as a whole, however, makes it difficult to determine pieces of information such as how many states are at a given frequency, and whether there is a direct or indirect band gap.
<h3>Thermal Expansion Coefficient</h3>
Using the equation <math>\frac{1}{V_0}\frac{\delta V}{\delta T}</math> gives the thermal expansion coefficient. This is a value that tells us how much a given crystal expands with an increase in temperature. We use the <math>\frac{1}{V_0}</math> term to normalise this value so it can be applied to a crystal of any size which is composed of the same crystal. Below, the different computational methods that can be used to investigate how these properties change with temperature are used.
<h3>Harmonic and Quasi-Harmonic approximations</h3>
The harmonic oscillator is simply a taylor expansion of INSERT EQUATION. The harmonic oscillator takes the third term (quadratic term) giving the easily the harmonic oscillator parabola. When looking at harmonic potential energies, as we change the temperature, the average bond length stays at an equilibrium point and therefore the volume of the crystal doesn't change. This means that although this approximation works well at lower temperatures as there is little change in volume at low temperatures in reality, at higher temperatures, harmonic and anharmonic (which is a much more realistic model) deviate a lot, giving unreliable data for higher temperatures. A compromise between harmonic which isn't a very reliable approximation as it only uses 3 terms from the taylor series expansion and using an anharmonic which uses many terms and will be very expensive is using a quasi-harmonic approximation. When we use the Quasi Harmonic approximation, although the approximation isn't as accurate as anharmonic, it still gives the point where the bond breaks at higher temperatures. INSERT HELMHOLTZ EQUATION. In the helmholtz free energy equations we use, the differences in the equations for harmonic and quasi-harmonic plots are in the frequency term. In a harmonic system, we would use equation 3 as the potential, however as we are looking at a quasi harmonic system, within the brillouin zone <math>-\pi < k < \pi </math> we use equation 1.
<br><math> w_k = (\frac{4J}{m})^{\frac{1}{2}} |sin \frac{ka}{2}| </math> <b>(Eq. 1)</b>
<br><math> w_k = (\frac{2J}{m})^{\frac{1}{2}} |sin \frac{ka}{2}| </math> <b>(Eq. 2)</b>
<br><math> w_k = (\frac{1}{2\pi}){\frac{k}{\mu}}^{\frac{1}{2}} </math> <b>(Eq. 3)</b>
<h3>Molecular Dynamics</h3>
The molecular dynamics method takes the values from the ideal MgO crystal, then when we change the parameters, generate random velocities. This process is continued at each time step until the average of the values converge to a single point. Once the values of the velocity has started to converge, can start to measure properties of the crystals e.g. Energy.
<h3>Helmholtz Free Energy</h3>
[[File:Screen_Shot_2018-02-27_at_10.12.39.png|thumb|centre|Eq. 4]]
A = U - TS
In order to look at the energy of the vibrational modes, we use the vibrational contribution to the entropy using statistical thermodynamics and the partition coefficient.
The first term of equation 4 is the internal energy, the second term is the zero point energy (energy at 0K) and the final term is the statistical thermodynamics contribution - the contribution to the energy from the different vibrations.
<h2>Method</h2>
For all of the calculations, the GUI DLVisualise was used. For all calculations, the execute GULP panel was brought up and in the GULP potential parameters, ionic.lib was used and include shells was not included.
<h3>Lattice vibrations -Computing the phonons</h3>
Single point was selected, and then run was pressed. Once the job has completed, this shows in the job list and the files were recovered. The energy could be found in the log file under Total Free Energy.
[[File:MgO_phonon_1.PNG|thumb|centre|Figure 3: Aide for bond lengths table]]
<h3>Calculating the free energy in the Harmonic approximation</h3>
Single point was selected. Under General Opts, phonon dispersion was selected with N points set to 50 and calculate phonon eigenvectors was selected, so that specific points in k-space (relating to the origin, x direction vectors being in phase etc.) are represented with greek letters. The files were recovered and the phonon dispersion was displayed using the 1D data display. This is simply one way of displaying all of the data - the same data can be displayed as a density of states. Instead of simply giving the angular frequency as a function of the wave vector, a plot of density of states can also show how many states have a given energy at each k vector point. This can show the energy difference in the band gap, and the most densely populated k-states more easily. In order to visualise the density of states, the phonon dispersion is 'chopped up' into sections. The number of sections is determined by the grid size (these parameters can be changed in the GULP execution panel. In this, the frequencies at the edge of the grids are taken and plotted (only one side as not to repeat the counting of the frequencies). The higher the grid size, the greater the accuracy of the density of states as this allows a greater number of wave vectors to be sampled. The value of the grid size is called the shrinking factor.
<h3>Thermal expansion of MgO</h3>
In order to determine a suitable shrinking factor, a number of grid sizes were investigated with the highest being 32x32x32. As having a higher shrinking factor takes longer and more processing power (leads to being more expensive), choosing the shrinking factor becomes a compromise between time and accuracy. As can be seen from the plot, in Figure [X], and a subsequent calculation to check this, the Helmholtz free energy stops changing after a grid size of 20x20x20 so this seemed like a reasonable value to use as not to compromise the accuracy for a quicker run time. The settings used for this calculation were, select phonon DOS, change shrinking factor to 20x20x20 and set the temperature from 0-3000 K in 100 K intervals. For each temperature, the free energy and BLAH were taken from the log file.
[[File:MgO_DOS_16x16x16.PNG|thumb|centre|Figure 3: Aide for bond lengths table]]
<h3>Molecular Dynamics</h3>
For Molecular dynamics, a super cell which is twice the length on every side of the cell. This is so the cell is larger and is more similar to 'an infinite lattice'. Molecular Dynamics option was selected on the execute panel, the NPT ensemble was chosen, the time step used was 1 femtosecond, equilibration steps of 500, production steps to 500 steps and sampling steps and trajectory write steps to 5. The temperature used went from 100-3000 K with 100 K interval steps. As the volume of the cell is different to the cell used in the quasi-harmonic approximation, the cell volume taken from the log file (the last average cell volume), the volume was divided by 8 to get the same volume as the conventional cell for the quasi-harmonic approximation.
<h2>Results and discussion</h2>
<h3>Phonon Modes of MgO</h3>
When the density of states for calculated for a single k point, it can be worked out by examining the phonon dispersion which k point this is. In this case, the point is L as there are only 4 peaks (2 that have double the intensity due to the double degeneracy) and at L there are only 4 branches where we would expect 6 (due to 3N to calculate number of vibrational branches expected). This can also be seen by comparing the frequencies of the peaks in the 1x1x1 DOS and the frequency of the 4 branches at L. This can be checked in the log file as although the character L is useful for helping us to determine the symmetry of the modes of vibrations, it's always useful to numerically quantify values. This is simply in the row 'K point' where the K point is 0.5. It's important that it can be viewed in the log file as on the phonon dispersion, only greek characters which represent different k points can be seen and it's not always obvious what numerical value goes with each character.
<table style="width:100%">
<tr>
<th>[[File:MgO_DOS_1.PNG|thumb|centre|Figure X: MgO Density of States 1x1x1 grid size]]</th>
<th>[[File:MgO_DOS_3x3x3.PNG|thumb|centre|Figure X: MgO Density of States 3x3x3 grid size]]</th>
<th>[[File:MgO_DOS_8x8x8.PNG|thumb|centre|Figure X: MgO Density of States 8x8x8 grid size]]</th>
<th>[[File:MgO_DOS_32x32x32.PNG|thumb|centre|Figure X: MgO Density of States 32x32x32 grid size]]</th>
</tr>
</table>
As can be seen from the progression in grid size that as the grid size increases, although the peaks with the greatest intensity stay in the same place, the shape of the density of states becomes much more defined and includes information on the states at k-points even where there is a lower number of states. When looking at the difference between grid sizes qualitatively, there isn't much of a difference between 8x8x8 and 32x32x32 apart from there is a more continuous curve with the same shape. It seems that a grid size of approximately 12x12x12 would be the minimum grid size as this is the grid size at which all k points are sampled and there aren't k-points with an angular frequency of 0 which then have non-zero values of angular frequencies at higher grid sizes. However, when deciding qualitatively, the convergence of the free energy can be used. At 20x20x20 and 6 decimal places, the free energy converges suggesting this is the optimum grid size to use. The density of states in terms of the phonon dispersion curve is basically plotting how many k points there are per value of angular frequency. This can be seen more easily by placing them side by side as seen in FIGURE ENTER HERE.
<h3>Free energy varied with grid size</h3>
INSERT TABLE HERE
As explained in the introduction, the grid size is simply how many pieces the phonon dispersion has been separated into for the analysis of the density of states. As can be seen in TABLE X, after the grid size is 20x20x20, the values have converged and don't change (when looking at 6 decimal places, it may be changing beyond 6 decimal places). This is optimal as it reduces the computing power needed for a calculation whilst retaining almost the same level of accuracy. Using a grid size of 3x3x3 will gives a suitable energy for calculations accurate to 1 meV, REMEMBER TO DO 0.5mEV!!!!!!!! and a grid size of 5 gives a suitable energy for calculations accurate to 0.1 meV.
<h3>Quasi-Harmonic Thermal expansion</h3>
<table style="width:100%">
<tr>
<th>[[File:Lattice_Parameter_Graph.png|thumb|centre|Figure X: MgO Density of States 32x32x32 grid size]]</th>
<th>[[File:Free_energy_vs_T.png|thumb|centre|Figure X: MgO Density of States 32x32x32 grid size]]</th>
</tr>
</table>
In Figure X, it can be seen that the free energy is constant, then almost linearly decreases. At very low temperatures, the only main contribution to the free energy is the zero point energy as this is temperature independent. At around 400 K (where the free energy starts to decrease). We can explain by looking at equation X as when the temperature increases, the entropic contribution increases. If we assume the internal energy doesn't change much with temperature, then the TS term will increase in magnitude and decrease the overall free energy term. By looking at equation X, the vibrational contribution term (also entropic term) is a ln term which will be negative, and again, decrease the overall free energy with an increase in temperature.
In the graph of lattice parameter vs temperature, the lattice parameter is initially constant and then increases linearly, then exponentially. After the exponential region, it can be assumed that the lattice parameter increases exponentially, presumably converging at the melting temperature, however this would need to be investigated to be confirmed. Using the quasi-harmonic approximation, this gives the potential a 'quasi-anharmonic' shape, meaning that the bond length will actually change with temperature but will just never reach the point where it breaks. Although this doesn't allow phase transitions to be investigated, it does allow the lattice parameter as a function of temperature to be investigated.

Using the Quasi Harmonic approximation, and fitting the data to get a line of best fit with a gradient of INSERT HERE which is equivalent to <math>\frac{\delta V}{\delta T}</math>. Multiplying by the reciprocal of the initial volume of the crystal, gives a thermal expansion coefficient for the region of 300-1300 K of <math>3.083\times 10^{-5}</math>. By taking the value from lower temperatures, this gives a higher accuracy compared to the 'real value'. This is because of the computational model that is being used. This model is quasi-harmonic which is based on the harmonic model. The harmonic oscillator is infinite and the bond between two atoms infinitely extends. However, at high temperatures this is not the reality. The bonds extend and break which is shown by the anharmonic model.
<table style="width:100%">
<tr>
<th>Temperature range (K)</th>
<th>Thermal expansion coefficient - Quasi-Harmonic</th>
<th>Thermal expansion coefficient - Molecular Dynamics</th>
<th>Thermal expansion coefficient - Literature</th>
</tr>
<tr>
<td>0-300</td>
<td><math>9.572\times 10^{-6}</math></td>
<td><math>1.942\times 10^{-5}</math></td>
<td><math>1.042\times 10^{-5}</math></td>
</tr>
<tr>
<td>300-1300</td>
<td><math>3.083\times 10^{-5}</math></td>
<td><math>3.140\times 10^{-5}</math></td>
<td><math>4.000\times 10^{-5}</math></td>
</tr>
</table>
[[File:QH-300-1300.png|thumb|centre|Figure X: ]]
[[File:MD-Thermal-expansion-100-300.png|thumb|centre|Figure X: ]]
[[File:QH-0-300.png|thumb|centre|Figure X: ]]
[[File:MD-300-1300.png|thumb|centre|Figure X: ]]
<h2>Conclusion</h2>