Hybrid-Functional Calculations

by Cecilia Vona, Ute Werner, & Dmitrii Nabok for exciting oxygen

Purpose: In this tutorial, you will learn how to perform ground-state calculations with the hybrid functionals PBE0 and HSE. We will use diamond as test material. Since Diamond has an indirect band-gap, instead of considering the gap, we will look at how the transition energies in specific symmetry point changes, depending on the functional employed. Specifically we will consider the transition Γ→Γ and Γ→X.

---

1. Theoretical background

The hybrid exchange-correlation functionals are constructed by substituting a fraction of local/semi-local exchange by a fraction of non-local exact-exchange:

(1)

\begin{align} E_{\mathrm{xc}}^\textrm{hyb}=E_{\mathrm{xc}}^{\mathrm{L}}+\alpha\, \left(E_{\mathrm{x}}^{\mathrm{NL}}- E_{\mathrm{x}}^{\mathrm{L}}\right). \end{align}

The non-local exact exchange is often computed employing Hartree-Fock (HF). The type of the DFT local/semi-local functional and the fraction of non-local exact exchange substituted ($\alpha$) vary depending on the hybrid functional type.

The hybrid functionals PBE0 and HSE are implemented in exciting. PBE0 has the following form:

(2)

\begin{align} E_{\mathrm{xc}}^{\mathrm{PBE0}}= E_{\mathrm{xc}}^{\mathrm{PBE}}+\alpha\,\left(E_{\mathrm{x}}^{\mathrm{HF}}-E_{\mathrm{x}}^{\mathrm{PBE}}\right), \end{align}

in which the semi-local functional is PBE and the fraction of exact-exchange substituted is uniquely determined by the mixing parameter $\alpha$, which has the value 0.25.

The hybrid functional HSE, which is a screened hybrid functional, has the following form:

(3)

\begin{align} E_{\mathrm{xc}}^{\mathrm{HSE}}= E_{\mathrm{xc}}^{\mathrm{PBE}}+\alpha\left[E_{\mathrm{x}}^{\mathrm{HF,SR}}(\omega)-E_{\mathrm{x}}^{\mathrm{PBE,SR}}(\omega)\right]. \end{align}

Also in HSE, PBE is used as semi-local functional, but only the short-range part of the non-local exact exchange is substituted. In fact, in HSE the amount of non-local exact exchange is determined by both the mixing parameter $\alpha$=0.25 and the screening screening parameter $\omega$, which defines the cutoff between the short-range and the long-range part of the Coulomb operator through the error function and its complementary:

(4)

\begin{align} v(r)=v^{\mathrm{SR}}(r,\omega)+v^{\mathrm{LR}}(r,\omega)= \frac{\mathrm{erfc}(\omega r)}{r}+\frac{\mathrm{erf}(\omega r)}{r}. \end{align}

The version of HSE which is implemented in exciting is HSE06, in which the PBE and HF contributions have the same value of $\omega$=0.11 a₀^-1.

---

2. PBE calculation

In order to investigate and compare the effect of hybrid functionals on the transition energies, we will start by performing a GGA calculation. For this purpose, we will use the PBE functional.

To prepare your calculation, create a new, empty directory named diamond-hybrids and move inside it.

$ mkdir diamond-hybrids
$ cd diamond-hybrids

i) Preparation of the input file

Inside the current directory, create the input file input.xml, to compute the ground-state of diamond using the semi-local functional PBE. The file can look like this:

<input>
 
   <title>Diamond PBE</title>
 
   <structure speciespath="$EXCITINGROOT/species" >
      <crystal scale="6.7425">
         <basevect> 0.0     0.5     0.5 </basevect>
         <basevect> 0.5     0.0     0.5 </basevect>
         <basevect> 0.5     0.5     0.0 </basevect>
      </crystal>
      <species speciesfile="C.xml">
         <atom coord="0.00 0.00 0.00" />
         <atom coord="0.25 0.25 0.25" />
      </species>
   </structure>
 
   <groundstate 
      ngridk="4 4 4"
      rgkmax="5.0"
      xctype="GGA_PBE">
   </groundstate>
 
</input>

ii) Execute the calculation

Before to start the calculation, remember to set the correct path for the exciting root directory ($EXCITINGROOT) to the one pointing to the place where the exciting directory is placed, using the command

$ SETUP-excitingroot.sh

Start the calculation by using the script EXECUTE-single.sh followed by the DIRECTORYNAME.

$ EXECUTE-single.sh PBE

The calculation will be executed inside the subdirectory xc-rundir which at the end of the run will be renamed PBE.

---

3. PBE0 calculation

Now we can repeat the calculation by using hybrid functionals. We will start with the hybrid functional PBE0.

i) Preparation of the input file

To compute the ground state of diamond with the hybrid functional PBE0 it is necessary to modify the groundstate element and, if you like, you can modify/adapt the content of the title element as in the example

<input>
 
<title>Diamond PBE0</title>
 
...   
 
   <groundstate 
      ngridk="4 4 4"
      rgkmax="5.0"
      nempty="20"
      xctype="HYB_PBE0">
      <Hybrid  
         excoeff="0.25" />
   </groundstate>
 
</input>

The first thing to notice is the presence of the new element Hybrid inside the element groundstate. This element is not mandatory to perform calculations with hybrid functionals. If it is not included, all the parameters specified by the element Hybrid will be initialized with their default values. Moreover, by comparing the input with the input of PBE, we observe that in the example is specified the number of empty states (nempty), which is crucial for having convergent hybrids calculations.

