Author Archives: Bryan William Brasile

Calculating the Barrier of Vinyl Alcohol Tautomerization to Acetaldehyde

Leave a reply

Tautomerization of Vinyl Alcohol

In this post we wish to determine the reaction energy, barrier, and transition state structure for the tautomerization reaction of vinyl alcohol converting to acetaldehyde using DFT. This reaction is shown below in Fig.1. Tautomers are constitutional isomers; having the same number of each atom, but with a different connectivity. Most often, as in this case, tautomerization reactions are intramolecular proton exchange reactions.

Fig. 1 Tautomerization reaction of vinyl alcohol (left) to acetaldehyde (right).

It is known that acetaldehyde is more stable than vinyl alcohol at ambient conditions [1]. Thus, we should expect to find that the reaction energy is negative (forward direction of Fig. 1), which we can verify after we have performed our calculations.

Setting Up the Calculation

We will use DMol3 to geometry optimize the reactant and product state as well as locate the transition state. In addition, we will calculate the vibrational modes and frequencies of each the reactant, product, and transition state for reasons we will return to later.

In order to begin a transition state search in DMol3, the user must supply initial geometries for the reactant and product. With DMol3, the transition state procedure includes a routine for optimizing the reactant and product states before finding the TS, so these initial structures don’t need to be very accurate. Below are figures showing the initial reactant and product structures used for this example.

Fig. 2 Initial geometries of reactant (left) and product (right) to DMol3.

As TS searches follow atoms through spatial translations, it is important to tell DMol3 which atoms from the reactant and product are “equivalent”. This can be done using the reaction preview tool. Once a reaction preview has been made between the initial and final states, we can enter the calculation parameters we wish to use.

Calculation Details

Calculations were performed using the DMol3 Module as available in Materials Studio. The TS search method used was the Complete LST/QST Method with conjugate gradient optimization of reactants and products using a RMS convergence of 0.00185 Ha/Å and a maxinum of 5 QST steps. Exchange and correlation were treated within the generalized gradient approximation using the PBE functional with the DNP basis set and explicit core electrons (all electron calculation) [2-5].  The SCF tolerance was set to 1.0E-5 Ha with linear smearing using a value of 0.005 Ha.

Results and Discussion

First we present the optimized structures and compare the geometries to experimental values available from NIST [6,7].

Fig. 3 Optimized reactant and product from DMol3.

Table 2 Calculated and Experimental Geometries of Vinyl Alcohol
CalculatedExperimental
C1-C21.3361.326
C2-O31.3711.372
C2-H71.0911.097
C1-H51.0881.079
C1-H61.0931.086
O3-H40.9740.960
C1-C2-O3126.7126.2
C1-C2-H7122.9129.1
C2-C1-H6122.0121.7
C2-C1-H5119.9119.5
C2-O3-H4108.6108.3
CalculatedExperimental
C1-C21.5021.501
C2-O31.2181.216
C1-H41.0971.086
C2-H71.1211.114
C1-C2-O3124.7123.9
H4-C1-H6110.3108.3
C1-C2-H7114.9117.5

Below is a plot showing the evolution of the TS search. This plot includes the LST, QST, and GC (conjugate gradient) paths along the reaction coordinate. Note the transition state is also labeled on this plot–indicating a transition state has been found with the calculation parameters.

Fig. 4 DMol3 TS Search Energy Trajectories

DMol3 also outputs the geometry of the TS. One should confirm that the TS found (by computational means) matches intuition about how the reaction proceeds. Below the TS is shown with atom labels.

Fig. 5 TS structure found using DMol3.

The transition state structure is what we might expect: somewhere in between vinyl alcohol and acetaldehyde. Here we can see that at the TS, the C-C double bond hasn’t quite broken while the O-H bond hasn’t quite formed. This leads us to believe we may have found the correct transition state, however, we must verify by examining the vibrational modes/frequencies of the TS. Transition states occur at saddle points in the potential energy surface (PES) and, therefore, should have one negative (imaginary) frequency corresponding the the vibrational mode along the reaction path.

