Theory

In this post we report our findings on our study to determine whether pure Li or LiH is more energetically favorable for different hydrogen pressures P\(_{H_2}\) and temperatures T. For this we need to determine which material minimizes the free energy of the a system containing gaseous H\(_2\) and Li at a given pressure P\(_{H_2}\) and temperature T. We use the grand potential of a given system as a proxy for the free energy since the grand potential of a system is lowest for a system with lowest free energy.

The grand potential for a system containing N\(_{Li}\) Li atoms and N\(_H\) hydrogen atoms is defined by,

\begin{equation}\Omega\left(T,\mu_{Li},\mu_H,N_{Li},N_H\right)=E\left(N_{Li},N_H\right)-TS-\mu_H N_H- \mu_{Li}N_{Li}.\label{gp_def}\end{equation}

Where \(E\left(N_{Li},N_H\right)\) is the internal energy of the material and is determined using a DFT calculation, \(S\) is the material entropy and respectively \(\mu_{species}\) and \(N_{species}\) define the chemical potential and number of atoms of the given species of atoms in the material.

Now we are interested in comparing the grand potentials for Li and LiH. Before proceeding with our calculations we make a few adjustments to equation \eqref{gp_def} to simplify our calculations. First the difference in internal energy between Li and LiH is much larger than the entropic differences between, hence the second term in the right hand side of equation \eqref{gp_def} can be neglected. Next we can eliminate the last term in equation \eqref{gp_def} by choosing to look at  a pure Li unit cell and  a LiH unit cell with equal number of Li atoms. Since \(\mu_{Li}\) is independent of P\(_{H_2}\) and we set \(\left(N_{Li}\right)_{pure\;Li}= \left(N_{Li}\right)_{LiH}\), the term \(\mu_{Li}N_{Li}\) is the same for both systems and will cancel out when considering the difference in grand potential between the two systems.

With above mentioned simplifications equation \eqref{gp_def} now reads as,

\begin{equation}\Omega\left(T,\mu_H,N_H\right)=E\left(N_H\right)-\mu_H N_H.\label{gp_red}\end{equation}

In equilibrium condition the chemical potential of molecular hydrogen and atomic hydrogen are related by,

\begin{equation}\mu_{H}=\frac{1}{2}\mu_{H_2}.\label{muH2_H}\end{equation}.

Treating H\(_2\) as an ideal gas the definition of it’s chemical potential is,

\begin{equation}\mu_{H_2}=E_{H_2}^{total}+\tilde{\mu}_{H_2}\left(T,p^0\right)+k_BTln\left(\frac{p}{p^0}\right).\label{H2_ideal_chem_pot}\end{equation}

Here \(k_B\) is the Boltzmann constant, \(E_{H_2}^{total}\) is the total energy of an isolated H\(_2\) molecule and can be obtained from a DFT calculation for an isolated H\(_2\) molecule. The second term in equation \eqref{H2_ideal_chem_pot}, \(\tilde{\mu}_{H_2}\left(T,p^0\right)\) is the difference between the \(T=0\;K\) and the temperature of interest chemical potential of molecular H\(_2\). \(\tilde{\mu}_{H_2}\left(T,p^0\right)\) can be calculated using data tabulated in NIST-JANAF Thermochemical Tables for molecular hydrogen. In terms of quantities defined in the tables referred to; by noting the chemical potential of an ideal gas is equal to the free energy \(G=H-TS\) of an ideal gas,

\begin{equation}\tilde{\mu}_{H_2}\left(T,p^0\right)=\left[H^0\left(T\right)-H^0\left(T_r\right)\right]-TS\left(T\right)-\left[H^0\left(0\right)-H^0\left(T_r\right)\right].\label{mu_tilde}\end{equation}

Calculations

Now we have the recipe to obtain the Li-LiH phase diagram. From DFT we need the internal energies of a pure Li unit cell, LiH unit cell and a Hydrogen molecule. All energy calculations were performed using DFT calculations using the DMol3 software package as implemented by Materials Studio. DMol3 employs numerical radial functions centered at atomic positions as basis sets. The GGA base PBE functional was used to treat exchange correlational effects and the double numerical precision – DNP 3.5 basis set was used.  All energy calculations were deemed to be converged once an SCF convergence of \(10^{-6}\) was achieved. Additionally for the periodic systems a k-point mesh with \(0.03\;\mathring{A}^{-1}\) maximum spacing between mesh points to have the internal energies calculated to an accuracy of \(1\;m\,eV\).

The systems for which we require the internal energies to perform our calculations are a pure Li primitive unit cell, a LiH primitive unit cell and the hydrogen molecule. First the primitive unit cells of Li and LiH and also the hydrogen molecule were geometry optimized to obtain the optimal structures. The optimal structures of the Li unit cell (primitive unit cell of a BCC lattice with lattice constant = 3.0155 Å), LiH unit cell (primitive unit cell of a FCC lattice with lattice constant = 2.9409 Å) and the hydrogen molecule (H-H bond length = 0.748444 Å)  are shown in figures 1, 2 and 3 respectively. The internal energies of the systems of interest are shown in table 01.

Figure 1 : Primitive unit cell of Li. This is a primitive unit cell of a BCC lattice with lattice constant = 3.0155 Å.

Figure 2 : Primitive unit cell of LiH. This is a primitive unit cell of a FCC lattice with lattice constant = 2.9409 Å.

Figure 3 : Hydrogen molecule. Hydrogen bond length = 0.748444 Å.

Table 01 : Internal Energies of systems of interest.

