Hybrid Functional Calculations

by Ute Werner for exciting boron-10

Purpose: In this tutorial, you will learn how to perform band-structure calculations using hybrid functionals. As an example, the band structure of diamond is calculated with the PBE0-functional.

---

[See here: Theoretical Background]

Theoretical Background

In hybrid functionals, the exchange potential $v_{\mathrm{x}}^{\mathrm{DFT}}$ is partially replaced by the exact-exchange counterpart $v_{\mathrm{x}}^{\mathrm{exact}}$. In the simpliest one parameter form, the resulting exchange-correlation potential $v_{\mathrm{xc}}^\mathrm{HYB}$ can be formulated as:

(1)

\begin{align} v_{\mathrm{xc}}^\mathrm{HYB}=\left(1-\alpha\right)v_{\mathrm{x}}^\mathrm{DFT}+\alpha_{\mathrm{x}} v_{\mathrm{x}}^\mathrm{exact}+v_{\mathrm{c}}^\mathrm{DFT}, \end{align}

where $\alpha_{\mathrm{x}}$ is the mixing coefficient for the exact exchange and $v_{\mathrm{c}}^{\mathrm{DFT}}$ the DFT correlation. In standard implementations, GGA functionals are employed for the DFT part while Hartree-Fock (HF) is used for the exact exchange. Alternatively, the local exact-exchange potential $v_{\mathrm{x}}^{\mathrm{EX}}$ as obtained from the optimized effective-potential (OEP) approach can be used for the construction of hybrid functionals.
In this tutorial, we will use the PBE0 functional defined as:

(2)

\begin{align} v_{\mathrm{xc}}^\mathrm{PBE0}=\left(1-\alpha\right)v_{\mathrm{x}}^\mathrm{PBE}+\alpha_{\mathrm{x}} v_{\mathrm{x}}^\mathrm{exact}+v_{\mathrm{c}}^\mathrm{PBE}. \end{align}

We will use the JCP exchange and the parameter $\alpha_{\mathrm{x}}=0.25$ following the original papers from Perdew et al. (JCP 1996) and Adamo et al. (JCP 1999).

---

1. PBE0 calculation

We start by creating the working directory for this tutorial and move inside of it.

$ mkdir Diamond_PBE0
$ cd Diamond_PBE0

i) Preparation of the input file

The inputfile input.xml for the ground-state and band-structure calculation of diamond using PBE0 could look like this

<input>
 
   <title>Diamond PBE0</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"
      nempty="50"
      xctype="HYB_PBE0">
      <Hybrid 
         exchangetype="HF" 
         excoeff="0.25" />
   </groundstate>
 
   <properties>
      <bandstructure>
         <plot1d>
            <path steps="200">
              <point coord=" 0.750   0.500   0.250" label="W" />
              <point coord=" 0.500   0.500   0.500" label="L" />
              <point coord=" 0.000   0.000   0.000" label="GAMMA" />
              <point coord=" 0.500   0.500   0.000" label="X" />
              <point coord=" 0.750   0.500   0.250" label="W" />
              <point coord=" 0.750   0.375   0.375" label="K" />
            </path>
         </plot1d>
      </bandstructure>
   </properties>
 
</input>

Let's have a closer look to the parameters in this file.

Parameter Description xctype = "HYB_PBE0" The calculation is performed using the hybrid functional PBE0 nempty = "50" The number of empty bands, crucial for the convergence exchangetype = "HF" You can choose between HF and DFT-exact exchange excoeff = "0.25" Value for the amount of exact exchange. For the example the standard value for the PBE0 functional is used

For the meaning of the different arguments and keywors of the properties section have a look to Electronic-structure calculation.
Since the PBE0 calculations are much more time consuming than corresponding LDA or GGA ones, computational parameters have been chosen in a way to speed up the calculation. For further details on the parameters see Input Reference.

ii) Execute the calculation

Be sure to set the correct path for the exciting root directory (indicated in this example by $EXCITINGROOT) to the one pointing to the place where the exciting directory is placed. In order to do this, use the command

$ SETUP-excitingroot.sh

Start the calculation by the shell command

$ time excitingser >& output.txt &

---

2. PBE calculation

In order to investigate the effect of hybrid functionals on the band structure, we will also perform a GGA calculation. For this purpose, we will use the PBE functional. Again, we start by creating a new working directory and move inside.

$ mkdir Diamond_PBE
$ cd Diamond_PBE

i) Preparation of the input file

You can copy the inputfile input.xml from the previous calculation. In order to use PBE, we modify the functional defined by the keyword xctype:

...
   <groundstate 
      ngridk="4 4 4"
      nempty="50"
      xctype="GGA_PBE">
   </groundstate>
...

ii) Execute the calculation

The calculation can be started exactly in the same way as the PBE0 one.

---

3. Post processing

For a graphical comparison of the band structures calculated with the PBE0 and PBE functionals you can use the script PLOT-PBE0bands.py. Move to the parent directory (the one where you can find the PBE0 and PBE subdirectories) and execute the script. It will produce the following output on the screen.

$ PLOT-PBE0bands.py

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

    Enter the working directories 

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

name of PBE0 directory: Diamond_PBE0
name of PBE directory: Diamond_PBE

$

In this example, the first input entry is preceded by "name of PBE0 directory: ". This entry value must be typed on the screen and must be the name of the directory which contains your PBE0 results. The second entry should then be the PBE directory. The script PLOT-PBE0bands.py produces a PostScript (bandstructures.eps) and a PNG (bandstructures.png) output file, as well as an image on the screen. Using the examples above you will obtain the following plot.

Note, that for PBE0 the band structur has to be interpolated. Thus, for a nice band structur plot you will need a larger k-grid. A PBE0 calculation for diamond with an 8$\times$8$\times$8 k-grid and 100 empty states will look like this:

---

4. Excercises

1. Try to modify some of the parameters such as nempty or ngridk in the input files and investigate the effect on the band structure and on the computational time.
2. Modify the keyword for the exact exchange type in order to use the exact exchange from DFT: exchangetype = "OEP". Compare the results to the previous PBE0 calculations. HINT: If you use the OEP potential, you can also modify the EXX/OEP parameters as explained in the tutorial Exact Exchange Calculations.

Literature

PBE0: Perdew, J. P. and Ernzerhof, M. and Burke, K., J. Chem. Phys. 22, 9982 (1996); Adamo, C. and Barone, V., J. Chem. Phys. 110, 6158, (1999).