**Part 1: An overview of available functionals**¶

By the end of this tutorial, you will be able to:

- explain what is a hybrid functional
- setup a PBE0 calculation
- compute the band gap

**1.1 Task**¶

*Compute the band gap of cubic diamond Si using PBE, and PBE0.*

Whether the Kohn–Sham (KS) orbitals, $\psi_{n\mathbf{k}}(\mathbf{r})$, reflect the right physics crucially depends on the choice of the exchange-correlation (xc) functional. The most basic possibilities are (i) the local density approximation (LDA), or (ii) the semi-local generalized gradient approximation (GGA). Yet, there are different parametrizations to construct LDA functionals, and some more freedom in constructing GGA functionals. The PBE functional is one of many implemented GGA functionals in VASP, and perhaps the most commonly used.

Hybrid functionals make up a special category of xc functionals, that mix some amount of Fock exchange into a density functional. This approach is inspired by the well-known Hartree–Fock (HF) method, which has an exact solution for uncorrelated systems. Regarding condensed matter, it aims to cure the infamous underestimation of band gaps by density functional theory (DFT). The nonlocal Fock potential reads
$$\tag{1}
V_{\mathrm{x}}\left(\mathbf{r},\mathbf{r}'\right)=
-\frac{e^2}{2}\sum_{n\mathbf{k}}f_{n\mathbf{k}}
\frac{\psi_{n\mathbf{k}}^{*}(\mathbf{r}')\psi_{n\mathbf{k}}(\mathbf{r})}
{\vert \mathbf{r}-\mathbf{r}' \vert},
$$
where $f_{n\mathbf{k}}$ is the occupation of band $n$ and $\vec{k}$ point $\mathbf{k}$, $\mathbf{r}$ is the position and $e$ the charge of an electron. Equation (1) does not depend on the density, but on the KS orbitals and thus, hybrid functionals are *not* density functionals.

PBE0 is one of the hybrid functionals available in VASP and has the following composition: $$\tag{2} E_{\mathrm{xc}}^{\mathrm{PBE0}} = \frac{1}{4} E_{\mathrm{x}}^{\mathrm{HF}} + \frac{3}{4} E_{\mathrm{x}}^{\mathrm{PBE}} + E_{\mathrm{c}}^{\mathrm{PBE}}. $$ Here, $E_{\mathrm{x}}^{\mathrm{HF}}$ is the exchange energy, $E_{\mathrm{x}}^{\mathrm{PBE}}$ is the PBE exchange energy and $E_{\mathrm{c}}^{\mathrm{PBE}}$ is the PBE correlation energy. Compared to others, PBE0 uses a nonempirical justification for the mixing ratio, and the Fock potential enters without range-dependent screening.

**1.2 Input**¶

The input files to run this example are prepared at `$TUTORIALS/hybrids-part1/e01_Si-gap`

. Check them out!

```
cubic diamond Si
5.46873
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
2
Direct
-0.125 -0.125 -0.125
0.125 0.125 0.125
```

```
SYSTEM = cd Si PBE
GGA = PE
ISMEAR = -5
LORBIT = 11
```

PBE0/INCAR

```
SYSTEM = cd Si PBE0
GGA = PE
LHFCALC = T
ISMEAR = 0
SIGMA = 0.01
ALGO = Damped
TIME = 0.8
```

```
K-Points
0
Gamma
7 7 7
0 0 0
```

*Pseudopotential of Si.*

The POSCAR file defines the cubic diamond structure of Si.

The INCAR files specify explicitly which xc functional is used. Check the default setting for the GGA tag! Is it necessary to set it in this calculation? PBE0 is activated by LHFCALC, as long as related tags, i.e., AEXX, AGGAC and ALDAC, have their default values. How do these tags relate to Equation (1)?

In the PBE calculation, ISMEAR=-5 switches on the tetrahedron method with Blöchl corrections to determine the partial occupations $f_{n\vec{k}}$. This smearing yields the smoothest density of states. LORBIT=11 enables that the DOS is written to the DOSCAR file and the **vaspout.h5** file.

In the PBE0 calculation, we use a damped algorithm for the electronic relaxation, where the total energy must be variational. That is, because generally direct optimization is more robust than iterative methods. Why do we not immediately compute the DOS, like in the PBE calculation?

## Click to see the answer!

The Blöchl corrections causes energies not to be variational, so that we use Gaussian smearing with ISMEAR=0 and SIGMA=0.01 instead. After convergence, a single-shot calculation with ICHARG=11 , ISMEAR=-5 and LORBIT=11 can be used to obtain the DOS.

The POTCAR file contains the pseudopotential generated for PBE, and the KPOINTS file defines an equally spaced 𝝘-centered grid of $\vec{k}$ points.

The overall settings are such that the precision lies within 0.1eV of the band gap. What would you need to do, to confirm this statement and how can you improve the precision?

## Click to see the answer!

You would need to perform convergence studies of the cutoff energy and the number of $\vec{k}$ points with respect to the band gap. In other words, you would loop over different settings and judge how much the value of the band gap changes.

**1.3 Calculation**¶

Open a terminal, navigate to this example's directory and run VASP by entering the following:

```
cd $TUTORIALS/hybrids-part1/e01_*
vasp_std
```

You can plot the DOS using py4vasp for the PBE calculation.

```
import py4vasp
pbe_calc = py4vasp.Calculation.from_path("./e01_Si-gap")
pbe_calc.dos.plot()
```

Where in the DOS plot can you see the band gap?

## Click to see the answer!

The electronic band gap is the energy difference between the lowest unoccupied molecular orbital (LUMO) and the highest occupied molecular orbital (HOMO), where in case of DFT molecular orbitals are substituted by KS orbitals. In the DOS plot, that is the energy range without any states in the vicinity of the Fermi energy, which is set to zero here.

Next, copy the WAVECAR file to the directory of the PBE0 calculation and run VASP there by entering the following:

```
cd $TUTORIALS/hybrids-part1/e01_*/PBE0
cp ../WAVECAR .
mpirun -np 2 vasp_std
```

The **stdout** printed during the calculation shows an additional line with information related to the damped algorithm at each step of the electronic relaxation.

## Click to see the answer!

No, generally total energies between different xc functionals cannot be compared. Only energy differences can be compared. So, it is possible to compare the band gap.

```
import py4vasp
import numpy as np
def gap(calc):
homo = np.amax(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] > 0.5])
lumo = np.amin(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] < 0.5])
return lumo - homo
pbe_calc = py4vasp.Calculation.from_path("./e01_Si-gap")
pbe0_calc = py4vasp.Calculation.from_path("./e01_Si-gap/PBE0")
print("Si bandgap\nPBE:", gap(pbe_calc))
print("PBE0:", gap(pbe0_calc))
```

Which band gap is larger?

## Click to see the answer!

The PBE0 band gap is larger! This is a general trend and aims to compensate for the underestimation of the band gap within DFT compared to the experimental observation in many materials.

To visualize the PBE0 gap, you can again compute the DOS. Open the INCAR file of your PBE0 calculation, change ISMEAR to `-5`

and add LORBIT=11, and ICHARG=11. Then, restart the calculation with the existing CHGCAR file by entering the following:

```
cd $TUTORIALS/hybrids-part1/e01_*/PBE0
mpirun -np 2 vasp_std
```

```
import py4vasp
pbe0_calc = py4vasp.Calculation.from_path("./e01_Si-gap/PBE0")
pbe0_calc.dos.plot()
```

Did the overall appearance of the DOS change much compared to the PBE calculation?

## Click to see the answer!

The DOS is renormalized, which affects the relative energy between bands, but also the band width of a single band. These subtle changes can crucially influence the properties of the material.

By the end of this tutorial, you will be able to:

- perform an HF-like calculation using exact exchange (EXX)
- use the B3LYP functional
- explain the limitations of HF and B3LYP

**2.1 Task**¶

*Compute the band gap of Ar using PBE, B3LYP and HF.*

B3LYP is an example of a semiempirical hybrid functional, that aims to describe molecular properties. Its parameters are chosen to reproduce experimental energies, such as the atomization energies, electron-proton affinities and ionization potentials, of molecules and atoms in a test set. This involves both levels, the underlying density functional and the amount of mixing with Fock exchange and reads as follows:

$$\tag{2 a} E_{\mathrm{xc}}^{\mathrm{B3LYP}} = E_{\mathrm{x}}^{\mathrm{B3LYP}} + E_{\mathrm{c}}^{\mathrm{B3LYP}} $$$$ \tag{2 b} E_{\mathrm{x}}^{\mathrm{B3LYP}} = 0.8 E_{\mathrm{x}}^{\mathrm{LDA}} + 0.2 E_{\mathrm{x}}^{\mathrm{HF}} + 0.72 \Delta E_{\mathrm{x}}^{\mathrm{B88}}. $$$$ \tag{2 c} E_{\mathrm{c}}^{\mathrm{B3LYP}} = 0.19 E_{\mathrm{c}}^{\mathrm{VWN3}} + 0.81 E_{\mathrm{c}}^{\mathrm{LYP}}. $$On the other hand, one could use the entire Fock exchange and perform an *HF-like calculation*. As the HF energy is evaluated on KS orbitals, this approach is termed exact exchange (EXX) method. Check out the corresponding tag: AEXX.

One fundamental flaw of the HF method, which is inherited by B3LYP and EXX, is that it fails in the limit of the homogeneous electron gas, as the Fourier transform of the Fock potential develops a singularity at $\mathbf{q}=0$. This causes a complete suppression of itinerant states. In other words, PBE0, B3LYP and EXX are not recommended for systems with significant itinerant character, such as metals and small gap semiconductors, see Paier *et al.* (2007).

Here, let us investigate the band gap of Ar, which is found to have a large band gap 14.4 eV in experiment!

**2.2 Input**¶

The input files to run this example are prepared at `$TUTORIALS/hybrids-part1/e02_Ar-gap`

. Check them out!

```
Ar Fm-3m
1.0
3.454253 0.000000 1.994314
1.151418 3.256701 1.994314
0.000000 0.000000 3.988628
Ar
1
direct
0.000000 0.000000 0.000000 Ar
```

```
SYSTEM = Ar
ISMEAR = -5
ENCUT = 300
LORBIT = 11
```

HF/INCAR

```
SYSTEM = Ar HF
LHFCALC = T
AEXX = 1
ISMEAR = 0
SIGMA = 0.01
ENCUT = 300
ALGO = Damped
TIME = 0.4
```

B3LYP/INCAR

```
SYSTEM = Ar B3LYP
GGA = B3
LHFCALC = T
AEXX = 0.2
AGGAX = 0.72
AGGAC = 0.81
ALDAC = 0.19
ISMEAR = 0
SIGMA = 0.01
ALGO = Damped
TIME = 0.4
```

```
K-Points
0
Gamma
4 4 4
0 0 0
```

*Pseudopotential of Ar.*

The POSCAR file defines the face centered cubic Ar crystal.

The first INCAR file, starts a PBE calculation. Check the meaning of unfamiliar tags on the VASP Wiki!

In the HF INCAR file, the evaluation of Fock exchange is switched on using the LHFCALC tag and the fraction of EXX set to AEXX. So only Fock exchange and no density functional exchange is used. Recall the meaning of the remaining tags! What are the default values of GGA, AGGAX, AGGAC and ALDAC here?

**2.3 Calculation**¶

Open a terminal, navigate to this example's directory and run VASP by entering the following:

```
cd $TUTORIALS/hybrids-part1/e02_*
vasp_std
```

You can visualize the DOS of the PBE calculation with py4vasp!

```
import py4vasp
pbe_calc = py4vasp.Calculation.from_path("./e02_Ar-gap")
pbe_calc.dos.plot()
```

```
import py4vasp
import numpy as np
def gap(calc):
homo = np.amax(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] > 0.5])
lumo = np.amin(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] < 0.5])
return lumo - homo
pbe_calc = py4vasp.Calculation.from_path( "./e02_Ar-gap")
hf_calc = py4vasp.Calculation.from_path( "./e02_Ar-gap/HF")
b3lyp_calc = py4vasp.Calculation.from_path( "./e02_Ar-gap/B3LYP")
for calc in [pbe_calc, hf_calc, b3lyp_calc]:
print(str(calc.path()).split("/")[-1], gap(calc))
```

## Click to see the answer!

After executing the following

```
cd $TUTORIALS/hybrids-part1/e02_*/B3LYP
cp ../WAVECAR .
mpirun -np 2 vasp_std
```

you can execute the python code above. While the PBE band gap of 8.5eV is much smaller than the experimental value of 14.4 eV, HF 17.6 eV and B3LYP 10.3 eV both also miss by more than 3 eV.

**2.4 Questions**¶

- Which systems does the PBE0 and B3LYP functional and exact-exchange method fail to describe well? Why?
- What does the
**AGGAC**tag control?

By the end of this tutorial, you will be able to:

- explain the basic assumptions taken by the HSE functional
- optimize the HSE functional with respect to the experimental band gap

**3.1 Task**¶

*Optimize the amount of Fock exchange included to reproduce the experimental band gap of MgO with HSE functionals.*

The HSE functional HSE06 is a semiemperical functional that decomposes the $1/r$ dependence in the Fock exchange into a short-range (SR) and a long-range (LR) contribution at the cost of a free screening parameter $\mu$:

$$\tag{3} E_{\mathrm{xc}}^{\mathrm{HSE}} = \frac{1}{4} E_{\mathrm{x}}^{\mathrm{HF,SR}}(\mu) + \frac{3}{4} E_{\mathrm{x}}^{\mathrm{PBE,SR}}(\mu) + E_{\mathrm{x}}^{\mathrm{PBE,S LR}}(\mu) + \frac{1}{4} E_{\mathrm{c}}^{\mathrm{PBE}}. $$In SR, it is identical to the PBE0 functional, but the LR singularity in the Fock exchange, that causes the suppression of all metallic states, is avoided by introducing orbital-dependent screening. Additionally, the computational cost can be reduced by downsampling. The screening parameter is the empirically determined to have a value of 0.2 - 0.3 Å$^{−1}$ and controlled by HFSCREEN in VASP. It renders the interaction range, 2/$\mu \approx$ 10Å, to be typically over a few nearest neighbors.

The gap size still crucially depends on the amount of Fock exchange introduced, and less strongly on the electronic screening. Let us take a practical approach here, and optimize the mixing to reproduce the experimental band gap of MgO, that is 7.7 eV! Strictly speaking, this is diverting from a first-principles philosophy, and also not following exactly the recipe prescribed by HSE.

**3.2 Input**¶

The input files to run this example are prepared at `$TUTORIALS/hybrids-part1/e03_MgO-gap-optimization`

. Check them out!

```
MgO Fm-3m (No. 225)
1.0
2.606553 0.000000 1.504894
0.868851 2.457482 1.504894
0.000000 0.000000 3.009789
Mg O
1 1
direct
0.000000 0.000000 0.000000 Mg
0.500000 0.500000 0.500000 O
```

```
SYSTEM = MgO
ENCUT = 500
ISMEAR = -5
LORBIT = 11
```

HSE06/INCAR

```
SYSTEM = MgO HSE06
ENCUT = 500
ISMEAR = 0
SIGMA = 0.01
LHFCALC = T
HFSCREEN = 0.2
AEXX =
ALGO = Damped
TIME = 0.4
```

```
K-Points
0
Gamma centered
4 4 4
0 0 0
```

*Pseudopotential of Mg, and O.*

**3.3 Calculation**¶

Open a terminal, navigate to this example's directory and run VASP by entering the following:

```
cd $TUTORIALS/hybrids-part1/e03_*
vasp_std
```

The band gap with PBE can be seen in the DOS. Plot it using py4vasp!

```
import py4vasp
pbe_calc = py4vasp.Calculation.from_path("./e03_MgO-gap-optimization")
pbe_calc.dos.plot()
```

Then, in the `HSE06`

directory run a loop over different values of AEXX by entering the following:

```
cd $TUTORIALS/hybrids-part1/e03_*/HSE06
bash loop.sh
```

**loop.sh**

```
for aexx in 0.45 0.50 0.55 0.60
do
vasp_rm
cp ../WAVECAR .
cat > INCAR << EOF
SYSTEM = MgO HSE06
ENCUT = 500
ISMEAR = 0
SIGMA = 0.01
LHFCALC = T
HFSCREEN = 0.3
AEXX = $aexx
ALGO = Damped
TIME = 0.4
EOF
mpirun -np 2 vasp_std
cat OUTCAR | grep "Total CPU time used"
mkdir aexx$aexx
cp vaspout.h5 aexx$aexx/.
done
```

```
import py4vasp
import numpy as np
def gap(calc):
homo = np.amax(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] > 0.5])
lumo = np.amin(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] < 0.5])
return lumo - homo
for aexx in ['0.45', '0.50', '0.55', '0.60']:
calc = py4vasp.Calculation.from_path("./e03_MgO-gap-optimization/HSE06/aexx" + aexx)
print("AEXX: ", aexx, ", Band gap (eV): ", gap(calc))
```

## Click to see the answer!

AEXX should optimally be 0.55.

Next, use the optimal value of AEXX and vary the screening HFSCREEN!

**loop_hfscreen.sh**

```
for hfscreen in 0.2 0.4 0.6
do
vasp_rm
cp ../WAVECAR .
cat > INCAR << EOF
SYSTEM = MgO screening
ENCUT = 500
ISMEAR = 0
SIGMA = 0.01
LHFCALC = T
HFSCREEN = $hfscreen
AEXX = 0.55
ALGO = Damped
TIME = 0.4
EOF
mpirun -np 2 vasp_std
cat OUTCAR | grep "Total CPU time used"
mkdir hfscreen$hfscreen
cp vaspout.h5 hfscreen$hfscreen/.
done
```

Is the result sensitive to the amount of screening?

```
import py4vasp
import numpy as np
def gap(calc):
homo = np.amax(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] > 0.5])
lumo = np.amin(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] < 0.5])
return lumo - homo
for hf_screen in [ '0.2', '0.4', '0.6']:
calc = py4vasp.Calculation.from_path("./e03_MgO-gap-optimization/HSE06/hfscreen" + hf_screen)
print("HFSCREEN: ", hf_screen, ", Band gap (eV): ", gap(calc))
```

## Click to see the answer!

Yes, less screening and more Fock exchange result in an increased band gap.

**3.4 Questions**¶

- What does the
**HFSCREEN**tag control? - Is
**PBE0**range-dependent? - What contributions are part of the HSE functional?

By the end of this tutorial, you will be able to:

- compute the band structure using a hybrid functional

**4.1 Task**¶

*Perform an HSE06 calculation for CaS, compute the band structure and the band gap.*

The band structure is computed along so-called high-symmetry lines that connect high-symmetry points in the first Brillouin zone. These are of special interest, e.g., due to symmetry imposed degeneracy. In many cases extrema form along these lines so that the $\vec{k}$ points which constitute the HOMO and LUMO can be identified. If it is the same $\vec{k}$ point, this is called *direct* band gap; otherwise it is an *indirect* gap.

There are a couple of tools that help you pick a $\vec{k}$ path. One of them is SeeK-path, which accepts POSCAR files as an input and is accessible via a web interface. You can download this example's POSCAR and try generating your own KPOINTS file!

Figure 1: First Brillouin zone of CaS generated by SeeK-path.

**4.2 Input**¶

The input files to run this example are prepared at `$TUTORIALS/hybrids-part1/e04_CaS-band`

. Check them out!

```
CaS Fm-3m (No.225)
1.0
3.500468 0.000000 2.020996
1.166823 3.300273 2.020996
0.000000 0.000000 4.041992
Ca S
1 1
direct
0.000000 0.000000 0.000000 Ca
0.500000 0.500000 0.500000 S
```

```
SYSTEM = CaS
ENCUT = 500
ISMEAR = -5
LORBIT = 11
```

HSE06-kpoints_opt/INCAR

```
SYSTEM = CaS HSE06
EDIFF = 1e-8
ENCUT = 500
ISMEAR = 0
SIGMA = 0.1
LHFCALC = T
HFSCREEN = 0.3
LORBIT = 11
KPOINTS_OPT_NKBATCH = 15
```

HSE06-explicit/INCAR

```
SYSTEM = CaS HSE06
EDIFF = 1e-8
ENCUT = 500
ISMEAR = 0
SIGMA = 0.1
LHFCALC = T
HFSCREEN = 0.3
LORBIT = 11
```

KPOINTS, HSE06-kpoints_opt/KPOINTS

```
Regular k-points mesh
0
Gamma-centered mesh
4 4 4
0 0 0
```

KPOINTS_OPT, HSE06-kpoints_opt/KPOINTS_OPT

```
Special k-points for band structure
10 ! intersections
line-mode
reciprocal
0.5000000000 0.2500000000 0.7500000000 W
0.5000000000 0.0000000000 0.5000000000 X
0.5000000000 0.0000000000 0.5000000000 X
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.5000000000 0.5000000000 0.5000000000 L
```

HSE06-explicit/KPOINTS

## Click to see the file!

```
Explicit k-points list for band structure
38 !Total number of k-points = 30 (band structure) + 8 (IBZKPT file)
reciprocal
0.00000000000000 0.00000000000000 0.00000000000000 1
0.25000000000000 -0.00000000000000 0.00000000000000 8
0.50000000000000 -0.00000000000000 0.00000000000000 4
0.25000000000000 0.25000000000000 0.00000000000000 6
0.50000000000000 0.25000000000000 0.00000000000000 24
-0.25000000000000 0.25000000000000 0.00000000000000 12
0.50000000000000 0.50000000000000 0.00000000000000 3
-0.25000000000000 0.50000000000000 0.25000000000000 6
0.500000 0.250000 0.750000 0
0.500000 0.222222 0.722222 0
0.500000 0.194444 0.694444 0
0.500000 0.166667 0.666667 0
0.500000 0.138889 0.638889 0
0.500000 0.111111 0.611111 0
0.500000 0.083333 0.583333 0
0.500000 0.055556 0.555556 0
0.500000 0.027778 0.527778 0
0.500000 0.000000 0.500000 0
0.500000 0.000000 0.500000 0
0.444444 0.000000 0.444444 0
0.388889 0.000000 0.388889 0
0.333333 0.000000 0.333333 0
0.277778 0.000000 0.277778 0
0.222222 0.000000 0.222222 0
0.166667 0.000000 0.166667 0
0.111111 0.000000 0.111111 0
0.055556 0.000000 0.055556 0
0.000000 0.000000 0.000000 0
0.000000 0.000000 0.000000 0
0.055556 0.055556 0.055556 0
0.111111 0.111111 0.111111 0
0.166667 0.166667 0.166667 0
0.222222 0.222222 0.222222 0
0.277778 0.277778 0.277778 0
0.333333 0.333333 0.333333 0
0.388889 0.388889 0.388889 0
0.444444 0.444444 0.444444 0
0.500000 0.500000 0.500000 0
```

*Pseudopotential of Ca, and S.*

The POSCAR file contains the face-centered-cubic structure of CaS.

The INCAR files should contain only familiar tags. Otherwise, search on the VASP Wiki!

The self-consistent–field (SCF) calculations of both, PBE and HSE06, are performed on a uniform $\vec{k}$ mesh, that is defined in the KPOINTS file. For both band structures, you can provide the high-symmetry path in *line mode* in a KPOINTS_OPT file. Line mode is an easy way to generate $\vec{k}$ points along a path and the format is described on the VASP Wiki for the KPOINTS file.

There is an important difference between band-structure calculations for DFT and hybrid functionals. Namely, for the band-structure calculation using a hybrid functional, you need both, the $\vec{k}$ points along a high-symmetry path and a uniform $\vec{k}$ mesh to evaluate the Fock-exchange energy. Read the VASP Wiki for detailed advice on how to perform a hybrid band-structure calculation.

For the hybrid band structure, you have two options to compuete the band structure: (i) Using the exlicit list of $\vec{k}$ points as given in HSE06-explicit, and (ii) providing the KPOINTS file with a $\vec{k}$ mesh and the KPOINTS_OPT file with a $\vec{k}$-points path.

The HSE06-explicit/KPOINTS file shows this explicit $\vec{k}$-points list. It comprises 8 irreducible $\vec{k}$ points with their corresponding weights to span the first Brillouin zone, and 30 $\vec{k}$ points with zero weight to define the same path as used for the PBE bands. To construct such an explicit $\vec{k}$-points list, you can copy the content of the IBZKPT files of the DFT SCF and band-structure calculations into one KPOINTS file.

**4.3 Calculation**¶

Open a terminal, navigate to this example's directory and run VASP by entering the following:

```
cd $TUTORIALS/hybrids-part1/e04_*
vasp_std
```

You can plot the band structure using py4vasp. Use the documentation to learn how to plot the character of a band and specify the KPOINTS_OPT file as the source!

```
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e04_CaS-band/band")
my_calc.band.plot?
```

Navigate to the `HSE06`

directory and run VASP by entering the following:

```
cd $TUTORIALS/hybrids-part1/e04_*/HSE06-explicit
cp ../WAVECAR .
mpirun -np 4 vasp_std
```

While this is running, compare the HSE06-explicit/KPOINTS file with the IBZKPT file and the kpoints_opt coordinates in the vasprun.xml file of the PBE calculation! What $\vec{k}$ points are specified? How are the weights determined.

## Click to see the answer!

The $\vec{k}$ points associated with the band structure are the same in both cases. For the hybrid functional, we additionally need to specify a uniform $\vec{k}$ mesh on which the Fock energy is evaluated also for the band-structure calculation. That is why the weight associated with the uniform $\vec{k}$ mesh are chosen in accordance with the irreducible Brillouin zone, and the weight of the band-structure $\vec{k}$ points is zero.

More conviniently, you can compute the HSE band structure by using a KPOINTS_OPT file by entering the following:

```
cd $TUTORIALS/hybrids-part1/e04_*/HSE06-kpoints_opt
cp ../HSE06-explicit/WAVECAR .
mpirun -np 4 vasp_std
```

Check the band structure by plotting it!

After the calculation finished, the HSE band structure is available. You can plot it again by specifying the KPOINTS_OPT file as the source.

```
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e04_CaS-band/HSE06-kpoints_opt")
my_calc.band.plot("1, 2", source="kpoints_opt")
```

For the explicit list of $\vec{k}$ points, you need to zoom in on the zero-weighted $\vec{k}$ points towards the end. Try to find the correct region by comparison with the band-structure plot above!

```
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e04_CaS-band/HSE06-explicit")
my_calc.band.plot("1, 2")
```

Finally, compute the band gap and compare the value to the band structure. Is the band gap direct or indirect? At which $\vec{k}$ point? How does it compare to the experimental value of 5.4 eV?

```
import py4vasp
import numpy as np
def gap(calc):
homo = np.amax(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] > 0.5])
lumo = np.amin(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] < 0.5])
return lumo - homo
pbe_calc = py4vasp.Calculation.from_path("./e04_CaS-band")
hse06_calc = py4vasp.Calculation.from_path("./e04_CaS-band/HSE06-explicit")
print("PBE band gap (eV): ", gap(pbe_calc))
print("HSE06 band gap (eV): ", gap(hse06_calc))
```

## Click to see the answer!

The band gap is indirect, since the LUMO is at the X point and the HOMO is at the 𝝘 point.

**4.4 Questions**¶

- What steps are necessary to compute the band structure of systems using hybrid functionals?
- Why do the $\vec{k}$ points associated with the high-symmetry points have zero weight in this example?
- What is an indirect band gap?