Now using the primitive unit cells of Li and LiH and comparing their respective grand potentials, we see that \(\Omega_{Li}=\Omega_{LiH}\) or Li and LiH are thermodynamic equilibrium when,

\begin{equation}\mu_H^{eq}=E_{LiH}-E_{Li}=-16.573\;eV;\label{req_mu_H}\end{equation}

Or when, (using equation \eqref{muH2_H})

\begin{equation}\mu_{H_2}^{eq}=-33.146\;eV.\label{req_mu_H2}\end{equation}

Now in order to determine the phase diagram we need to determine the pressure at which \(\mu_{H_2}=-33.146\;eV\) for a given temperature. To do this we look at equation \eqref{H2_ideal_chem_pot}. From the NIST-JANAF Thermochemical Tables for molecular hydrogen we have obtained the plot , shown in figure 4, of \(\tilde{\mu}_{H_2}\left(T,p^0\right)\) vs \(T\) for a reference pressure \(p^0=0.1\;MPa\).

Figure 4 : Plot of \(\tilde{\mu}_{H_2}\left(T,p^0\right)\) vs \(T\) for a reference pressure \(p^0=0.1\;MPa\).

Now for a given temperature we can solve for the pressure \(p\) that makes \(\mu_{H_2}=\mu_{H_2}^{eq}\), by rewriting equation \eqref{H2_ideal_chem_pot} as,

\begin{equation}\ln\left(\frac{p^{eq}}{p^0}\right)=\frac{\mu_{H_2}^{eq}-E_{H_2}^{total}-\tilde{\mu}_{H_2}\left(T,p^0\right)}{k_BT}.\label{solve_for_p}\end{equation}

The curve \(\ln\left(\frac{p^{eq}}{p^0}\right)\) vs \(T\) which is the phase boundary between Li and LiH is shown in figure 5. For any point in phase space, corresponding to a point above and/or to the left of the \(\ln\left(\frac{p^{eq}}{p^0}\right)\) vs \(T\) curve, LiH is the thermodynamically preferred structure. Conversely, any point in phase space, corresponding to a point below and/or to the right of the \(\ln\left(\frac{p^{eq}}{p^0}\right)\) vs \(T\) curve, pure Li is the thermodynamically preferred structure.

Figure 5 : Plot of \(\ln\left(\frac{p^{eq}}{p^0}\right)\) vs \(T\).

Introduction

This post will focus on applying transition state theory to hydrogen atoms adsorbed onto the surface of Cu(111). We first use geometry optimization to find the preferred position of hydrogen on the FCC and HCP sites on the Cu(111) surface. After this, transition state calculations using the complete linear synchronous transit and quadratic synchronous transit (LST/QST) method are performed in order to find the energy barriers [1].

The Slab Model and Geometry Optimization

We use a slab model consisting of 4 layers of Cu and a single hydrogen atom adsorbed to the surface. The vacuum is chosen to be 10 Å thick, and the lattice vectors are (2.21 Å, -1.28 Å, 0), (0, 2.56 Å, 0), and (0, 0, 16.26 Å).

Fig 1: The slab model for Cu(111) with a hydrogen adsorbed onto the surface.

All calculations were performed with the CASTEP package in materials studio with the GGA Perdew-Burke-Ernzerhof functional [2-3]. The pseudopotential was solved with the Koelling-Harmon atomic and 1s2 2s2 2p6 3s2 3p6 for the inner shells, while 3d10 4s1 were used for the outer shells.

Geometry optimization calculations used convergence tolerances of 10^-5 eV/atom, 0.03 eV/Å for the max force, 0.05 GPa for the max stress, and a max displacement of 0.001 Å. The k-point set is increase from 6x6x1 to 10x10x1, and then from 10x10x2 to 12x12x2 to test for convergence. In changing to k-point set from 12x12x2 to 13x13x2, there is less than a 0.01 eV change in the system energy. Convergence for the cutoff energy is tested by starting at 440 eV and increasing until 570 eV. This test shows that there is less than a 0.01 eV change in the energy when going from 480 eV to 490 eV.

Fig 2: Test of the convergence for the cutoff energy. An identical trend is observed for hydrogen at the HCP site.

Transition State Calculations

The trajectory for the hydrogen transition was predicted using the Materials Studio “reaction preview” tool.

Fig 3: Side view of hydrogen moving from the HCP site to the FCC site on the surface of Cu(111) as determined by the reaction preview tool.

Fig 4: Top-down view of hydrogen moving from the HCP site to the FCC site on the surface of Cu(111) as determined by the reaction preview tool.

The transition state and energy barrier are then calculated using the LST/QST method. The RMS convergence is chosen to be 0.05 eV/Å with the maximum number of QST steps set to 5. Calculations were done to show that the energy barriers are converged well for changes of 20 eV to the energy cutoff for cutoffs past 480 eV. The reaction pathway and the energies are shown in fig 4. The transition state is found to have an energy of -6736.2477 eV.

Table 1: Variation of the energy barriers with increasing energy cutoff.

Fig 4: Reaction pathway for a hydrogen jumping from the HCP site to the FCC site. The blue star indicates the position of the transition state.

We have previously found that the hydrogen atom at the HCP site has vibrational frequencies of 36.9 THz, 32.7 THz, and 33.4 THz, while at the FCC site it has vibrational frequencies of 32.1 THz,  31.8 THz, and  31.1 THz. We can once again apply the Hessian matrix method to find the vibrational frequencies for the transition state. These frequencies give zero-point energy corrections of 0.213 eV and 0.196 eV, respectively.