The vibrational analysis tool can be used on the TS to verify these requirements are met. Indeed, we find the TS has one negative vibrational frequency: -2000. To verify this mode is along the reaction path we may animate it and confirm visually. This animation is reproduced below as a series of frames.

Fig. 6 Imaginary mode for the TS found using DMol3.

Thus we can be confident that we have found a real TS corresponding to the tautomerization of vinyl alcohol. Of special importance are the reaction energy and reaction barrier which DMol3 calculates during the search procedure. The energy for this reaction was -11.2 kcal/mol while the barrier was 51.9 kcal/mol. Note that our prediction of a negative reaction energy was correct. We also notice this reaction has a relatively high barrier (inaccessible at room temperature).

So, we have found a TS structure for the reaction of vinyl alcohol to acetaldehyde and have calculated the reaction energy and barrier. Next we need to compare the structure and energies with literature values.

Juan Andres et. al. did a comprehensive ab inito study on this very reaction and their results are readily available online [8]. In their results they report that DFT methods found the reaction energy to be in the range 8.5-14.3 kcal/mol while the barrier was found to fall within the range 87.9-93.5 kcal/mol. Thus there is great agreement with our reaction energy but the barriers are in severe disagreement. This, however, is most likely a symptom of basis set selection and the differences in DFT methods for finding TS between our attempt and those of Andres.

To verify that our structures are the same we reproduce their table of geometry data and include our own. The labels in the tables below are consisted with Fig. 4.

Table 1 TS Geometry
This WorkLiterature
C1-C21.407
1.428
C1-O32.2282.299
C1-H41.5001.541
O3-H41.3111.305
C1-C2-O3111.093110.500
H4-C1-C266.73767.200
C2-O3-H475.82778.600
C1-C2-H7130.421131.100
H5-C1-C2-O3155.067147.4
H6-C1-C2-O363.97677.400

From the table above we see there is great agreement in the TS geometry, with the greatest errors showing in the dihedral angles. This leads us to conclude that we have indeed found a TS and that the departures from literature values (especially the reaction barrier) are purely from differences in the level of theory/basis set used in the calculation.

References

[1]      R.D. Johnson III. “CCCBDB NIST Standard Reference Database”. Retrieved 2014-08-30.

[2]     B. Delley, J. Chem. Phys. 92, 508 (1990).

[3]     B. Delley, J. Chem. Phys. 113, 7756 (2000).

[4]     J. P. Perdew, K. Burke, and M. Ernzerhof. Generalized gradient approximation                                 made simple. Phys. Rev. Lett., 77:3865, 1996.

[5]     J. P. Perdew, K. Burke, and M. Ernzerhof. Erratum: Generalized gradient                                             approximation made simple. Phys. Rev. Lett., 78:1396, 1997.

[6]    Kuchitsu (ed.), Landolt-Bornstein: Group II: Atomic and Molecular Physics Volume 15:                     Structure Data of Free Polyatomic Molecules. Springer-Verlag, Berlin, 1987.

[7]    H. Hollenstien, Hs. H. Gunthard, “Solid State and gas infrared spectra and normal                           coordinate analysis of 5 isotopic species of acetaldehyde” Spec. Acta 27A, 2027 (1971)

[8]      Andrés, J. , Domingo, L. R., Picher, M. T. and Safont, V. S. (1998), Comparative                                   theoretical study of transition structures, barrier heights, and reaction energies for the                 intramolecular tautomerization in acetaldehyde/vinyl alcohol and                                                       acetaldimine/vinylamine systems. Int. J. Quantum Chem., 66: 9-24.

 

Considerations for DFT Frequency Calculations

Leave a reply

Using DFT to Calculate Vibrational Modes/Frequencies

In this post we explore the effect of numerical method as well as displacement magnitude on the results of vibrational mode calculations using DMol3.