To the parameters introduced in the ground-state calculation with hybrid functionals will be discussed later. Before we do that, execute the calculation with PBE0 and create the HSE input file.

ii) Execute the calculation

As well as in the case of PBE, to perform the calculation use the script EXECUTE-single.sh but this time the DIRECTORYNAME will be PBE0 .

$ EXECUTE-single.sh PBE0

---

4. HSE calculation

Now we repeat everything using the hybrid functional HSE.

i) Preparation of the input file

Start by modifying again the groundstate element in the file input.xml, inside the directory diamond-hybrids, as in the example:

<input>
 
<title>Diamond HSE</title>
 
... 
 
   <groundstate 
      ngridk="4 4 4"
      rgkmax="5.0"
      nempty="20"
      xctype="HYB_HSE">
      <Hybrid 
         excoeff="0.25" 
         omega="0.11" />
   </groundstate>
 
</input>

You have probably noticed, that the running time slightly increases when HSE is employed instead of PBE0. This can seem weird, since HSE is known for being faster. For a fix number of k-points this is not the case in exciting, due to the LAPW+lo basis. However, convergent HSE calculations are usually faster than convergent PBE0 calculations, since the latter converge considerably slower with respect to the number of k-points.

---

☛ Plot k-points convergence

In these plots are shown the difference of the transition energies between Γ→Γ and Γ→X, computed with the PBE, PBE0 and HSE xc-functionals, for different k-points grid. The reference calculation is computed for $N_{\mathbf{k}}=10^3$. From this you can have a better understanding of the k-points convergence of the different functionals.

To repeat the calculations with PBE0 and HSE, for higher number of k-points, requires too much time for the purpose of this tutorial. You can try to reproduce these plots as exercise.

---

We can now have a closer look to the most relevant parameters in ground-state calculations with hybrid functionals.

Parameter Description xctype = "HYB_PBE0" or "HYB_HSE" The type of (hybrid) functional used needs to be specified in xctype. In exciting the hybrid functionals available are PBE0 or HSE and the one shown are the respectively keywords. nempty = "20" The number of empty states is crucial for the convergence and is dependent from rgkmax. excoeff = "0.25" Keyword, inside the element Hybrid, for the mixing parameter $\alpha$. The default value is 0.25. omega = "0.11" Keyword, inside the element Hybrid, for to screening parameter $\omega$ in units of a0-1 (only with HSE). The default value is 0.11 a0-1.

For further details on the parameters see Input Reference.

ii) Execute the calculation

Also in this case to run the calculation in the directory HSE use EXECUTE-single.sh .

$ EXECUTE-single.sh HSE

---

5. Post processing

To observe some of the effects that the different functionals produce on the electronic structure, we compare the transition energies in specific high-symmetry points. Specifically, we will consider the following transition Γ→Γ and Γ→X.

To do it, after verified to be in the parent directory (diamond-hybrid), run the script COMPARE-transition-energies.py. The script will ask you to entry the names of the directories where the calculations have been performed and the output will look like this

$ COMPARE-transition-energies.py

################################################

 Enter the names of the working directories and 
 enter "(Q)uit" or "(q)uit" to terminate input.

 Directory name ==> PBE
 Directory name ==> PBE0
 Directory name ==> HSE
 Directory name ==> q

------------------------------------------------

 Transition energies in eV:

                    PBE       PBE0       HSE
 Gamma -> Gamma:    5.632     7.886     7.098
 Gamma -> X    :    4.812     6.842     6.047

################################################
$

Do not compare the values obtained in this tutorial with other references, since due to the considerable amount of time needed to compute the ground-state with hybrid functionals, which is much more in comparison with the time needed with LDA or GGA functionals, the parameters used have been chosen to speed up the calculations and are far from the convergent one.

However, we can still notice that the transition energies computed with the hybrid functionals are higher than the one computed with PBE, which is known to highly underestimates them. Moreover, PBE0 slightly overestimates the gap of semi-conductors, and therefore the transition energies, while HSE is able to better reproduce them. Looking at the values obtained, we notice that the transition energies computed with HSE are between the one computed with PBE and PBE0, and closer to the latter.

If you want to analyse more transition energies than the one taken into account in this tutorial, you can extract them from the output file EIGVAL.OUT. Remember that the energies in the file are in Ha.

---

6. Band structure and density of states

When hybrids functionals are employed, to compute properties, such as band structure and density of states, the Wannier-functions interpolation schema needs to be adopted. To use the Wannier interpolation schema is not straightforward and therefore there is a tutorial dedicated Wannier Functions for Interpolation in Reciprocal Space.

In the tutorials the properties are computed only for the hybrid functional PBE0, as exercise you can try to repeat the Wannier-functions tutorial for the hybrid functional HSE.

---

7. Exercise

Try to modify some of the parameters such as rgkmax, nempty, or ngridk in the input files. Remember that the optimal value of nempty depends on the value of rgkmax. For more complex systems nempty can be quite big.

In order to get an hint about the maximum number of empty states that can be considered, use the following command.

$ grep "Maximum Hamiltonian size" INFO.OUT

---

Literature

• C. Adamo and V. Barone, J. Chem. Phys. 110, 6158, (1999)
• M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. 110, 5029 (1999)
• A. V. Krukau et al., J. Chem. Phys. 125, 224106 (2006)
• M. Schlipf et al., Phys. Rev. B 84, 125142 (2011)