Table 2: Perturbations in the energy for small changes in the hydrogen atom position at the transition state. These are used to calculate the Hessian matrix.

Doing this, we find that the transition state has real frequencies of 63.2 Thz and 32.0 THz. There is one imaginary frequency of 23.0 Thz, which corresponds to the hydrogen moving from the HCP site to the FCC site or vice-versa. The corresponding zero-point energy correction is then 0.197 eV.

Using the zero-point energy corrections, we can then find the corrections to the energy barriers. This correction is defined as

\[\Delta E_{ZP} = \bigg( E^{\dagger} + \sum_{j=1}^{N-1} \frac{h \nu_{j}^{\dagger}}{2} \bigg) – \bigg( E_{A} + \sum_{i=1}^{N} \frac{h \nu_{i}^{\dagger}}{2} \bigg) \]

Here, \(E^{\dagger} \) is the energy of the transition state and \(E_{A} \) is the energy in either the HCP site or the FCC site. The sum with j sums over the N-1 real modes at the transition state, while the sum over i sums over the normal modes in either the FCC or HCP site. These two sums simply give the zero point energy corrections that we have calculated above.

For a transition with hydrogen starting in the HCP site, the energy barrier will be corrected to 0.230 eV. For a transition from the FCC site, the correction to the energy barrier is 0.219 eV.

Conclusions

Using transition state theory, we find that hydrogen hopping from the HCP to FCC vacancy site on the surface of Cu(111) has an energy barrier of 0.247 eV, while hopping in the other direction has an energy barrier of 0.220 eV. Using the zero-point energy correction, we find that these barriers are corrected to 0.230 eV and 0.219 eV, respectively. These correspond to 4.3% and 0.5% changes from the original energy barriers.

References

[1] Govind, N., Petersen, M., Fitzgerald, G., King-Smith, D., & Andzelm, J. (2003). A generalized synchronous transit method for transition state location. Computational Materials Science,28(2), 250-258.

[2] First principle methods using CASTEP Zeitschrift fuer Kristallographie 220(5-6) pp. 567-570 (2005).

[3] Perdew, J.P., K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple.Physical review letters, 1996. 77(18): p. 3865.

Introduction

In this project the phase diagram for LiH and Li will constructed.  To do this three single point energy calculations must be preformed, one for the LiH crystal, the Li crystal and molecule. Before the single point energies could be calculated, the geometry of the crystals and molecule were optimized.

The calculations were performed with CASTEP [2] and used the GGA XC-functional PBE [1], OTFG ultrasoft pseudopotentials, the Koelling-Hamon relativistic treatment, a k-point grid of 11x11x11 for the  crystals and 2x2x2 for for the and a energy cutoff of 600eV. The SCF tolerance was set to eV/atom with a convergence window of three steps, For the geometry optimizations the energy convergence was set to eV/atom, the max force was set to 0.01 eV/A, the max stress was set to 0.02 GPa, and the max displacement was set to 0.0005A. For the two crystals a full cell optimization was used with hard comprehensibility.

The k-point grid used converged the energy of the Li and LiH crystals to within 0.005eV and energy cutoff used converged the energy of the Li and LiH crystals to within 0.02eV.

Pseudo atomic configuration Li: 1s2 2s1

Pseudo atomic configuration H: 1s1

Optimization

The geometries were optimized using the settings above and the results are shown in Figure 1, 2 and 3

Figure 1: Primitive bcc lattice of Li, a = 2.983Å

Figure 2: Primitive NaCl type lattice of LiH, a = 2.835Å

Figure 3: hydrogen molecule with bond length 0.752Å in a 15Å sided cube of vacuum.

The Grand Partition Function

To find the equilibrium point between LiH and Li the grand partition functions, [3] of the two were set to be equal and the chemical potential of atomic hydrogen, , was solved for

     (eq 1)

     (eq 2)

where and are the internal energies of Li and LiH crystals respectively,  is the number of hydrogen atoms in the LiH crystal and is an additive constant that is the same for all materials.  When those two equations are combined the equilibrium chemical potential of atomic hydrogen can be expressed as,

     (eq 3)

In addition the, chemical potential of atomic hydrogen can be expressed in terms of the chemical potential of molecular hydrogen ,

     (eq 4)

Assuming that the hydrogen behaves like an ideal gas one can write the chemical potential of molecular hydrogen as,

     (eq 5)

where is the total energy of an isolated molecule at T=0K, is the difference in the chemical potential of between T=0K and the temperature of interest at the reference pressure, is the reference pressure, p is the pressure of the system, and T is the temperature of the system.

Chemical Potential Difference

The values for the difference in the chemical potential of  between T=0K and the temperature of interest at the reference pressure can be evaluated using data tabulated in the NIST-JANAF Thermochemical Tables and the following equation

     (eq 6)

Where ;  , the difference in the enthalpy at temperature T and the enthalpy at the reference temperature;  and S(T), the entropy at temperature T, are values given in the reference table. Selected values are listed in Table 1, along with the calculated value of .

Table 1: Values from NIST-JANAF Thermochemical Tables and the calculated difference in chemical potential.

Results

Single point calculations were performed on the optimized geometries for the Li crystal, the LiH crystal and the molecular hydrogen, the results of which are listed in Table 2.  Using equations 3 and 4 the value of the chemical potential of atomic hydrogen and molecular hydrogen were found, the resulting values can be found in Table 2.

Table 2: Results of single point energy calculations for the three systems.

As we know the value of , and equation 5 can be rewritten to solve for p.

     (eq 7)