The vibrational modes of a system or collection of atoms/molecules are determined from the Hessian Matrix, H, which has elements [1]:

\begin{equation}H_{ij}=\frac{1}{\sqrt{m_im_j}}\frac{{\partial^2{E}}}{\partial{q_i}\partial{q_j}}\end{equation}

In this notation qi are the coordinates of atom i. The eigenvectors of H are the vibrational modes themselves (displacements of atoms along the normal mode) while the vibrational frequencies are given as the square roots of the eigenvalues.

Numerical Methods

Most DFT packages (DMol3 included) cannot calculate an analytic function for the energy, E with respect to all of the degrees of freedom (coordinates). Thus, numerical approximations to the derivatives in the above expression need to be used.

DMol3 offers two options for the numerical method used in estimating such derivatives; these are the Single-Point Difference (Forward Difference) Method and the Two-Point Difference (Central Difference) Method.

The Single-Point expression is [1]:

\begin{equation}\frac{{\partial^2{E}}}{\partial{q_i}\partial{q_j}} \approx \frac{\frac{\partial{}}{\partial{q_i}}E(q_i+\delta)-\frac{\partial{}}{\partial{q_i}}E(q_i)}{\delta}\end{equation}

The Two-Point expression is [1]:

\begin{equation}\frac{{\partial^2{E}}}{\partial{q_i}\partial{q_j}} \approx \frac{\frac{\partial{}}{\partial{q_i}}E(q_i+\delta)-\frac{\partial{}}{\partial{q_i}}E(q_i-\delta)}{2\delta}\end{equation}

In these expressions, δ is the displacement length (magnitude) which is a calculation parameter that can be adjusted by the user.

While the documentation for DMol3 suggests the Central Difference method may reduce rounding-errors (while doubling the computational load), it gives little explanation as to any major differences in results obtained by each method.

Similarly, other than mentioning using smaller (larger) values of δ for steeper (flatter) potential energy surfaces, little insight is given as to the effect of the displacement length on the results of these vibrational mode calculations.

Thus, in this post we will use DMol3 to calculate the vibrational modes and frequencies of methane (CH4) using both the Single-Point (SP) method and the Two-Point (TP) method as well as the effect of the magnitude of displacement.

Vibrational Modes of Methane

As methane is a non-linear molecule with 5 atoms, there are 9 vibrational modes (3N-6, N=5), however; due to symmetry (Td point group for methane) there are only 4 unique vibrational frequencies. For reference, these modes are visualized below with their frequencies taken from the NIST Webbook [2].

Fig. 1: A1 Symmetry Mode, f=2917 1/cm

Fig. 2: E Symmetry Mode, f=1534 1/cm

Fig. 3: F2 Symmetric Stretch Mode, f=3019 1/cm

Fig. 4: F2 Symmetric Bend Mode, f=1306 1/cm

 

Calculation Details

All calculations (geometry optimization, frequency analysis) were performed using the DMol3 module as implemented in Materials Studio using the PBE functional for exchange and correlation with the DNP+ basis set [3-6]. SCF cycles were converged to 1.0-5 eV and the forces on the atoms were converged to 0.001 eV/Å. Core electrons were treated explicitly (all electron).

Results & Discussion

Before calculating the vibrational modes, we must first geometry optimize the methane molecule. This was done in DMol3 using the DNP+ basis set (for non-periodic systems). Below are images of the initial configuration as well as the optimized structure.

Fig. 5 Initial methane configuration

Fig. 6 Optimized methane structure

