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 a0-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 .



### 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)