Since the relationship of temperature to pressure is known at the equilibrium point between LiH and Li, the phase diagram can be constructed and is presented in Figure 4.

Figure 4: Phase diagram of Li and LiH

References

1. Perdew, J. P.; Burke, K.; Ernzerhof, M. Physical Review Letters 1996, 77 (18), 3865–3868.

2. Clark, S. J.; Segall, M. D.; Pickard, C. J.; Hasnip, P. J.; Probert, M. I. J.; Refson, K.; Payne, M. C. Zeitschrift für Kristallographie – Crystalline Materials 2005, 220 (5/6).

3. Sholl, D.; Steckel, J. A.; Sholl. Density Functional Theory: a Practical Introduction; Wiley: Somerset, 2011.

This project aims to apply Kinetic Monte Carlo to study dynamic H2 dissociation on a homogeneous surface and to approach traditional Langmuir adsorption model. Fichthorn et al. presented details of theoretical basis for dynamic Monte Carlo and gave an example of non-dissociative adsorption [1]. Based on their research, we utilized Kinetic Monte Carlo to further explore dissociative adsorption on surface.

Kinetic Monte Carlo Parameters Set Up
We assume each H atom adsorbed on surface does not interact with another H atom nearby and occupy one adsorption site. And there is no site difference in our model. Adsorption equilibrium will be achieved when adsorption rate is equal to desorption rate.There are basically two criterion we should choose carefully. The first one is to decide event happening probability. In H2 dissociatively adsorption, there are two events in total, namely, H2 dissociatively adsorb on surface, with rate of rads and 2H associatively desorb from the surface to form H2, with rate of rdes. Therefore, we choose adsorption probability as Pads=rads/(rads + rdes). If randomly selected number i is larger than Pads, then we select two neighboring H atoms to desorb from the surface and if i is smaller than Pads, then we select two neighboring vacancy sites on surface to adsorb 2H atoms. The second criterion is the time increment, which can be described as 1/(rads + rdes), based on first criterion. rads and rdes can be calculated as the following equations:

where theta is the surface coverage, kads is the adsorption rate constant and kdes is the desorption rate constant, PH2 is H2 pressure, which is set as 1 atm in our case. Since kads and kdes vary among different surfaces, we randomly choose kads=3 and kdes=2 in our case. But specific value can be obtained through calculating activation barrier and transition state theory depending on different surfaces. In the following section, we will discuss H2 dissociation on Pd(111) surface.

Results
Surface lattice and coverage are shown in Figure 1 and Figure 2.

Figure 1. The surface lattice at equilibrium state with blue dots representing H atoms

Figure 2. Kinetic Monte Carlo derived surface coverage of H and Langmuir adsorption analytical solution

H2 dissociation on Pd(111) surface
As mentioned above, transition state theory can be applied to obtain H2 dissociation rate constant. We take 5 layer 3X3 Pd (111) surface slab as an example. The Vienna ab initio simulation program (VASP) was used for all DFT calculations [2-4]. The Perdew−Burke−Ernzerhof (PBE) generalized gradient approximation [5] was applied with a plane wave basis set with energy cutoff of 450 eV. A 4 × 4 × 1 Monkhorst−Pack k-point mesh for surface slab was used [6]. Structural optimization was carried out until the force on all atoms was less than 0.05 eV/Å. The climbing image nudged elastic band method (CI-NEB) was used to find transition state [7,8]. Optimized initial state, transition state and final state are shown in Figure 3. It is notable that the key point of this part is to illustrate how to obtain rate constant based on DFT transition state search. Exploration about the most stable final structure is not involved in this work. Usually, hollow site is the most favorable for H adsorption. But for the sake of time and illustrating the key point, we simply chose atop site as the final state. Activation barrier obtained from DFT calculation is 0.015 eV and corresponding rate constant under 298 K is 1.1X10^13/s.

Figure 3. H2 dissociation on Pd (111) initial state, transition state and final state

In conclusion, we utilized Kinetic Monte Carlo to explore dynamic change of surface coverage during H2 dissociation and equilibrium coverage is in consistency with Langmuir adsorption model. We further illustrate an example about how to obtain H2 dissociation rate constant based on transition state theory.

Reference
[1] K. Fichthorn and W. Weinberg. J. Chem. Phys 95 1090 (1991)
[2] G. Kresse and J. Furthmüller. Comput. Mater. Sci. 6 (1996)
[3] G. Kresse and J. Furthmüller. Condens. Matter Mater. Phys. 54 (1996)
[4] G. Kresse and J. Hafner. Phys. Rev. B: Condens. Matter Mater. Phys. 47 (1993)
[5] J.P. Perdew; K. Burke and M. Ernzerhof. Phys. Rev. Lett. 77 (1996)
[6] H.J. Monkhorst and J. D. Pack. Phys. Rev. B. 13 (1976)
[7] G. Henkelman; B. P. Uberuaga and H. Jonsson. J. Chem. Phys. 113 (2000)
[8] G. Henkelman; and H. Jonsson. J. Chem. Phys. 113 (2000)

Abstract

In this paper, we report a systematic method of calculating phase diagram for metal hydride compound using DFT. While our calculation is limited to perfect crystal phase, it is not hard to generalize this method to more complicated systems that could even contain defects, surface phases, etc.

Structure of Li, LiH and Hydrogen Molecule

The lithium is body-centered cubic (bcc) with lattice constant 3.509Å[4] and its hydride LiH is NaCl type face-centered cubic (fcc) with lattice constant 4.007Å [1] . DFT calculations were performed to get the optimized lattice constant for LiH and the total internal energy of the lattice for both crystals, shown in figure 1.

