Lattice dynamics and Molecular dynamics methods were used to compute structure and to calculate the thermal expansivity of MgO. For the two techniques, this was done over a variety of temperatures, from 0K to 1000K. An optimal 32x32x32 grid size was established to run the simulations. For the Lattice Dynamic simulation, used in the quasi-harmonic approximation, was found to be 2.6709 x 10-5 K-1. For the Molecular Dynamic simulations, was found to be 3.3068 x 10-5 K-1. Both values differed from the literature value, due to the limitations of the methods. Lattice Dynamics and Dynamics Methods were compared, it was determined that LD works better at low temperature and that MD give more accurate results at high temperature.
A crystal can be defined as an effectively infinite repetition of groups of atoms.[1] The smallest building block that can be used to represent a crystal structure is called a primitive unit cell, which corresponds to a single lattice point. From this primitive cell several properties of the crystal can be determined, such as density and thermal expansion.
The density and thermal expansion of those solids are crucial properties. By studying the temperature dependence behaviour of such systems, a better understanding of their other properties, like dielectric constant, thermal conductivity and heat transfer, can be achieved.[2] This is particularly important when it comes to construction. Indeed, when choosing a material for an edifice, thermal expansion needs to be taken into account. This is due to the fact that fluctuations in temperature might lead to the collapse of the construction.
Thermal expansion is caused by the vibrational response of a solid to the thermal stresses when temperature is changed.[3]
Two types of simulation, Lattice Dynamics (with the quasi-harmonic approximation) and Molecular Dynamics, were used to analysis the thermal expansion behaviour of MgO with DL visualize (GULP command). Molecular Dynamics is based on classical mechanics equations and recorded the atoms' vibrations as kinetic energy is added to the system. Lattice Dynamics recorded the phonon dispersion and the density of state (DOS) as the temperature was increased. The experiment started by simulating the density of state for a range of grid sizes in order to find the optimal condition. This optimal grid size, 32x32x32, was then used to run the following simulations. A higher grid size would have been more accurate but was more time consuming. By recording the Helmholtz free energy and volume of the system over a temperature range (from 0K to 1000K) the coefficient of the thermal expansion of the system was determined. The 32 unit MgO cell (Figure 2), which is much bigger than the primitive cell, had to be run using the Molecular Dynamics simulation. Whereas results for the primitive cell (Figure 1) were obtained only by using the Lattice Dynamics simulation, with quasi-harmonic approximation. In both simulations, changing the temperature of the system triggered a change in volume, which when plotted in equation 1, gave the thermal expansion coefficient.
Eqn. 1 , where V0 is the initial volume and where the change in volume, , over the change in temperature, , is pressure constant
Both calculations presented different values, which made them interesting to compare between them and between a literature value.
Experimental
General
Each calculation was run using the software DL Visualize, a graphical user interface for modelling, which allows visualization of materials structures and properties. All calculations were performed in isobaric ensemble. All the simulations had ionic.lib selected for the potential model. The MgO, primitive unit cell structure, was loaded and calculations were done by executing the GULP (General Utility Lattice Program) command. The phonon dispersion curves were obtained by setting a N points of 50. The phonon density of state was computed for several shrinking factors until reaching the optimal grid size: 1x1x1, 2x2x2, 3x3x3, 4x4x4,...,32x32x32. This optimal grid size (32x32x32) was then used in the following simulations.
Lattice Dynamics, with the Quasi-Harmonic Approximation
The Quasi-Harmonic Approximation was used to compute the vibrational energy levels for the MgO, primitive cell, and to understand its phonon dispersion. Using a 32x32x32 grid, the Free Energy, Volume and Lattice parameter of the system where recorded.
The Quasi-Harmonic Approximation follows a Taylor expansion (equation 4), in which the higher terms are neglected (called anharmonic terms). By neglecting the anharmonic terms, it allows the energy of the lattice to be the same as the set of harmonic oscillator. [4]
Eqn.4
Molecular Dynamics
For 32 MgO a different type of calculation was used: Molecular Dynamics. The reason this type of simulation was not used to calculate the MgO primitive cell is that a larger number of atoms is required to obtain higher vibration degrees of freedom. All the simulations for Molecular Dynamics were done on a 32x32x32 grid at a variety of temperature: 0K to 1000K. The different parameters were adjusted: time Step was set at 1.0 femtoseconds, the number of equilibration steps and production steps were both 500. Finally the sampling steps and Trajectory write steps were picked to be 5. The ensemble was step to a fix number of particles, a fix external pressure and a fix temperature (NPT).
Molecular Dynamics is based on solving the Newton's equation of motion, , for a set of atoms in the anharmonic interatomic potential. [5] This is a stepwise process; the velocities and positions of atoms are determined from the previous velocities and positions obtained. From those data computed, the Molecular dynamics method can then calculate an accurate classical trajectory for this set of atoms.
Results & Discussion
Computing the phonon density and determining an ideal grid size for the DOS
Figure 3: Band structure of MgO, also referred as Phonon Dispersion CurvesFigure 4: Density of State ; from top to bottom, left to right; grid size: 2x2x2, 3x3x3, 4x4x4, 16x16x16, 24x24x24 32x32x32
The phonon dispersion of MgO was simulated for 50 k points. The graph obtained plotted frequency (cm-1) against symmetry points; W, L, Γ, X, K (where Γ is the origin). The phonon dispersion graph presents 6 vibrational bands, this is due to the fact that the study was done on MgO, a 2 atoms molecule, in 3D (). Those 6 branches can be divided into two categories: acoustics and opticals. The acoustic branches refer to the ones leading to Γ value of 0, whereas the optical ones represent the other branches. By looking at the Density of state for the 1x1x1 grid size, which is defined as the number of states per unit energy at a given energy, it has possible to determine the k-point. This k-point was found to be L (0.5,0.5,0.5). Indeed, the DOS graph presents 4 peaks at 288.49, 351.76, 676.23, 818.82 cm-1, which match with the frequency of L found on the phonon dispersion graph. The degeneracy of the states is reflected in the intensity of the peaks. The peaks at 288.49 cm-1 and 351.76 cm-1 are twice the size of the two others, at 676.23 and 818.82 cm-1 . This is due to the fact that there are degenerate, and thus the frequencies overlap. By looking at Table 1, it can be determined that for a calculation accuracy to 1meV, 2x2x2 grid size is appropriate. Whereas for an accuracy of 0.5 meV and 0.1meV a 3x3x3 grid size is sufficient.
Variation of Helmholtz Free Energy with grid size
Grid Size
A (eV)
ΔA (eV)
1x1x1
-40.930301
N/A
2x2x2
-40.926609
0.003692
3x3x3
-40.926432
0.000018
4x4x4
-40.926450
0.000028
8x8x8
-40.926478
0.000040
16x16x16
-40.926482
0.000040
24x24x24
-40.926483
0.000001
32x32x32
-40.926483
0.000000
64x64x64
-40.926483
0.000000
Figure 7: Graph of Helmholtz Free Energy against Grid sizes of the Lattice Dynamic simulation
Unfortunately, a single K point does not allow a realistic representation of the DOS. A wide range of grid sizes was computed in order to find an optimal condition (Figure 2). This optimal condition was determined by looking at the Helmholtz Free Energy of the different grid sizes. Table 1 and Figure 6 show an increase in the Free Energy as the grid sizes increases. This is due to the fact that a bigger grid is synonym of more vibrations and therefore more energy. The Helmholtz Free Energy (Table 1) presents a stable value at 32x32x32, which was determined to be the minimum for a reasonable approximation to the density of states. This stabilization of the grid size can by visualized on the DOS graphs. The smaller grids, like 1x1x1 , 2x2x2, 3x3x3 and 4x4x4, only presented a few peaks, which means less frequencies. As the grid size was increased more and more frequencies overlaped leading to a more complete picture of the DOS. A grid size bigger than 32x32x32 could have been used for the following calculations (e.g: 64x64x64) for more accuracy however it was more time consuming (13.5 seconds for 32x32x32 against 820 seconds for 64x64x64).
CaO presents a very similar structure to MgO. They are both in the same group, however Ca is slightly larger than Mg. Thus, CaO has a larger lattice structure and therefore a larger lattice parameter (a=0.480 Å) , a, than MgO (a=0.421 Å).[6] From this can be deduced that the same grid size could be used to run a calculation on CaO. On the other hand, Faujasite, a zeolite with covalent bonding character, has a much bigger lattice parameter (a = 24.638–24.65 Å) than MgO. According to equation 2, this means that the the reciprocal lattice parameter in much smaller than MgO. Therefore, a smaller grid size could be used for the simulations on Faujasite. Finally a metal, such as lithium, can be considered. The main aspect of Lithium is that it present a metallic bonding character.This leads to a low dispersion due to the 'sea of electron' surrounding the metal center, which acts as a shield. Thus a much smaller grid would be used for lithium.
Eqn. 2
Lattice Dynamics Calculation on MgO, primitive unit cell
In order to study the thermal expansion of MgO, the free energy of the system was computed at different temperatures, starting from 0K to 1000K in steps of 100K. The Lattice method was used here, estimated in the quasi-harmonic approximation. By plotting the graph of the Helmholtz free-energy against temperature, it was observed that the free energy decreases as the temperature increases (figure 5). This can easily be explain using the Helmholtz free-energy equation (equation 3). The entropy, S, which represents the disorder of the system, increases as the temperature escalates. Therefore, when the temperature was raised, the TS component led to an overall decrease of the Helmholtz free-energy. By looking at the graph for the lattice parameter against temperature, it was possible to observe an increase in a as the temperature was increased. The origin of the thermal expansion comes from the rise of the atomic vibrations. The thermal expansion occurs because of the asymmetric character of the interatomic potential (figure 8).[7] Indeed, as the temperature increases the atoms are able to access higher energy levels. On those energy levels, the equilibrium distance is greater than the original one, which leads to more vibrations and thus to an increase in the lattice vibrations of the crystal. This is then reflected in an increase of the lattice parameter. However, although very rare, some materials presents a negative thermal expansion, such as Cubic Zirconium Tungstate (ZrW2O8). This means that the materials contract with increasing temperature. Therefore, the normal potential energy diagram can no longer describe this phenomenon.[8] It is believe that it originates due to the transverse vibrations. [9] In a diatomic molecule with an exact harmonic potential, an increase in bond length would not be expected. As it was mentioned previously, the increase in bond length is due to the fact that the higher energy levels present a larger equilibrium internuclear separation. In the case of a purely harmonic potential, the equilibrium distance stays the same as we go on higher energy levels, which would explain why an increase in bond lenght is not expected as the temperature increases.
Figure 9: Volume against Temperature of the Lattice Dynamic simulation
From the plot showing the lattice parameter increasing with Temperature (Figure 6) can be determined the thermal expansion coefficient. Indeed, as the lattice parameter increases, the volume of the cell increases. The change in volume can be plotted in equation 1 to find . By looking at the graph above (figure 9), volume against temperature, was found to be 0.0005. Note that the first three points, 0 K, 100K and 200K, were not included in the the trendline as they do not follow the linear relationship. Values between 300 K and 1000 K follow the Quasi-Harmonic Approximation trend. However at much higher values, neglected anharmonic terms have a significant role in making the thermal expansion coefficient diverge.[10] Thus, as the temperature approaches the melting point of MgO, which is 3100 K, the phonon modes no longer represent the actual motion of the ions.[11] Therefore, the main limitation of Lattice Dynamics is that it cannot be used at high temperature. The thermal expansion coefficient,, was calculated by dividing the gradient by the initial volume, which was 18.7200 Å-3. was found to be 2.6709 x 10-5 K-1. This value is slightly smaller than the experimental one: =1.34 x 10 -5 K-1 run at room temperature.[12]
Molecular Dynamics Calculation on 32 MgO
For 32 MgO, a 2x2x2 supercell, a different type of calculation was used: Molecular Dynamics. The reason this type of simulation was not used to calculate the MgO primitive cell is that a larger number of atoms is required to obtain higher vibration degrees of freedom. In order to compare later the Lattice dynamic and Molecular dynamic simulations, a 32x32x32 grid was used. The simulations were computed on the same temperature range as the lattice dynamic simulations. By looking at the graph below (figure 10), volume against temperature, was found to be 0.0197. Note that the first two points, 100K and 200K, were again not included in the the trendline as they do not follow the linear relationship. Since Molecular Dynamics is based on classical mechanics equations that recorded the kinetic energy applied to the system, it was not possible to calculate 0 K as there is no kinetic energy at that point. The thermal expansion coefficient,, was calculated by dividing the gradient by the initial volume, which was 595.74 Å-3. was found to be 3.3068 x 10-5 K-1. Even though Molecular Dynamic simulations allow the use of larger models, its main limitation is that it cannot be used at low temperatures. Indeed, MD is based on equations of classical mechanics, which cannot be used at low temperature, where the quantum effects dominates. [13]
Figure 10: Graph of Volume against Temperature for the Molecular Dynamic simulations
Comparison between the two calculation methods
Figure 11: Graph of Volume against Temperature for the Lattice Dynamic and Molecular Dynamic simulations
In order to compare the volume expansion of the system between Lattice dynamic and Molecular dynamic simulations, the volume values found for the Molecular dynamic calculations were divided by 32, which is the number of MgO molecule in the supercell. For both type of simulation a normal thermal expansion was predicted. As it was stated previously, the experimental thermal expansion coefficient was found to be 1.39 x 10-5 K-1. This value is much smaller than the two computed values. This difference is due to approximations that were done on the system. The program DL visualize assumes that the crystal is perfect and infinite, in other words defects-free. This would probably increase the thermal expansion coefficient as the ionic bonds would be weaken. Because the simulations simplify the situation within the crystal, they tend to neglect the coloumbic interactions. This is reflected in the expansion of the system, that would be more difficult if the coloumbic interactions were taken into account. [14]
By comparing the LD and MD, it was determined that they work best at a different temperature range. Lattice Dynamics, with the quasi-harmonic approximation, presented good results at low temperature but gave inaccurate ones at very high temperature. This is due to the fact that as temperature increases, the anharmonic effects increases as well, which make them no longer negligible. Molecular Dynamics cannot be used at low temperature, where the quantum effect dominates. Indeed, Molecular Dynamics are based on classical mechanics equations, thus work better at high temperature.[15] This explains the difference in thermal expansion coefficient between the two methods. Upon comparison of the thermal expansion coefficients obtained, it was possible to determine that the Lattice dynamic calculation are more accurate than the Molecular dynamic calculations.
Conclusion
In this experiment, two types of simulations were compared in order to study the thermal expansion of MgO: Lattice Dynamics and Molecular Dynamics. It was determined that both values obtained were bigger than the one obtained experimentally. It was discovered that both types of simulations present limitations. On one hand, Lattice dynamics were useful to run calculations over a wide range of temperature but broke down at very high temperature. On the other hand, Molecular Dynamics worked better at high temperature, where the quantum effect is smaller, but turned out to give poor results at low temperature. Further experiments would study the effect of change in pressure with those two methods.
References
↑T.Albrecht, Electronic Structure of Solids, Lecture 1, 2017, 13-15