Below are plots of the various vibrational frequencies determined via the two methods for various δ (frequency as a function of 1/ δ is used for clarity). The frequencies (f#) are grouped according to symmetry (e.g. frequencies 1, 2, and 3 should be the same so they appear on the same plot). The choice of “which frequencies are identical” is  The displacement values used were 0.005, 0.0075, 0.01, 0.015, 0.02, 0.025, 0.05, 0.1, 0.25, and 0.5 Å.

Fig. 7: Vibrational Frequencies for Single Point Method various displacements.

Fig. 8: Vibrational Frequencies for Two Point Method at various displacements.

In the table below, the frequency for each mode in each method are listed. The frequencies were determined by averaging the frequencies for which the change in the frequency with displacement was ~10 cm-1 (the level of precision in acceptable experimental results). In the case of the both the SP and TP methods, all of the frequencies were converged for δ < 0.02 Å, though individual modes may have converged more quickly.

Table 1: Vibrational Frequencies for SP and TP Methods
ModeSP Frequency (1/cm)TP Frequency (1/cm)
11311.741315.76
21326.901315.76
31326.901315.76
41531.281535.83
51549.901535.83
62990.582996.03
73108.943116.29
83112.343116.29
93112.363116.29

First, it should be noted that while both the SP and TP predict 9 vibrational modes, the single-point method does not quite capture the reality of the situation. As mentioned earlier, due to belonging to the Td point group, there should be 4 unique modes. Using the numbering above, frequencies 1, 2, and 3 should be identical as should the pair 4, 5 and the trio of frequencies 7, 8, and 9. These modes are equivalent for all calculations done using the two-point method, as is observed by the overlapping of the points in Fig. 8. The SP method eventually converges to have approximately equal frequencies for the symmetric modes though using larger displacement values increases this slight inaccuracy.

In these cases, the single-point method can either predict a larger or smaller frequency than the two-point method with the former having some variation in the convergence patterns. This highlights how drastically calculation parameters can affect results and conclusions drawn from such calculations.

From the data above, we notice that for either the SP or TP method, as the displacement magnitude is lowered, the frequency of each mode eventually “converges” to a stable value. This is to be expected as smaller displacements are closer to the equilibrium bond lengths and the harmonic approximation is valid, while at large separations there are sever deviations from the simple harmonic approximation.

It is worth noting that had we used looser convergence criteria (or even a smaller basis set) we might expect the behavior of the frequency for very small displacements to be different. This is due to the fact that eventually the displacement size, and thus the changes in energy, would be comparable to or less than the machine precision possible using the “looser” calculation parameters. That is to say, using a smaller basis set or less demanding criteria for the energy and force convergence, the energy differences calculated numerically using the SP and TP methods would be non-significant and could produce unphysical results.

Given these arguments, our calculations using TP and SP methods are still in decent agreement, however, the above exercise is insightful in that it shows that convergence criteria for frequency calculations can have a large impact on the results, especially when using larger displacements or for larger molecules having more vibrational modes.

The most important conclusion we can draw from this work is that the single-point and two-point methods differ mostly in their accuracy of predicting symmetric modes, with the TP method outperforming the SP method. This suggests that when calculating the vibrational modes and frequencies of high-symmetry molecules, the TP method is better suited.

Finally, we can see that the calculated results for both are in decent agreement with the experimental values but with errors around 5%- for sufficiently small displacements. This may serve as a reminder that computation and simulation can help corroborate experimental findings but generally cannot predict exact values.

References

[1]     Accelrys Software Inc., Materials Studio Release Notes, Release 7.0, San Diego: Accelrys                      Software Inc., 2014.

[2]     Shimanouchi, T., Tables of Molecular Vibrational Frequencies Consolidated Volume I,                      National Bureau of Standards, 1972, 1-160.

[3]     B. Delley, J. Chem. Phys. 92, 508 (1990).

[4]     B. Delley, J. Chem. Phys. 113, 7756 (2000).

[5]     J. P. Perdew, K. Burke, and M. Ernzerhof. Generalized gradient approximation made                        simple. Phys. Rev. Lett., 77:3865, 1996.

[6]     J. P. Perdew, K. Burke, and M. Ernzerhof. Erratum: Generalized gradient approximation                 made simple. Phys. Rev. Lett., 78:1396, 1997.

 

Effect of Coverage on Adsorption Energies

Leave a reply

Binding Energies on Surfaces

DFT is routinely used to determine the adsorption energies of different atoms and molecules on metal surfaces. The adsorption energy is simply the change in energy when an atom or molecule is brought from (infinitely) far away from a surface to it’s equilibrium adsorption configuration.

In the case of a single atom X (of a diatomic molecule X2) adsorbing on a surface it is typical to evaluate the adsorption energy as [1]:

\begin{equation}E_{ads} = E_{surf+X} –  0.5E_{X_{2}}-E_{surf}\end{equation}

Adsorption Sites on FCC (111) Metal Surfaces

The (111) surface of FCC crystals (metals) have 4 unique adsorption sites. They are called the atop, fcc hollow, hcp hollow, and bridge site for the surface. They are pictured below convenience.

Fig. 1 Adsorption sites on Pt(111); (left to right) bare surface, atop, fcc, hcp, bridge).