DFT calculation is also performed to obtain the total internal energy of hydrogen molecule. This is done by fix the hydrogen molecule at a very big cubic lattice(10Å by 10Å by 10Å) and perform CASTEP calculation (figure 2). Details of calculation can be found in appendix.

Figure 1 Structure of Lithium and Lithium hydride. Li is bcc lattice while LiH is NaCl structure.

Figure 2 Hydrogen molecules. The energy is calculated by performing DFT in simple cubic lattice. The inter-molecule distance is far more greater than the size of the molecule, and the bond distance between hydrogen atoms is 0.7414Å[1].

Grand Potential for LiH and Li

In general, the grand potential of metal hydride obains the following form[2]:

\begin{equation}\Omega=E_i-T S_i-\mu_{O,i}N_{O,i}-\Omega^M\end{equation}

The first term is the energy of the metal or metal-hydride. The second term represents the effect of entropy of crystals. The lowest order approximation is to take the entropy of the crystals zero. The third term is the chemical potential of hydrogen, which depends on the number of hydrogen atoms in the crystal. The last term is a constant for two phases in our report, which is the same since both are lithium atoms. So it could be taken as zero. To be explicit, the grand potential for LiH and Li are

\begin{equation}\Omega_{LiH}=E_{LiH}-\mu_{H,i}N_{H, LiH}\end{equation}

\begin{equation}\Omega_{Li}=E_{Li}\end{equation}

In this report, we compute the grand potential per metal atom. Hence . The internal energy given by DFT are

\begin{equation}E_{Li}=-202.44eV\end{equation}

\begin{equation}E_{LiH}=-219.13eV\end{equation}

At equilibrium, the grand potential for LiH and Li are equal. By setting two equations equal to each other, one obtains the following chemical potential relationship:

\begin{equation}\mu_{H}=\frac{1}{2}\mu_{H_2}=\frac{E_{LiH}-E_{Li}}{N_{H,LiH}}=-16.69eV\end{equation}

The chemical potential of hydrogen molecule could be evaluated by standard thermodynamic method:

\begin{equation}\mu_{H_2}=E_{H_2}^{total}+\tilde{\mu}_{H_2}(T, p^o)+k_B T \ln{(p/p^o)}\end{equation}

Where is the total energy of hydrogen molecule and is the standard chemical potential at standard pressure and temperature, which could be obtained by calculating Gibbs free energy[2]:

\begin{equation}\tilde{\mu}_{H_2}(T, p^o)=[H^o(T)-H^o(T_r)]-TS(T)-[H^o(0)-H^o(T_r)]\end{equation}

The standard value could be found in the NIST-JANAF Thermochemical Tables[3]. 

Phase Diagram

The phase diagram could be obtained by solving for :

\begin{equation}\ln{p/p^o}=\frac{1}{k_B T} (\mu_{H_2}-E_{H_2}^{total}-\tilde{\mu}_{H_2}(T, p^o))\end{equation}

The results are shown in Table 1. It is noticeable that at low temperature and high pressure, the system tends to stay at LiH while at high temperature with low pressure, the system is more stable in Li phase. At pressure close to atmospheric pressure, the transition temperature is around 1200K, this result agrees with the experimental experience that LiH can be synthesize by Li and hydrogen at several hundreds of degrees.

Figure 3 Phase diagram of LiH/Li.

Appendix

DFT Calculation Details of Li

xc_functional : GGA PBE
spin_polarized : false
cut_off_energy : 600.0eV
calculate_stress : false
Atomic calculation performed for Li: 1s2 2s1
Pseudo atomic calculation performed for Li 1s2 2s1
number of bands : 6
k-Points For BZ Sampling: 10*10*10

DFT Calculation Details of LiH

xc_functional : GGA PBE
page_wvfns : 0
cut_off_energy : 600.0eV
Atomic calculation performed for H: 1s1
Atomic calculation performed for Li: 1s2 2s1
Pseudo atomic calculation performed for H： 1s1
Pseudo atomic calculation performed for Li： 1s2 2s1
k-Points For BZ Sampling: 6*6*6

DFT Calculation Details of Hydrogen molecule

Key parameters:
xc_functional : GGA PBE
spin_polarized : false
cut_off_energy : 450.0eV
Atomic calculation performed for H: 1s1
k-Points For BZ Sampling: 1*1*1
Lattice constants: 10.0A*10.0A*10.0A

Reference

[1] https://en.wikipedia.org/wiki/Hydrogen

[2] Sholl, David, and Janice A. Steckel. Density functional theory: a practical introduction. John Wiley & Sons, 2011.

[3] Stull, Daniel Richard, and Harold Prophet. JANAF thermochemical tables. No. NSRDS-NBS-37. National Standard Reference Data System, 1971.

[4] Dassault Systèmes BIOVIA, BIO, BIOVIA Materials Studio 2018 (18.1.0.2017), San Diego: Dassault Systèmes, 2017.

1. Project Description

