Rep:Mod:JPS1121
The Free Energy and Thermal Expansion of MgO (James Simpson CID:00733493)
Introduction & Theory
In this computation numerous exercises will be done and analysed. Firstly the optimum grid size of an MgO lattice will be found in order to compute an accurate density of states function and to calculate the Helmholtz free energy of the system. Also in this experiment will be the calculation of the thermal expansion coefficient using the quasi-harmonic approximation and molecular dynamics. These two methods will then be analysed to literature values and compared to each other.
Lattice
Magnesium oxide is an solid made of an ionic lattice. Within this lattice the ions form a regular repeating pattern. The structure of MgO is face-centered cubic (Cubic close packed) and so contains four repeat units of MgO and therefore eight atoms overall in the conventional cell. This experiment uses this structural arrangement of oxide and magnesium ions and assumes that the lattice contains no defects. During the course of the computational exercises there will be interconversion between the lattice parameters of the primitive cell, which contains one oxide ion and one magnesium ion, and the volume of the primitive cell. The 3D shape of the primitive cell is rhombohedral. The volume of it is below, a is the cell parameter:
The coordinate system used in this experiment is not as expected direct space, but reciprocal space, or k-space. This is done by conducting a fourier transform on the spatial wavefunction of the original lattice. The unit of length in k-space is k which is related to the frequency of direct space. If vibrations between atoms in the lattice are considered then the value of k is related to the wavelength of the vibration as follows:
These vibrations can be considered as waves, as is usually done classically, but they can also considered to be particles in a similar way that light can be. These vibration particles are called phonons. A k-point is simply a point in k-space. There can for one k-point be multiple branches of phonon dispersions. This means that there are different energy phonons for a particular k-point. This is because in a diatomic lattice there are vibrational modes that in a monoatomic lattice were degenerate and are no longer. This leads to two vibrational energies where previously there was one. This of course can happen in each of the three spatial dimensions, leading to six branches in the lattice of MgO.
Free Energy
Gibbs Free Energy
In chemistry the Gibbs free energy is usually used to determine if a particular reaction will be spontaneous. Typically the quantity wanted is not the Gibbs free energy itself but the change in Gibbs free energy. It is a thermodynamic property which measures the amount of non-expansive that a closed system can do at constant temperature and pressure. The change of the Gibbs free energy ΔG, in a reversible reaction, equals the work exchanged by the system with its surroundings.
Helmholtz Free Energy
The Helmholtz free energy is a measure of the useful work that a closed system can do at a constant temperature. Like the Gibbs free energy it is a thermodynamic quantity. When the system in question is at constant volume then the Helmholtz free energy is exactly equal to the amount of work transferred or taken from the surroundings. This is because no work is performed on the surroundings. In terms of k-space the Helmholtz free energy is the summation of all phonons of every branch of every k-point.
Thermal Expansion
When a substance is heated it increases in volume. This is due to phonons interacting with the molecules in the substance. Temperature is a measure of average kinetic energy and so when temperature is increased the kinetic energy of the molecules increases, leading to an increase of vibrational quantum number. Therefore the bonds of the molecule contract and expand and the molecule moves up and down the Leonard Jones potential. In a purely harmonic oscillator the average bond position would not change. However the Leonard Jones potential is an anharmonic oscillator and so, as temperature increases the average bond length will increase also. The rate of this increase can be calculated and is shown below:
This quantity is known as the thermal expansion coefficient and as seen in the rate of change of substance, or cell, volume compared to change in temperature. This is then divided by the volume of the cell at absolute zero. The value of this constant is related to temperature and increases with increasing temperature.
Computational Programmes
DL Visualize is the programme that all the computations will be accessed with. It is a free to use package used for visualization of the properties and structures of materials and provides the user with the ability to make and see crystals, molecules and even surfaces. There are many tools within the programme that can be used to learn about the MgO but the one used for this exercise is GULP.
GULP is the programme that will actually do all the calculations within this experiment. The programme can perform many different types of simulation on substances using different numbers of spatial dimensions, in this case three. The programme's code focuses on analytical solutions, through the use of lattice dynamics where possible, rather than on molecular dynamics. However in the experiment both molecular dynamics and lattice dynamics (quasi-harmonic approximation) will also be used.
Techniques
The quasi-harmonic approximation is the main method used in this experiment. It is an improvement upon the harmonic approximation and allows for thermal expansion by using phonons. This is done by assuming that the harmonic approximation is correct at every bond length, or cell parameter length, so that from the point of view of phonons for each volume the harmonic approximation holds, so the energies of vibrations are in fact dependent on the volume. The other technique used is molecular dynamics, which uses Newton's Second Law to allow systems to evolve in time. This means that atoms simply follow the trajectories that they would have followed in reality and so properties have to be modelled as time averages of their behaviour. Using this technique involves using a larger number of particles than the quasi-harmonic approximation does and so the properties of such complex systems cannot be solved analytically. This is the reason that the properties must be time averaged; there are large errors involved with the calculations made using molecular dynamics.
Computational
Calculating the phonon modes of MgO
Phonon Dispersion
By using the Quasi-Harmonic Approximation, it is possible to show the vibrational modes of MgO. Figure 1 shows the phonon dispersion of the material. There are many vibrations that the atoms in the crystal can make: every k-point is associated with a vibration and therefore a wavelength and an energy as well, which can range from 0cm-1 to about 1200cm-1. There are two types of atom in the lattice, magnesium and oxygen, and each of the atom types can also vibrate in three spatial dimensions, leading to the six branches seen in the phonon dispersion graph for each k-point. The density of states function can be predicted from the phonon dispersion graph. The DoS shows the number of degenerate vibrational states at each energy level. At 400cm-1 there are many k-points that have vibrational modes around this energy level. This is why the DoS function, Figure 3, has many vibrational states around 400cm-1. Figure 1 also shows that there are no phonons with an energy above about 1200cm-1, and so the DoS shows that there are no states which have an energy above 1200cm-1
Density of States
The Density of States (DoS) function shows the possible vibrational modes and the degeneracy of these modes for a number of k-points. Figure 2 shows the DoS function for a 1x1x1 grid, meaning that only one k-point is sampled. The energies of the phonons are; 286cm-1, 351cm-1, 676cm-1 and 806cm-1. From looking at Figure 1 the k-point can only be either L or X as those are the only k-points for which four phonons are to be expected. There are six vibrations in total as there are six branches but two of these vibrations will be degenerate in energy. In order to determine if the k-point is L or X then the energies of the phonons must be looked at. As stated earlier the energies of the phonons are found from looking at the phonon dispersion graph (Figure 1) and the energies are; 286cm-1, 351cm-1, 676cm-1 and 806cm-1. With this information it is clear that the k-point is L. This is because the two phonons, which are close in energy, are 286cm-1 and 351cm-1: if the k-point used was X, then the two phonons that were close in energy would be found at around 500cm-1.
To compute a DoS function of use many more k-points must be sampled; there is an equilibrium that must be reached between the accuracy of the DoS function when compared to experiment and the time taken for the computation. It is impossible to compute a DoS function as found by experiment as this would require an infinite number of k-points and therefore infinite computational time. However, as the grid of k-points is increased more k-points are sampled, although not the number that perhaps might be expected. For a 2x2x2 grid eight k-points should be sampled, however only four have been sampled - this is because four k-points describe all the unique k-points in this grid. The other four k-points can be described as one of the four already sampled, due to the symmetry of k-space so, in order to use computation time more effectively, k-points that have already been sampled are not re-sampled. As shown in Figure 1, a 1x1x1 grid shows four phonons, a 2x2x2 grid gives seven non-degenerate phonons.
| Grid Size | Computation Time |
|---|---|
| 1x1x1 | 0.0024s |
| 2x2x2 | 0.0034s |
| 4x4x4 | 0.0113s |
| 8x8x8 | 0.0852s |
| 16x16x16 | 0.7606s |
| 32x32x32 | 11.4884s |
For smaller grid sizes, only isolated phonons energies can be seen. However, once the grid is large enough then a continuum of energies is found. Figure 3 shows a continuum formed by a 16x16x16 grid with the majority of modes can be seen to be around 400cm-1. The DoS generated from a 16x16x16 grid, was chosen as the best approximation. The time taken for the computation was much longer than that for a 8x8x8 grid, however the time taken for the computation is still small and will not inconvenience the person who is doing the computation. If a larger gird was chosen then the DoS function would be smoother but this would require a much longer computational time of 11.4884s. A 8x8x8 grid will also give a continuum of energies but as the extra time needed to compute a 16x16x16 grid is still small a 16x16x16 is the optimum grid size needed in order to generate a DoS function quickly.
Speculation
The grid determined as the most optimum could also be used for different lattices. For CaO this approach will most likely be valid1 as the calcium and magnesium ions are similar and both the oxides have the same lattice structure, face-centered cubic. The main difference between them is size and thus while effecting the phonon energies, the shape of the DoS should not change much. However, for a zeolite with a face-centered cubic space group2, for example Faujasite, the molecule contains many more atoms than MgO and so there will be many more branches and phonon dispersions will not appear the same as for MgO. Therefore DoS function will also look different: it is quite possible that as there are so many vibrational modes for each Faujasite molecule, a continuum of energies will be found with a relatively small grid size. However as phonon dispersions will look very different to MgO, the DoS function for this small grid might not accurately represent phonon dispersions for Faujasite. For a crystal made of only atom type, in this case Lithium, the grid would also not work. One of the reasons for this is that the crystal structure is not face-centered cubic, but body-centered cubic. The lattice also contains only one atom type and so there would only be three branches expected. The grid would probably need to be larger than that for MgO as metals typically need more k-points than semi-metals and insulators in order to have an accurate DoS function.
Calculating the Helmholtz Free Energy of MoG
Using the quasi-harmonic approximation the Helmholtz free energy can be computed. The Helmholtz free energy is simply the sum of all the vibrational modes for all of the branches for every k-point. However, just as with the DoS function computed earlier, in order to do this every k-point an infinite computational time would be required and thus a smaller grid must be used, just as before, in order to find the Helmholtz free energy for a required precision. The grid size was gradually increased and so values closer to the true values were found empirically. These results are found below:
| Grid Size | Helmholtz Free Energy (eV) | Difference in energies (eV) | Helmholtz Free Energy (kJmol-1) | Difference in energies (kJmol-1 |
|---|---|---|---|---|
| 1x1x1 | -40.9303(01) | - | -3949.1480(06) | - |
| 2x2x2 | -40.9266(09) | 0.0036(92) | -3948.7918(10) | 0.3561(96) |
| 3x3x3 | -40.9264(32) | 0.0001(77) | -3948.7747(13) | 0.0170(97) |
| 4x4x4 | -40.9264(50) | 0.0000(18) | -3948.7764(19) | 0.0017(06) |
| 5x5x5 | -40.9264(63) | -0.0000(13) | -3948.7777(49) | -0.0013(30) |
| 6x6x6 | -40.9264(71) | -0.0000(08) | -3948.7784(75) | -0.0007(16) |
| 7x7x7 | -40.9264(75) | -0.0000(04) | -3948.7788(82) | -0.0004(07) |
| 8x8x8 | -40.9264(78) | -0.0000(03) | -3948.7791(23) | -0.0002(41) |
| 10x10x10 | -40.9264(80) | -0.0000(02) | -3948.7793(73) | -0.0002(50) |
| 12x12x12 | -40.9264(81) | -0.0000(01) | -3948.7794(87) | -0.0001(14) |
| 14x14x14 | -40.9264(82) | -0.0000(01) | -3948.7795(49) | -0.0000(62) |
| 16x16x16 | -40.9264(82) | -0.0000(00) | -3948.7795(80) | -0.0000(31) |
| 32x32x32 | -40.9264(83) | -0.0000(01) | -3948.7796(41) | -0.0000(61) |
| 64x64x64 | -40.9264(83) | -0.0000(00) | -3948.7796(48) | -0.0000(07) |
Optimisation
It is clear that the calculated value for the Helmholtz free energy increases from a 1x1x1 grid to a 3x3x3 but then decreases for any grid size larger than 3x3x3. The difference in energies also decreases, suggesting that larger grids in fact calculate a free energy closer to the true value than a smaller grid. The largest grid used was a 64x64x64 grid: the energy difference between that grid and a 32x32x32 grid was only 0.000007kJmol-1. Increasing the grid size any further after that would be pointless as the computational time would be relatively large and the Helmholtz free energy calculated would be almost identical. To illustrate the trend the results were also plotted: see Figure 4.
Depending on the precision of the Helmholtz free energy needed, different grid sizes would be required. If the Helmholtz free energy was required to a precision of 1meV then a 3x3x3 grid would be the smallest grid size needed to calculate the free energy to the required standard of -40.9264. However, for a precision of 0.5meV then a grid size of 2x2x2 would be the smallest grid size required to calculate the free energy as requested. However this would not be reasonable as it would mean that a smaller grid size would be the minimum standard if a higher precision was wanted. Therefore to calculate the Helmholtz free energy to 0.5meV, a 3x3x3 grid should still be used. In the case of 0.1meV precision a 4x4x4 grid would be the minimum grid size needed as it rounds to -40.92648. The most precision that the Helmholtz free energy can be reported to would be 0.1meV, so the free energy can therefore be stated as -40.9265eV.
However there is another way to decide the smallest grid size needed for a certain precision i.e.by using the differences between a calculated Helmholtz free energy for a particular grid size and the Helmholtz free energy calculated with the grid size one smaller as it is known that the energies converge as the grid size increases. The largest energy difference that is smaller than 1meV is 0.000177eV, being the energy difference between a 2x2x2 grid and a 3x3x3 grid; this means that a 2x2x2 grid size is the smallest grid in which the calculated Helmholtz free energy is less than 1meV from the actual value. It is also accurate to 0.5meV while the smallest grid for which the Helmholtz free energy is less than 0.1meV from the true value is a 3x3x3 as the energy difference between it and a 4x4x4 grid is 0.000018eV.
Speculation
As in the previous computational exercise it may be that the optimal grid size required for the most precision possible is also the optimal grid size for other molecules. The Helmholtz free energy is the summation of all vibrational modes for all of branches for every k-point so the the size of the grid needed for CaO, Faujasite and Lithium has to be smaller than those in 3.1.3. Once an accurate DoS function has been found it would be expected that the Hemlholtz free energy calculated would be accurate. However for MgO the grid size needed for an accurate and precise Helmholtz is smaller than that needed for an accurate DoS function as it would be expected that the relationship of the size of the grids compared with MgO will be the same as in the previous section. However, the grids will be smaller than those needed in part one for an accurate DoS function.
The Thermal Expansion of MgO
The quasi-harmonic approximation can be used to work out the Gibbs free energy. This is done by optimising the structure of the primitive cell of MgO and therfore can be done for mulitple termperatures. The values for the free energy and the cell volume, along with the parameters of the primitive cell, at each temperature are found below:
| Temperature (K) | Cell Volume (Å3) | Cell Parameter (Å) | Gibbs Free Energy (eV) |
|---|---|---|---|
| 0 | 18.8365 | 2.9866 | -40.9019 |
| 100 | 18.8365 | 2.9866 | -40.9024 |
| 200 | 18.8562 | 2.9876 | -40.9094 |
| 300 | 18.8901 | 2.9894 | -40.9281 |
| 400 | 18.9325 | 2.9916 | -40.9586 |
| 500 | 18.9800 | 2.9941 | -40.9994 |
| 600 | 19.0312 | 2.9968 | -41.0493 |
| 700 | 19.0851 | 2.9996 | -41.1071 |
| 800 | 19.1413 | 3.0026 | -41.1719 |
| 900 | 19.1997 | 3.0056 | -41.2430 |
| 1000 | 19.2601 | 3.0088 | -41.3198 |
| 1250 | 19.4205 | 3.0171 | -41.5339 |
| 1500 | 19.5953 | 3.0261 | -41.7749 |
| 1750 | 19.7878 | 3.0360 | -42.0391 |
| 2000 | 20.0032 | 3.0470 | -42.3237 |
| 2250 | 20.2504 | 3.0595 | -42.6268 |
| 2500 | 20.5478 | 3.0746 | -42.9472 |
| 2750 | 20.9451 | 3.0941 | -43.2843 |
| 2900 | 21.3230 | 3.1126 | -43.4948 |
| 2950 | 21.5522 | 3.1237 | -43.5667 |
Calculation
After this temperature the Quasi-Harmonic Apporximation breaks down and so the computation can no longer accurately be done. This is because the melting point of MgO, 3100K3 , is being approached. From Figures 5 and 6 it is clear that there is a quadratic relationship between increasing temperature and the lattice parameter. There is also an inverse quadratic relationship between the Gibbs free energy and an increase in temperature. This however means that it is harder to calculate the thermal expansion coefficient for MgO as the coeficient is not just simply the gradient of a straight line. From inspection of the graphs it is noticeable however that the data points of the lattice parameter between 400K and 1000K can be approximated to be linear, so that an estimate of the thermal expansion coefficient can be calculated. Doing so gives the gradient of the line, using the volume of the primitive cell instead of the cell parameter, to be 5.602x10-4 Å3K-1. If this is divided by the primitive cell volume at absolute zero, 18.8365Å3, then the thermal expansion coefficient of MgO is found to be 2.9740x10-5 K-1. In comparison a literature value of two points in this range are 1.13x10-5 K-1 at about 525K to 1.31x10-5 K-1 at about 925K, so the calculated value is about twice as expected.
A more accurate measurement therefore would be to use the quadratic trendline seen on Figure 5 in order to calculate the thermal expansion coefficient. The equation of the trendline when using cell volumes and not cell parameters is 2x10-7x2+10-4x+18.846. This gives the equation of the gradient to be 2x10-8x+9x10-6 and from this the gradient and therefore the thermal expansion coefficient of MgO made without a linear approximation can be made. At 525K the gradient is 3.1x10-4 Å3K-1 and this gives the value of the thermal expansion coefficient at 525K to be 1.6457x10-5 K-1. This value is much closer to the accepted value of 1.13x10-5 K-1 4. At the the other value for which literature values have been taken, the gradient is 2.85x10-4 Å3K-1. This gives the thermal expansion coefficient at 925K of 2.4952x10-5 K-1 and the literature value found at this temperature is 1.31x10-5 K-1. While these calculated values are still inaccurate they are closer to the actual values, showing that a linear approximation decreases the accuracy of the calculated thermal expansion coefficient.
Speculation
There is still one question to be answered i.e. why does MgO expand when heated? As the oxygen and magnesium ions in the lattice increase in energy they vibrate more and more. In a simple harmonic oscillator this would have no overall effect as although the ions would travel further up the potential energy surface the equilibrium bond length would not change. However in the quasi-harmonic approximation the equilibrium bond length can change, this is because for every bond length the harmonic approximation is assumed to hold. Therefore as the temperature increases, by following multiple harmonic oscillator paths, the equilibrium bond length also increases. When the temperature is very high, meaning approaching the melting point, the approximation breaks down. This is because the quasi-harmonic approximation does not take phonon- phonon interactions into account and these interactions occur at higher temperatures.
Molecular Dynamics
The final computational exercise involves calculating the thermal expansion coefficient using a different method to the quasi-harmonic approximation. The method used is molecular dynamics and in order to compute an accurate value for the thermal expansion coefficient a cell larger than the primitive cell must be used. The cell used is a 2x2x2 super cell, so-called because the super cell is eight times larger than the MgO conventional cell and contains 32 MgO units. The results are below:
| Temperature (K) | Cell Volume (Å3) | Primitive Cell Volume (Å3) |
|---|---|---|
| 100 | 599.3687 | 18.7302 |
| 200 | 600.8710 | 18.7772 |
| 300 | 602.8994 | 18.8406 |
| 400 | 603.5871 | 18.8621 |
| 500 | 605.9801 | 18.9369 |
| 750 | 610.2941 | 19.0717 |
| 1000 | 613.9842 | 19.1870 |
| 1250 | 618.2643 | 19.3208 |
| 1500 | 623.0140 | 19.4692 |
| 1750 | 626.0555 | 19.5642 |
| 2000 | 636.3543 | 19.8861 |
| 2250 | 638.9943 | 19.9686 |
| 2500 | 651.6952 | 20.3655 |
| 2750 | 654.9254 | 20.4664 |
Calculation
It is clear from Figure 7 that the relationship between the volume of the primitive cell can be modelled as either linear or quadratic. The linear case will be considered first. In the previous section it was seen that the literature values of the coefficient change depending on the temperature. Moreover, when a linear approximation was used for cell volumes by using the quasi-harmonic approximation it was seen that the calculated coefficient was less accurate than when the quadratic relationship was used. However molecular dynamics is a different approach and so must be considered: the equation of the line of best fit through this data is 0.0006x + 18.598. Molecular dynamics cannot be carried out on the MgO system at 0K and so the value for Vo was taken to be the y-intercept of the equation as this is the expected volume if the computation was able to be carried out at 0K. The coefficient using this data was calculated to be, for all temperatures, 3.2262x10-5 K-1. This value, while the same order of magnitude as literature, is much less accurate than even the linear approximation of the quasi-harmonic approximation.
However for the quadratic case the equation of the line of best fit is 10-7x2 + 0.0003x + 18.721. This gives the gradient to be 2x10-7x + 0.0003, giving the value of thermal expansion coefficient to be 2.1633x10-5 K-1 at 525K and 2.5907x10-5 K-1 at 925K. The values calculated using the quasi-harmonic approximation were 1.6457x10-5 K-1 at 525K and at 925K the value was calculated to be 2.4952x10-5 K-1. The literature values found are 1.13x10-5 K-1 at 525K and 1.31x10-5 K-1 at 925K. The calculated values can now be compared. Only the values calculated using a quadratic fit will be compared as these values are more accurate than those calculated using a linear approximation.
Comparision
| Temperature (K) | Literature Value 10-5 K-1 | Quasi-Harmonic Value 10-5 K-1 | Molecular Dynamics Value 10-5 K-1 |
|---|---|---|---|
| 525 | 1.13 | 1.6457 | 2.1633 |
| 925 | 1.31 | 2.4952 | 2.5907 |
It is clear from the table above that the Quasi-Harmonic approximation has given thermal expansion coefficient values that are closer to the literature values than those calculated using molecular dynamics. However the temperature used to calculate the thermal expansion coefficient only goes to about a third of the range of temperatures for which MgO is a solid. It can also be seen that the difference in the thermal expansion coefficient values for the quasi-harmonic approximation is about double the difference for the values calculated using molecular dynamics. This means that as the temperature increases the values calculated using the quasi-harmonic approximation increase more than in the molecular dynamics calculations. Therefore it can be predicted that at a sufficiently high temperature the thermal expansion coefficients calculated using molecular dynamics will be more accurate than those calculated using the quasi-harmonic approximation. One of the major reasons for this is that the programme used for the computational exercise, GULP, focuses on analytical solutions using the quasi-harmonic approximation. However properties that are being measured using molecular dynamics cannot be solved using analytical solutions. Therefore the using molecular dynamics itself for anything involving GULP will inherently inaccurate.
Speculation
It is also possible to calculate the free energy using molecular dynamics instead of using the quasi-harmonic approximation as was done previously. When using quasi-harmonic approximation to do this a 2x2x2 grid size was required but it is known that a larger grid is needed when using molecular dynamics. As a 2x2x2 grid was needed for the thermal expansion coefficient a 2x2x2 super cell should be used for molecular dynamics, being being a 4x4x4 cell relative to the conventional MgO cell. At high temperatures, for instance the melting point, both methods should break down - this is seen with the quasi-harmonic approximation at high temperatures around 2900K, as when the computation is run at higher temperatures the simulation fails. It is difficult to say what would happen using molecular dynamics. The data points show that the relationship is either linear or quadratic: if the relationship is linear then molecular dynamics will never predict a melting point, but if it is a quadratic relatonship, which is more logical, then with an increase in temperature the increase in the thermal expansion coefficient will increase, and so eventually there will be a melting point. However it is clear by 2750K that the increase in thermal expansion coefficient is small enough that the melting point predicted by molecular dynamics is much higher than the actual melting point.
Conclusion
In this computation numerous exercises were conducted and analysed. The optimum grid size of MgO was found in order to compute an accurate density of states function, 16x16x16, and also to calculate the Helmholtz free energy of the system. The grid size needed depended on whether precision or accuracy was required. For precision the grid size wanted was a 3x3x3 grid for 0.5meV precision and a 4x4x4 grid was needed to get a 0.1meV precision. If accuracy was wanted then a 2x2x2 grid was required to be within 1meV of the actual Helmholtz free energy and a 3x3x3 grid needed to be within 0.1meV of the true value.
The other major part of this experiment was the comparison of the thermal expansion coefficient for MgO using both the quasi-harmonic approximation and molecular dynamics and both of these methods were also compared to literature. It was found that using a linear approximation of the result at lower temperature made the values of the coefficients less accurate so that most accurate approach was using a quadratic fit. For both temperatures for which numerical values of the coefficent were calculated, the quasi-harmonic approximation gave the result closest to literature. However these results were still not close to those found in literature.
References
1. M.J.L. Sangster, T.G. Peckham, D.H. Saunderson, J. Phys. C., 1970, 3, 1026-1036
2. IZA Commission on Natural Zeolites, http://www.iza-online.org/natural/Datasheets/Faujasite/faujasite.htm, 03.11.2014
3. L. Vofiadlo, G.D. Price, Phys. Chem. Miner., 1996, 23, 42-49
4. J.B. Austin, J. Am. Ceram. Soc, 1931, 14, 798