Atoms and molecules tend to preferentially bind to certain sites. This is something we would like to be able to determine. Luckily, (1) is valid for all of these binding sites so we can simply calculate the adsorption energies directly to determine what site an atom/molecule of interest may bind to.

Coverage Effects

When using plane-wave DFT codes, ones must always be aware of mirror images interacting. The distance between mirror images in such DFT calculations depends on the size of the supercell chosen as well as the number of adsorbates in the supercell.

It is common practice to place only 1 adsorbate within a supercell, meaning the supercell size determines the distance between mirror images. Below are figures showing how mirror images of adsorbates might “see” one another and how the distance between mirror images changes with supercell size.

Fig 2. Two 2×2 supercells of an O atom adsorbed on Pt(111). (The indicated distance is in Angstrom)

Fig 3. Two 3×3 supercells of an O atom adsorbed on Pt(111). (The indicated distance is in Angstrom).

The above figures help us infer that adsorption energies might decrease (adsorption more favorable) with increasing supercell size. This is consistent with the idea that most interactions between atoms fall off pretty rapidly with distance (e.g. vdW).

We investigate this by comparing the adsorption energies of an O atom adsorbed on the atop, fcc, hcp, and bridge sites of  Pt(111) using two different sized supercells, (2×2) and (3×3). These correspond to 1/4=0.25 ML coverage (O:Pt = 1:4) and 1/9 = 0.11 ML coverage (O:Pt = 1:9) respectively.

Calculation Details

All calculations were performed using the plane-wave Vienna Ab Initio Software Package (VASP) with the PBE exchange-correlation functional [1-4,7-8]. Core electrons were treated using the Projector Augmented Wave approach [5,6]. 1x1x1 Monkhorst-Pack mesh was used to sample k-space for the isolated O atom whereas 12x12x1 and 8x8x1 Monkhorst-Pack meshes were used to sample k-space for the 2×2 and 3×3 supercells, respectively. The plane wave cut-off energy was set to 550 eV and the structural optimizations considered complete when the magnitude of the forces on each atom was less than 0.02 eV. Dipole corrections were included in all surface calculations. Surface calculations used a 4-layer slab model wherein the bottom two layers were frozen during optimization. The lattice constants used was that determined using DFT instead of experiment; a = 2.78.

For discussions on convergence with respect to k-points and energy cut-off follow this link.

Results

Below are presented all energies calculated using VASP for the purposes of this exercise. First we present the energies of the bare surfaces as well as the isolated oxygen molecule followed by the calculated energies of the adsorbed oxygen at different binding sites.

SystemEnergy (eV/atom)
O2-9.864
(2x2) Surf-5.762
(3x3) Surf-5.762
Table 1. Energies of Oxygen and Bare Surfaces

Fig. 4 Adsorption energies of atomic Oxygen on Pt(111) at 0.25 and 0.11 ML.