In this project, we aim to find the most stable configuration of the defective hematite (\(Fe_2\)\(O_3\)) by comparing the energy of each configuration. We modeled a 3*3*1 hematite supercell first, took one of the iron atoms out of the structure, and added three hydrogens with stoichiometry. The iron occupancy of the structure was therefore 0.991. After modeling the initial configuration, Density Functional Theory (DFT), with the Vienna Ab-initio Simulation Package (VASP) code[1], was used to optimize lattice parameter and calculate the energy. The exchange-correlation was treated in the generalized gradient approximation (GGA), as parameterized by Perdew, Burke, and Ernzerhof (PBE) [2]. GGA+U was applied to correct self-interaction of the d-electrons in Fe atoms. The following valence electron configurations were employed in the potentials: 3d6 3d7 4s1 for Fe, 2s2 2p4 for O, and 1s1 for H. Psedupotentials were used in place of other electrons, with core radii of 2.200Å for Fe, 1.100Å  for O, and 0.800Å for H. The Brillouin zone integration was performed using a 1*1*1 Monkhorst-Pack k-point grid, and plane-wave kinetic energy cutoff was 800eV [3]. The SCF tolerance is 1.0e-4 eV/atom. The spin states of iron atoms in the adjacent layers were treated with opposite spin states.

2. Crystal Structure and Initial Defective Configurations

The unit cell of \(Fe_2\)\(O_3\) used Blake et al. (1970) bulk hematite structure that a=b=5.038Å, c=13.772Å, α=β=90 degrees, γ=120 degrees, and the space group is R-3c [4]. The equivalent fractional coordination in bulk hematite for iron is (0, 0, 0.35), for oxygen is (0.31, 0, 0.25). Iron occupies two-thirds of slightly distorted octahedron sites, and oxygen locates at the tetrahedron site. The Fe-O3-Fe unit repeated along the z-axis and stacked oxygens form a hexagonal close-packed (HCP) structure (Figure 1).

Figure 1. Crystal structure of hematite unit cell (\(Fe_2\)\(O_3\), Left) from Blake et al. (1970) and hematite 3*3*1 supercell (\(Fe_{108}\)\(O_{162}\), Right) after geometry optimization. View from the [110] direction. Red is oxygen; blue is iron.

Figure 2. FeO6 octahedron (left), and O6H3 structure (right) by removing one iron and adding three hydrogens

We began with modeling one iron defect (iron occupancy 0.991) in the hematite supercell (\(Fe_{108}\)\(O_{162}\)). To compensate one defective iron (Fe3+) and reach a stoichiometry in the structure, three hydrogens were modeled and bonded with three oxygens individually (Figure 2). Because the octahedron site surrounded by six oxygens is slightly distorted in hematite structure (Figure 3), the hydrogens attached to different oxygens will generate different energy. Our goal is to find the preferred oxygen position for the hydrogen to bond with by comparing the energies of 8 different configurations.

.

Figure 3. The geometry of the distorted octahedron site in hematite structure. Blue and pink atoms are iron but spin up and down respectively. Red atom is oxygen.

3. Result and Discussion

After geometry optimization and energy convergence, the energy of the 3*3*1 hematite supercell is -2002.174eV. The lattice parameters of the supercell showed that a=b=15.231Å, c=13.855Å, α=β=90 degrees, and γ=120 degrees (Figure 1). Compared to the experimental hematite structure, the computational error is 0.772%.

Figure 4. Energy convergence of configuration O 1-2-3.

The energies of eight configurations are summarized in table 1, and one of the energy convergence plots is shown in Figure 4. The lowest energy of different configurations is -2156.654eV. The most stable structure is therefore formed where hydrogens bonded with O3, O4, and O5 (Figure 5). Accordingly, the lattice parameter of this configuration is a=15.090Å, b=15.088Å, c=13.806Å, α=β=90 degrees, and γ=120 degrees. Compared to the experimental lattice parameter of hematite, the error of the computational lattice parameter is within 0.2%.

Table 1. The energy of each configuration showed that hydrogens preferred to bond with O3, O4, and O5.

Figure 5. The most stable structure for hydrogens to fit in the octahedron before (left) and after (right) geometry optimization.

The energy differences among various configurations can be caused not only by the relative atom positions but also the hydrogen bonds. The predicted OH bond length in our DFT calculation ranges from 0.93-1.0Å, which agrees with Pinney et al. (2009) [3]. For our most stable configuration O3-4-5, the OH bond length is 0.982Å (O3H), 0.991Å (O4H), and 0.992Å (O5H). Two weak hydrogen bonds are O3H—-O2 and O3H—-O1, with bond length 2.2000Å and 2.366Å respectively. One stronger hydrogen bond formed between O5H and O1, with bond length 1.824Å.

Reference

[1] Kresse and J. Furthmuller, Phys. Rev. B, vol. 54, 1996, pp. 11169-11186.

[2] P. Perdew, B. Kieron and E. Matthias, “Generalized gradient approximation made simple,” Physical review letters, vol. 77, no. 18, p. 3865, 1996.

[3] Pinney, N., Kubicki, J. D., Middlemiss, D. S., Grey, C. P., & Morgan, D. (2009). Density functional theory study of ferrihydrite and related Fe-oxyhydroxides. Chemistry of Materials, 21(24), 5727-5742.

[4] L. Blake, R. E. Hessevick, T. Zoltai and L. W. Finger, “Refinement of hematite structure,” American Mineralogist, vol. 51, p. 123, 1966.

Bi2Se3, known as a topological insulator, has attracted much attention recently. In this material, due to the strong spin-orbit coupling (SOC), the conduction band and the valence band of Bi2Se3 are inverted in momentum space, which shows a Dirac-corn like surface state. For this project, I want to calculate the bulk band structure of Bi2Se3 with and without SOC.

1.The crystal structure of Bi2Se3