From the above results, particularly Fig.4 , a few observations can be made (elaboration to these observations is given in the following section):

  1. At 0.25 ML coverage, the adsorption energies for O at the fcc and bridge sites are identical, and the lowest out of all sites (meaning O appears to preferentially bind to both the fcc and bridge sites at this coverage).
  2. At 0.11 ML coverage the adsorption ebergies for O at the fcc and bridge sites are identical, and the lowest out of all sites (meaning O appears to preferentially bind to both the fcc and bridge sites at this coverage).
  3. The difference in adsorption energies between 0.25 and 0.11 ML coverage is somewhat inconsistent: ignoring the bridge site calculations and comparing only hcp, fcc and atop sites, the adsorption energy is slightly lower for 0.25 ML in the case of  the fcc site while lower for 0.11 ML in the case of the hcp site and again lower for 0.11 ML in the case of the atop site.
  4. Overall there is little difference in the magnitude of the adsorption energy between coverages for the same site.

Conclusions

We now try to reason reason with our results/observations from above.

In regards to point 1, a simple look at the optimized geometries reveals that initially placed bridge oxygen “fell” into the more stable fcc site. If one is careful about the choice in calculation parameters (specifically the maximum ionic displacement), it is possible to recover the actual bridge site adsorption energy. We leave this discussion here as it is theoretically and experimentally predicted that O will not bind to bridge sites, though calculating the bridge site adsoption energy would make a nice exercise in understanding how different calculation parameters affect one’s results. We can conclude that at 0.25 ML O adsorbs at the fcc site.

Point 2, similar to point 1, simply reveals that the 0.11 ML bridge site calculation “fell” to the more stable fcc, reinforcing the idea that the bridge site equilibrium geometry is sensitive to the calculation parameters. Regardless, the bridge sites for both 0.25 and 0.11 ML coverage failed to converge to the desired geometry and instead relaxed to other adsorption sites. From this we may conclude that the bridge site is not the preferred binding site of atomic oxygen on Pt(111).

Below we show the geometry of the system as built as well as after convergence for the bridge site calculation at 0.11 ML to show what we mean be the Oxygn “falling” into the more stable fcc adsorption site.

Fig. 5 0.11 ML bridge site geometry as built (left) and after convergence (right).

Comparing the two figures we can see the “guessed” (initial) position of the oxygen is rather close to both the equilibrium fcc and hcp sites. Due to this, during optimization when the atoms move, it is possible the O explores a region in space that is “too close” to the minimum associated with the fcc site. Since VASP is searching for a minimum and not a specific minimum, this means once the O explores regions of space that fall in the fcc minimum, the calculation will continue to allow the O atom to relax into the fcc site.

Moving to point 2 we may concisely conclude that at 0.11 ML O adsorbs at the fcc site.

Our last two observations are a little more nuanced. In this calculation scheme, we have chosen not to apply zero-point energy corrections, nor have we included any entropic effects. In this case, we might expect out entropic effects to be roughly the same since both the 0.25 and 0.11 ML cases have 1 O atom and the same number of metal atoms “before” and after “adsorption”. Zero-point energy corrections (ZPE) can be significant, at least relative the the energies we are considering.

Thus, for now, we can conclude that O binds to the fcc site of Pt(111) surfaces at 0.25 and 0.11 ML coverage and that without ZPE and entropy corrections, there is a negligible difference in adsorption energies with respect to the coverage.

While the effect of coverage is not clear using the methods outlines above, it is reassuring that our calculations predict O to preferentially bind to the fcc site at both 0.25 and 0.11 ML coverage, as is found in experiment and predicted by computation [9].

Future Work

Naturally one would like to see how the ZPE and entropy corrections influence the results. One would expect there to be some increase in the difference between adsorption energies at 0.25 and 0.11 ML coverage. On a similar note, in this work we have used an asymmetric slab model with 4 layers (as is common for efficiency). One may also consider a symmetric slab model or perhaps a 5 layer slab to see if there is any splitting between these two coverages. This will be explored in future posts.

References

 

[1]     G. Kresse and J. Hafner. Ab initio molecular dynamics for liquid metals. Phys. Rev. B,                      47:558, 1993.

[2]     G. Kresse and J. Hafner. Ab initio molecular-dynamics simulation of the liquid-metal-                      amorphous-semiconductor transition in germanium. Phys. Rev. B, 49:14251, 1994.

[3]     G. Kresse and J. Furthmüller. Efficiency of ab-initio total energy calculations for metals and            semiconductors using a plane-wave basis set. Comput. Mat. Sci., 6:15, 1996.

[4]     G. Kresse and J. Furthmüller. Efficient iterative schemes for ab initio total-energy                            calculations using a plane-wave basis set. Phys. Rev. B, 54:11169, 1996.

[5]     D. Vanderbilt. Soft self-consistent pseudopotentials in a generalized eigenvalue                              formalism. Phys. Rev. B, 41:7892, 1990.

[6]      G. Kresse and J. Hafner. Norm-conserving and ultrasoft pseudopotentials for first-row                   and transition-elements. J. Phys.: Condens. Matter, 6:8245, 1994.

[7]      J. P. Perdew, K. Burke, and M. Ernzerhof. Generalized gradient approximation made                      simple. Phys. Rev. Lett., 77:3865, 1996.

[8]     J. P. Perdew, K. Burke, and M. Ernzerhof. Erratum: Generalized gradient approximation                 made simple. Phys. Rev. Lett., 78:1396, 1997

[9]     Chen, M., Bates, S. P., Santen, van, R. A., & Friend, C. M. (1997). The chemical nature of                    atomic oxygen adsorbed on Rh(111) and Pt(111): a density functional study. Journal of                  Physical Chemistry B, 101(48), 10051-10057.

 

Determining the Lattice Parameters of Hf

Leave a reply

DFT and Lattice Parameters

For many metals (simple cubic, body-centered cubic, and face-centered cubic) using DFT to calculate the lattice parameter for a metal or crystal is rather straightforward as there is only one parameter to vary. One may begin by assuming the total energy of the system is a Taylor expansion of the lattice parameter, a, as below:

\begin{equation}E_{tot}(a) = E_{tot}(a_{0}) + \alpha (a-a_{0}) + \beta (a-a_{0})^2\end{equation}

Following Sholl (1), we may reduce this equation to:

\begin{equation}E_{tot}(a) = E_{tot}(a_{0}) + \beta (a-a_{0})^2\end{equation}

Thus, the lattice parameter for many metals and crystals can be determined by sampling various values of a. The energy of the system can then be calculated at each value of a and should look like a quadratic with a minimum at the true value of the lattice parameter.

Lattice Parameters for hcp Crystals

For hexagonal-close packed (hcp) structures, determining the lattice parameters is not so straightforward. For sc, fcc, and bcc metals there is only one parameter to optimize the energy. When considering an hcp metal, there are two lattice parameters on which the total energy depends: a, and c.

It is still possible, however, to manage this multi-variable minimization using our technique from above. First, we fix the ratio c/a to some value, r and sample various values of a at which to calculate the energy. With c/a fixed we will also know the value of c and we may construct our energy curves at varying values of r. The curve with the lowest total energy at its minima will be considered the “theoretical values” of the lattice parameters.

Lattice Parameters of Hf

Hafnium, element 72 on the periodic table, is a d-block transition metal with an hcp crystal structure. We wish to use DFT, as outlined above, to determine the equilibrium (ground state) lattice parameters of Hf.

Below are results obtained from CASTEP single point calculations for r = 1.40, 1.48, 1.58, 1.72, and 1.85 with an energy cut-off of 290 eV and a k-point grid of 8x8x6.

From this plot we can determine that the optimum value of r is about 1.58. Using this, a refined set of calculations may be performed at r = c/a=1.58 for various a to calculate an accurate estimate for both a and c. The results of these calculations are found below.

Fig 2. Refined determination of a.