The crystal structure of Bi2Se3 is shown in Figure 1. The yellow atoms are Se atoms and the purple atoms are Bi atoms. Bi2Se3 has a rhombohedral structure of the space group D5 3d with R3m (#166). The lattice parameter are a=b=4.113 Å, c=28.5859 Å and the bond angle are a=b=90 degrees and c=120 degrees.

Figure 1, the crystal structure of Bi2Se3

2. Calculation of band structure without SOC

After the construction of the crystal, the geometry optimization is performed and the band structure is calculated, which is shown in Figure 3 and 4. As shown in these two figures, without considering SOC, the Bi2Se3 has a direct band gap, which is about 0.505 eV.

Figure 2. The structure of Bi2Se3 with the Brillouin zone (BZ). The blue lines show the Brillouin zone (the primitive unit cell of the reciprocal lattice) with the high symmetry points G (Center of the BZ), A (Center of a hexagonal face), H (Corner point), K (Middle of an edge joining two rectangular faces), M (Center of a rectangular face), L (Middle of an edge joining a hexagonal and a rectangular face).

Figure 3 The band structure of Bi2Se3 without SOC

Figure 4 The band structure (near G point) of Bi2Se3.

3. Calculation of band structure with SOC

The band structure with SOC is shown in Figure 5 and Figure 6. By including SOC, Bi2Se3 has an indirect band gap, which is about 0.281 eV. By comparing figures in Section 2 and Section 3, it can be seen that strong SOC has a significant impact on the band structure, and changes the band gap at the G point. Also since there is no direct band gap, it is possible that SOC induces the inversion between conduction and valence band at the G point, which could be a signal of a non-trivial topology.

Figure 5 The band structure of Bi2Se3 with SOC

Figure 6 The band structure of Bi2Se3 (near G point) with SOC

4.Calculation details

The pseudo atomic calculation performed for Bi: 6s2 6p3

The pseudo atomic calculation performed for Se: 3d10 4s2 4p4

Functional: GGA PBE

Cut off energy: 500eV

k-point set: 3*3*1

SCF tolerance: 1.0*10^ (-5) eV

Pseudopotential: Norm conserving

5. Conclusion and Discussion

In this project, the band structure of Bi2Se3 with and without SOC is calculated and compared. From the band structure, it can be seen that without changing other parameters, SOC can significantly change the band structure. So it is important to include SOC, especially for the elements with the atomic number being greater than the value of 50. Also by considering SOC, there is a signal of band inversion at G point in Bi2Se3, which could be an evidence of a non-trivial topology.

In this project, you can see that the parameters I choose are not very accurate. For one reason, I just want to point out the importance of spin-orbital coupling, so the accuracy of the results is not very high. For another reason, since the structure of Bi2Se3 are very complicated which needs huge calculation time and physical memory to calculate, it is impossible to improve the accuracy just using my computer.

For some future works, using the material studio, we can also get some other properties of Bulk Bi2Se3, such as the optical properties. Also, it is possible that we can build a slab with just a few layers to study the properties of the Bi2Se3 surface.

Purpose of the Calculation

One of the quantities of interest when doing DFT calculations is the electronic density of states (DOS). The DOS is defined as the number of electron states with an energy in the interval E + dE, and is one method for determining the electronic state of the material without examining the full band structure. Our test case for studying the convergence of the DOS will be silver, a metal. The indicative feature of a metal in the DOS is that there are a continuum of available states from below the Fermi energy to above the Fermi energy, which is defined as the energy of the highest energy occupied electron state at zero temperature. This post addresses the topic of properly converging the DOS calculation with respect to k-point mesh, and how that convergence varies from the convergence of the system energy with respect to k-point mesh.

Fig: 1 – The silver unit cell used in this convergence study

Theory of Our Calculation

As with any result that we might report, if we choose to compute the DOS, we must ensure that the result is appropriately converged with respect to the input parameters. To calculate the DOS, we integrate the converged electron density of the system in k-space [1]. Naively, we might then expect that the convergence of the electronic DOS would be similar to the system energy. However, we need to consider the fact that the DOS is a continuous function of energy that needs to be well-resolved for the entire range which we are covering.

The only points at which we calculate the electron density are the points on our k-point mesh. Computationally, we evaluate the integral by doing a weighted sum over the points that we have, using Simpson’s Rule or a similar algorithm. However, when doing something like a Simpson’s Rule algorithm, the function is interpolated in between the points with known values. With a function that rapidly varies in narrow regions, Simpson’s Rule will not give at accurate depiction of the true integral, unless the region is well-sampled. As a result, we should expect that the DOS will not be well-converged without many more k-points than is necessary for the convergence of the energy of the system. We should also expect that this convergence should take the form of the density of states becoming more smoothly varying in the regions where it is non-zero, as the under-sampled regions become appropriately represented.

Calculation Methodology

We use CASTEP [2], a plane-wave basis set package for computing the energies and DOS for our system. Our system is the system depicted in Figure 1, a unit cell of silver atoms on a fcc lattice, with lattice parameter a = 4.086 angstroms. The GGA PBE functional [3] was used for all calculations. The SCF tolerance is 5 x 10^-7 eV, the Koelling-Harmon treatment is used for the atomic solutions [4], and ultrasoft pseudopotentials are generated on the fly. The atoms used in the pseudopotential are 4s2 4p6 4d10 5s1. The integration for the DOS is performed using the 0.2 smearing.

Results

First, a convergence study was performed for the cutoff energy. The results are shown in Figure 2.

Figure 2: Convergence of the system energy with respect to cutoff energy

Above a cutoff energy of 600 eV, the system energy is converged to the order of tenths of an eV, specifically to 0.15 eV. We select 600 eV as the cutoff energy for the rest of our calculations due to being the most computational efficient value converged to this degree. We then calculate both the system energy and the density of states for different k-point meshes. Because our unit cell is cubic, all k-point meshes are of the form NxNxN, where N is an integer, and we will characterize each mesh simply by the number N. Calculations were performed for N=6 to N=20, covering all possible values.

Figure 3: Convergence of the system energy with respect to k-point number

The relevant feature of the convergence results for the system energy with respect to the k-point mesh is that the variance from the 7x7x7 mesh results is never greather than 0.05 eV for higher N, which is already more converged than our results for the cutoff energy. So we can say that our system energy is sufficiently converged at N = 7.

Before comparing the convergence of the DOS, we need to determine a metric for the convergence. It is clear that the behavior of the DOS should qualitatively change in a certain way with an increasing number of k-points, but it is harder to quantitatively measure this change than the change in system energy. We can no longer compare the change in a single value, as the calculated DOS is a curve for each k-point mesh, as show in Figures 4 and 5.

To compare these curves, we take the difference of the DOS value at every point on the curve. We then square these differences, and take the average for each curve. This quantity, defined as (1/N)*Sum(DOS(k, i) – DOS(k-1,i)) is the mean squared deviation between the two DOS curves. Here, N is the number of points on the curve, k or k-1 is the number of points in the k-point mesh used in that DOS calculation, and i is the eV value of a specific point on the curve. The mean squared deviation is a measure of how different two curves are, and minimizing this quantity will imply that the DOS curves are more similar.

Figure 4: The DOS for silver, calculated with a 6x6x6 k-point mesh at 600 eV. The 0 of the x-axis is the Fermi energy.

Figure 5: The DOS for silver, calculated with a 20x20x20 k-point mesh at 600 eV. The 0 of the x-axis is the Fermi Energy.

Because of the extent of the DOS, much of which has no states, we restrict our analysis from eight eV below the Fermi energy to eight eV above the Fermi energy. This portion of the DOS curves are plotted in figures 6 and 7.

Figure 6: The DOS around the Fermi energy for silver, calculated with a 6x6x6 k-point mesh

Figure 7: The DOS around the Fermi energy for silver, calculated with a 20x20x20 k-point mesh

As expected, when the DOS is calculated with a denser sampling of k-space, the curve is smoother and more detailed, rather than sharply varying due to poor sampling.

Quantitatively, as we move from one value of N to the next, the average change in a single point on the curve decreases. It is noteworthy for the some steps in N, there is no change at any point on the DOS curve. This is not indicative of the DOS being fully converged, as shown by later changes to the average deviation. This is representative of the sampling of the new k-point mesh not being significantly different from the previous mesh. The results of this analysis are shown in figures 8 and 9.

Figure 8: Average Squared Difference for the DOS curves with k-mesh specified by N and N-1

Figure 9: Data from figure 8, re-scaled to be more visible

It is noteworthy that these curves are not quite the same as those given in figures 2 and 3. In the case of the system energy, we can simply plot the resulting system energy against k-point number and see the overall trend of the values. We can measure the convergence and fit for behavior at large numbers of k-points. In contrast, we cannot make a truly analogous comparison for the DOS curves.

When we speak about the convergence of the system energy with respect to k-point number, we are addressing the variation between calculation results at a given number of k-points and the result we would obtain with an infinite number of k-points. For a variational parameter, such as the variation of system energy with respect to cutoff energy, this is a relatively easy comparison, since the change in the system energy will always decrease for every additional increase of the energy cutoff by a given amount. For the case of the variation of system energy with respect to k-points, we often test a computationally efficient range of k-points and show that after a certain number of k-points, the system energy does not vary beyond a certain interval, such as 0.05 eV.

In the case of the DOS in this post, we use the mean squared deviation between subsequent curves as a method analogous to comparing variation in the system energy from k-point number to k-point number. The analogous convergence condition is then that the mean squared deviation above a certain number of k-points is below a certain threshold. To facilitate this comparison, we plot the sum of the mean squared deviation in figure 10.

Figure 10: Sum of mean squared deviations for DOS curves. Summed from N=20 back to N=6

Here, we compare the summed mean squared deviation against N=20, our most computationally intensive calculation. This figure makes it clear how the convergence of our calculation proceeds, but we cannot extract predictions about the behavior of the DOS curve at arbitrarily high number of k-points from this figure.

Conclusions

As predicted by the theory, we need many more k-points to appropriately converge the DOS results than we do for the system energy. This should be clear from the clear qualitative differences between figures 6 and 7. More quantitatively, the sum of the mean squared deviations of the DOS curves is approximately 0.005 from N=13 onward. From N=18 onward, the sum of the mean squared deviations is 0.001. If the DOS curve itself is being reported, then the calculation performed at 13x13x13 k-points would give a well-converged result, whereas a lower number of k-points would leave the the sum of the mean squared deviations above 0.18. Meanwhile, the convergence of the system energy with respect to k-points needs only a 6x6x6 k-point mesh.

Additional work in this topic area can investigate the area in an even narrower region around the Fermi energy, which would require an even denser sampling of k-space in that region. This part of the density of states is of particular interest for thermodynamic properties, but is beyond the scope of this post.

References

[1]  D. Sholl and J. Steckel, Density Functional Theory: A Practical Introduction. (Wiley 2009)

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

[3]  John P. Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized Gradient Approximation Made Simple”, Phys. Rev. Lett. 77, 3865 – Published 28 October 1996; Erratum Phys. Rev. Lett. 78, 1396 (1997)

[4]  D D Koelling and B N Harmon 1977 J. Phys. C: Solid State Phys. 10 3107