In the above plot the blue points are CASTEP results and ther orange line is the harmonic approximation of the total system energy around the true lattice parameter. From this data we can observe a couple of things. First, we can say that the lattice parameters of Hf are approximately:

$$a_{0} = 3.22$$ $$c_{0} = 5.10$$

We may also see that the behavior of the energy is not truly quadratic with respect to the lattice parameter(s). The harmonic potential is a very good approximation for a near a0, but for a > a0 the harmonic potential is an overestimate of the energy and an underestimate for a < a0. This is because at larger separations (a > a0) the energy decreases as it should approach the dissociation energy of Hf (a → ∞) while at shorter separations (a < a0) the energy increases due to strong repulsive forces.

Fig. 3 Optimized Hf Cell

 

 

WebElements, an online reference for chemical elements, reports the lattice constants for Hafnium as (2):

$$a_{0}^{ref} = 3.20$$

$$c_{0}^{ref} =  5.05$$

Thus our calculations agree quite well with available data and we are satisfied.

Convergence

The above results were found using a relatively small energy cut-off (290 eV) and k-point grid (9x9x6) to allow for quick calculations that give an idea of the behavior or the energy with respect to the lattice constant. We now wish to see if we were converged with respect to the energy cut-off and k-point grid.

For the energy cut-off:

Fig 4. Energy Cut Off Convergence

While it appears that a rather high energy cut-off (~600 eV) is needed for convergence, it is important to realize that the energy differences, even between Ecut=250-270 are within chemical accuracy (~0.04 eV) and therefore using an Ecut of 290 should produce rather reliable results (as we have confirmed with literature values). Chemical accuracy is a standard used by computational chemists as a benchmark for making reliably accurate chemical predictions. Essentially, if energies are within ~0.04 eV then we can make confident predictions; this is exactly the case for our Ecut convergence.

Similarly for the k-points (where we have fixed the ratio kx/kz = ky/kz = 4/3) we find the following. Here our convention for “# of KPoints” is just to add the # of KPoints in each direction (e.g. 8x8x6 = 22 KPoints in the plot)

Fig. 5 K-Point Convergence

Again, we see that not very many k-points are needed (~20) before our energies are within chemical accuracy. This means that our use of an 8x8x6 k-point grid is sufficient for our needs.

Single Point Calculations

For the single point calculations performed for Hf with varying lattice constants the following (ultra-fine) calculation settings were used:

Exchange Correlation Functional Type: Generalized Gradient Approximation (GGA)

Exchange Correlation Functional: PBE

Plane-Wave Energy Cut-Off: 290.0 eV

K-Point Grid: 9x9x6

Pseudopotentials: Ultrasoft

Electronic Energy Convergence Criteria: 5.0 10-7 eV/atom

Geometry Optimization

For the geometry optimization of Hf, the following (ultra-fine) calculation settings were used:

Exchange Correlation Functional Type: Generalized Gradient Approximation (GGA)

Exchange Correlation Functional: PBE

Plane-Wave Energy Cut-Off: 290.0 eV

K-Point Grid: 9x9x6

Pseudopotentials: Ultrasoft

Electronic Energy Convergence Criteria: 5.0 10-7 eV/atom

Ionic Energy Convergence Criteria: 5.0 10-6 eV/atom

Ionic Force Convergence Criteria: 0.01 eV/Å

Future Work:

In the future it would be of interest to investigate the effect of the type of pseudo-potential used (e.g. ultrasoft, soft, etc.) as well as the psuedo-potential cut-off for treating the core and valence electrons. Also of interest would be how our results change using different XC Functionals.

References:

(1) Sholl, D. & Steckel, J. A. Density Functional Theory: A Practical Introduction. (John Wiley & Sons, 2011).

(2) https://www.webelements.com/hafnium/crystal_structure.html

(3) First principles methods using CASTEP”, Zeitschrift fuer Kristallographie 220(5-6) pp.                  567-570 (2005)  S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K.                    Refson, M. C. Payne