Part 1: An overview of available functionals


1 Band gap of Si with the Perdew-Burke-Ernzerhof (PBE) and PBE0 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.

The computationally most efficient functionals are the local density approximation (LDA), the generalized gradient approximations (GGA) and the meta-GGA functionals. The GGA and meta-GGA functionals are of the semilocal type and a plethora of them have been proposed. The GGA PBE functional, implemented in VASP, is the most commonly used in solid-state physics. However, note that the LDA and GGA functionals do not provide a correct description of the band gap. Whether the Kohn-Sham (KS) orbitals $\psi_{n\mathbf{k}}(\mathbf{r})$ reflect the right physics crucially depends on the type of the used exchange-correlation (xc) functional. In this respect meta-GGA and hybrid functionals are more appropriate for band gap calculation.

Hybrid functionals make up a special category of xc functionals that mix some fraction of Fock exchange into a semilocal functional. This approach is inspired by the well-known Hartree–Fock (HF) method, which is the exact solution for uncorrelated systems. Regarding extended systems, it aims to cure the infamous underestimation of band gaps by the LDA and GGA. The nonlocal Fock potential reads $$\tag{1} V_{\mathrm{x}}^{\mathrm{HF}}\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 at a specific $\mathbf{k}$ point, $\mathbf{r}$ is the position and $e$ the charge of an electron. Equation (1) depends on the KS orbitals instead of the density 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 Fock exchange energy (the corresponding potential is Equation (1)), while $E_{\mathrm{x}}^{\mathrm{PBE}}$ and $E_{\mathrm{c}}^{\mathrm{PBE}}$ are the PBE exchange and correlation energies, respectively. Compared to many other functionals, PBE0 contains no parameter that were fitted to experimental data.

1.2 Input

The input files to run this example are prepared at $TUTORIALS/hybrids/e01_Si-gap. Check them out!

POSCAR


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

INCAR


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

KPOINTS


K-Points
 0
Gamma
 7 7 7
 0 0 0

POTCAR


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, AEXX, AGGAC and ALDAC, have their default values. How do these tags relate to Equation (2)?

Click to see the answer!

The POTCAR file already defines PBE as the density functional. Therefore, it is not necessary to explicitly set the GGA tag. AEXX, AGGAC and ALDAC correspond to the prefactors in Equation (2).

In the PBE calculation, ISMEAR = -5 switches on the tetrahedron method with Blöchl corrections to determine the partial occupations $f_{n\mathbf{k}}$. This smearing yields the smoothest density of states. With LORBIT = 11, VASP writes the projected DOS 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. This direct optimization algorithm is generally more robust than iterative methods. In many systems, hybrid functionals require this robustness to achieve convergence. Why do we not immediately compute the DOS, like in the PBE calculation?

Click to see the answer!

The Blöchl corrections cause energies not to be variational, thus we use Gaussian smearing with ISMEAR = 0 and SIGMA = 0.01 instead. With the Gaussian smearing, the energies are variational as required for the damped algorithm. After convergence, we perform a single-shot calculation without updating the orbitals (ALGO=None) to get the DOS with the tetrahedron method (ISMEAR = -5) and orbital projections (LORBIT = 11).

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

The overall settings are such that the precision for the band gap lies within 0.1 eV. What would you need to do in order 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 $\mathbf{k}$ points 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/e01_*
vasp_std

Plot the DOS using py4vasp for the PBE calculation.

In [1]:
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 valence band maximum and the conduction band minimum. In the DOS plot, this corresponds to the energy range without any states in the vicinity of the Fermi energy. Note that the energy axis is aligned to the Fermi energy.

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

cd $TUTORIALS/hybrids/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.

Can you compare the total energy of PBE and PBE0?

Click to see the answer!

No, generally it does not make sense to compare total energies obtained with different exchange-correlation functionals. Only physical quantities like the band gap can be compared.

Let us now compare the band gap predicted by PBE and PBE0!

In [2]:
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")

for calc in (pbe_calc, pbe0_calc):
    print(f"{str(calc.system):10} {gap(calc):6.3f}")
cd Si PBE   0.616
cd Si PBE0  1.842

Which band gap is larger?

Click to see the answer!

The PBE0 band gap (1.84 eV) is larger than the PBE one (0.62 eV). As a general trend in many materials, hybrid functionals lead to larger band gap and partially correct for the underestimation obtained with GGA functionals compared to experiment.

To visualize the PBE0 gap, you can again compute the DOS. Open the INCAR file of your PBE0 calculation and set ISMEAR = -5, ALGO = None, and LORBIT = 11. You may also remove the tag SIGMA and TIME which are not used. Then, restart the calculation with the existing WAVECAR file by entering the following:

cd $TUTORIALS/hybrids/e01_*/PBE0
mpirun -np 2 vasp_std
In [3]:
import py4vasp
pbe_calc = py4vasp.Calculation.from_path("./e01_Si-gap")
pbe0_calc = py4vasp.Calculation.from_path("./e01_Si-gap/PBE0")

pbe_calc.dos.plot().label("PBE") + pbe0_calc.dos.plot().label("PBE0")

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

Click to see the answer!

The PBE0 DOS is different. Not only is the band gap enlarged, but also the relative energy between bands and the bands width are different. These subtle changes can crucially influence the properties of the material.

1.4 Questions

  1. What is a hybrid functional?
  2. How can you compute the band gap of a system?
  3. What does the LHFCALC tag control?
  4. Why should you not combine ISMEAR = -5 with ALGO = Damped?

2 Band gap of Ar with the PBE and B3LYP functionals and the Hartree-Fock (HF) method

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

  • use the B3LYP functional
  • use the HF method
  • explain the limitations of HF and B3LYP

2.1 Task

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

B3LYP was originally designed to describe molecular properties and is a semiempirical hybrid functional. Semiempirical implies that parameters were fitted to reproduce experimental results for properties like the atomization energies, electron-proton affinities and ionization potentials for atoms and molecules of a test set. The empirical parameters are the ones that multiply the various components of the B3LYP functional:

$$\tag{1 a} E_{\mathrm{xc}}^{\mathrm{B3LYP}} = E_{\mathrm{x}}^{\mathrm{B3LYP}} + E_{\mathrm{c}}^{\mathrm{B3LYP}} $$$$ \tag{1 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{1 c} E_{\mathrm{c}}^{\mathrm{B3LYP}} = 0.19 E_{\mathrm{c}}^{\mathrm{VWN3}} + 0.81 E_{\mathrm{c}}^{\mathrm{LYP}}. $$

In a HF calculation, the exchange consists of the full Fock exchange (AEXX = 1) and no correlation effects are considered ($E_{\mathrm{c}}=0$).

B3LYP and other unscreened functionals inherit a fundamental flaw of the HF method—its failure in the limit of the homogeneous electron gas. This failure originates from a singularity of the Fourier transform of the Fock potential at $\mathbf{q}=0$

$$\tag{2} V_{\mathrm{x}}^{\mathrm{HF}}\left(\mathbf{q}\right)= -\frac{4\pi e^2}{2}\sum_{n\mathbf{k}}f_{n\mathbf{k}} \frac{\psi_{n\mathbf{k}}^{*}(\mathbf{r}')\psi_{n\mathbf{k}}(\mathbf{r})} { \mathbf{q} }~. $$

This causes a complete suppression of itinerant states. In other words, the HF method and hybrid functionals 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 of 14.4 eV in experiment!

2.2 Input

The input files to run this example are prepared at $TUTORIALS/hybrids/e02_Ar-gap. Check them out!

POSCAR


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

INCAR


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

KPOINTS


K-Points
 0
Gamma
4  4  4
0  0  0

POTCAR


Pseudopotential of Ar.


The POSCAR file defines the face-centered-cubic Ar crystal.

The first INCAR file runs 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 HF is set by AEXX. So only Fock exchange and no semilocal functional exchange-correlation is used. Recall the meaning of the remaining tags! What are the default values of GGA, AGGAX, AGGAC and ALDAC here?

Click to see the answer!

GGA = PE, AGGAX = 1 – AEXX = 0. Since AEXX = 1, AGGAC and ALDAC are automatically set to zero. See AEXX on the VASP Wiki!

For the B3LYP INCAR, the semilocal functional is changed using the GGA tag and the appropriate fraction for each exchange and correlation component is set. Compare the fractions to Equation (1)!

2.3 Calculation

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

cd $TUTORIALS/hybrids/e02_*
vasp_std

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

In [4]:
import py4vasp
pbe_calc = py4vasp.Calculation.from_path("./e02_Ar-gap")

pbe_calc.dos.plot()

Next, copy the WAVECAR to the HF directory and run VASP by entering the following:

cd $TUTORIALS/hybrids/e02_*/HF
cp ../WAVECAR .
mpirun -np 2 vasp_std

Finally, do the same for B3LYP and then compare the band gaps! Which method yields the band gap closest to the experimental value?

In [5]:
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(f"{str(calc.system):10} {gap(calc):6.3f}")
Ar PBE      8.504
Ar HF      17.612
Ar B3LYP   10.313
Click to see the answer!

While the PBE band gap of 8.5 eV 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

  1. Which class of systems do the PBE0 and B3LYP functionals and HF method fail to describe well? Why?
  2. What does the AGGAC tag control?

3 Band gap optimization for MgO for screened hybrid functionals

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

  • understand the fundamentals of screened hybrid functionals
  • explain the basic assumptions taken by the HSE06 functional
  • optimize the fraction of exchange and the screening length for the band gap

3.1 Task

Optimize the fraction of Fock exchange in a screened hybrid functional in order to reproduce the experimental band gap of MgO.

Screened hybrid functionals decompose the $1/r$ dependence of the Fock exchange into a short-range (SR) and a long-range (LR) contribution introducing a screening parameter $\mu$:

$$\tag{1} E_{\mathrm{xc}}^{\mathrm{screened}} = \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,LR}}(\mu) + E_{\mathrm{c}}^{\mathrm{PBE}}. $$

The SR limit is identical to the PBE0 functional. But by truncating the LR contribution, screened hybrid functionals avoid the singularity in the Fock exchange. This singularity is responsible for the suppression of all metallic states. Additionally, screened hybrid functionals allow to reduce the computational cost by downsampling.

In solid-state physics, one commonly used functional is HSE06. In this functional, the screening parameter was empirically optimized for a value of μ = 0.2 Å$^{−1}$. The corresponding interaction range $2 / \mu \approx 10$ Å typically reaches a few nearest neighbors. In VASP, you control the screening length by the HFSCREEN tag.

The gap size still crucially depends on the fraction of Fock exchange introduced. Let us take a practical approach here and optimize the mixing to reproduce the experimental band gap of MgO of 7.7 eV! Note that you sacrifice the predictive power of a first-principles approach here. Instead, this approach is closer to model physics where you want to get qualitative insights into the system. Please do not refer to this fitted functional as HSE as it is occasionally done in the literature. HSE06 corresponds only to AEXX = 0.25 and HFSCREEN = 0.2.

3.2 Input

The input files to run this example are prepared at $TUTORIALS/hybrids/e03_MgO-gap-optimization. Check them out!

POSCAR


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

INCAR


SYSTEM = MgO

ENCUT  = 500
ISMEAR = -5

LORBIT = 11

screened_hybrid/INCAR


SYSTEM   = MgO HSE06

ENCUT    = 500
ISMEAR   = 0
SIGMA    = 0.01

LHFCALC  = T
HFSCREEN = 0.2
AEXX     = 0.25

ALGO     = Damped
TIME     = 0.4

KPOINTS


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

POTCAR


Pseudopotential of Mg, and O.

Check the meaning of all tags used in the INCAR file and make comments in the file! Other related tags are HFRCUT and HFALPHA.

3.3 Calculation

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

cd $TUTORIALS/hybrids/e03_*
vasp_std

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

In [6]:
import py4vasp
pbe_calc = py4vasp.Calculation.from_path("./e03_MgO-gap-optimization")

pbe_calc.dos.plot()

Then, in the screened_hybrid directory run a loop over different values of the fraction of HF exchange AEXX by entering the following:

cd $TUTORIALS/hybrids/e03_*/screened_hybrid
bash loop_aexx.sh

loop_aexx.sh


for aexx in 0.40 0.45 0.50 0.55
do
vasp_rm
cp ../WAVECAR .
cat > INCAR << EOF
SYSTEM   = MgO $aexx

ENCUT    = 500
ISMEAR   = 0
SIGMA    = 0.01

LHFCALC  = T
HFSCREEN = 0.2
AEXX     = $aexx

ALGO     = Damped
TIME     = 0.4
EOF

mpirun -np 2 vasp_std

mkdir aexx$aexx
cp vaspout.h5 aexx$aexx/.
done

Compute the band gap for each value of AEXX to extract the optimal value of AEXX that reproduces the experimental band gap of 7.7eV for this material!

In [7]:
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.40', '0.45', '0.50', '0.55']:
    calc = py4vasp.Calculation.from_path(f"./e03_MgO-gap-optimization/screened_hybrid/aexx{aexx}")
    print(f"AEXX:  {aexx},  Band gap (eV): {gap(calc):6.3f}")
AEXX:  0.40,  Band gap (eV):  7.219
AEXX:  0.45,  Band gap (eV):  7.580
AEXX:  0.50,  Band gap (eV):  7.946
AEXX:  0.55,  Band gap (eV):  8.313
Click to see the answer!

The optimal AEXX should be between 0.45 and 0.50.

Next, we use a value of AEXX close to the optimum and vary the screening HFSCREEN!

bash loop_hfscreen.sh

loop_hfscreen.sh


for hfscreen in 0.15 0.20 0.25 0.30
do
vasp_rm
cp ../WAVECAR .
cat > INCAR << EOF
SYSTEM   = MgO $hfscreen

ENCUT    = 500
ISMEAR   = 0
SIGMA    = 0.01

LHFCALC  = T
HFSCREEN = $hfscreen
AEXX     = 0.5

ALGO     = Damped
TIME     = 0.4
EOF

mpirun -np 2 vasp_std

mkdir hfscreen$hfscreen
cp vaspout.h5 hfscreen$hfscreen/.
done

Is the result sensitive to the fraction of screening?

In [8]:
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 hfscreen in ['0.15', '0.20', '0.25', '0.30']:
    calc = py4vasp.Calculation.from_path(f"./e03_MgO-gap-optimization/screened_hybrid/hfscreen{hfscreen}")
    print(f"HFSCREEN: {hfscreen}, Band gap (eV): {gap(calc):6.3f}")
HFSCREEN: 0.15, Band gap (eV):  8.291
HFSCREEN: 0.20, Band gap (eV):  7.946
HFSCREEN: 0.25, Band gap (eV):  7.643
HFSCREEN: 0.30, Band gap (eV):  7.372
Click to see the answer!

Yes, a shorter screening length increases the amount of Fock exchange and consequently the band gap.

3.4 Questions

  1. What does the HFSCREEN tag control?
  2. Is PBE0 range dependent?
  3. What contributions are part of the HSE06 functional?

4 Band structure of CaS with the PBE and HSE06 functionals

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 $\mathbf{k}$ points which constitute the valence band maximum and conduction band minimum (also called HOMO and LUMO) can be identified. If they are at the same $\mathbf{k}$ point, the band gap is called direct, otherwise it is indirect.

There are a couple of tools that help you pick a $\mathbf{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/e04_CaS-band. Check them out!

POSCAR


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

INCAR


SYSTEM = CaS

ENCUT  = 500
ISMEAR = -5

LORBIT = 11

band/INCAR


SYSTEM = CaS band

ENCUT  = 500
ISMEAR = 0
ISIGMA = 0.01

ICHARG = 11
LORBIT = 11

HSE06/INCAR


SYSTEM   = CaS HSE06

ENCUT    = 500
ISMEAR   = 0
SIGMA    = 0.01

LHFCALC  = T
HFSCREEN = 0.2

ALGO     = Damped
TIME     = 0.4

HSE06/band/INCAR


SYSTEM   = CaS HSE06 band

ENCUT    = 500
ISMEAR   = 0
SIGMA    = 0.01
NELMIN   = 4

LHFCALC  = T
HFSCREEN = 0.2
LORBIT   = 11


KPOINTS


K-Points
 0
Gamma
 6 6 6
 0 0 0

band/KPOINTS


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

    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/KPOINTS


K-Points
 0
Gamma
 6 6 6
 0 0 0

HSE06/band/KPOINTS


Click to see the file!
Explicit k-points list for band structure
53 ! Total number of k-points = 37 (band structure) + 16 (IBZKPT file)
reciprocal
    0.00000000000000    0.00000000000000    0.00000000000000             1
    0.16666666666667   -0.00000000000000    0.00000000000000             8
    0.33333333333333   -0.00000000000000    0.00000000000000             8
    0.50000000000000   -0.00000000000000    0.00000000000000             4
    0.16666666666667    0.16666666666667    0.00000000000000             6
    0.33333333333333    0.16666666666667    0.00000000000000            24
    0.50000000000000    0.16666666666667    0.00000000000000            24
   -0.33333333333333    0.16666666666667    0.00000000000000            24
   -0.16666666666667    0.16666666666667    0.00000000000000            12
    0.33333333333333    0.33333333333333    0.00000000000000             6
    0.50000000000000    0.33333333333333    0.00000000000000            24
   -0.33333333333333    0.33333333333333    0.00000000000000            12
    0.50000000000000    0.50000000000000    0.00000000000000             3
    0.50000000000000    0.33333333333333    0.16666666666667            24
   -0.33333333333333    0.33333333333333    0.16666666666667            24
   -0.33333333333333    0.50000000000000    0.16666666666667            12
  0.500000  0.500000  0.500000      0  L
  0.500000  0.472222  0.527778      0
  0.500000  0.444444  0.555556      0
  0.500000  0.416667  0.583333      0
  0.500000  0.388889  0.611111      0
  0.500000  0.361111  0.638889      0
  0.500000  0.333333  0.666667      0
  0.500000  0.305556  0.694444      0
  0.500000  0.277778  0.722222      0
  0.500000  0.250000  0.750000      0  W
  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  X
  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  Gamma
  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  L

POTCAR


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, please look up their meaning on the VASP Wiki!

The self-consistent calculations of both, PBE and HSE06, are performed on a uniform $\mathbf{k}$ mesh, that is defined in the KPOINTS file. Then, the band/KPOINTS file for the PBE band-structure is written using the line mode. This is an easy way to generate $\mathbf{k}$ points along a path.

For the band-structure calculation using a hybrid functional, you need both, the $\mathbf{k}$ points along a high-symmetry path and a uniform $\mathbf{k}$ mesh to evaluate the Fock exchange energy. This is why the HSE06/band/KPOINTS file comprises 16 irreducible $\mathbf{k}$ points with their corresponding weights to span the first Brillouin zone, and 37 $\mathbf{k}$ points with zero weight to define the same path as used for the PBE bands.

4.3 Calculation

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

cd $TUTORIALS/hybrids/e04_*
vasp_std

Next, compute the band structure of the PBE calculation by entering the following:

cd $TUTORIALS/hybrids/e04_*/band
cp ../CHGCAR .
vasp_std

You can plot the band structure using py4vasp. Use the documentation to learn how to plot the character of a band! Then, replace the character ? of my_calc_band.plot by what is appropriate.

In [9]:
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e04_CaS-band/band")

my_calc.band.plot?
Signature: my_calc.band.plot(*args, **kwargs)
Docstring:
Wrapper around the :py:meth:`to_graph` function.

This will merge multiple graphs if you specify different sources with the
selection arguments. All arguments are passed to the wrapped function.
File:      /opt/tljh/user/lib/python3.9/site-packages/py4vasp/_third_party/graph/mixin.py
Type:      method

Navigate to the HSE06 directory and run VASP by entering the following:

cd $TUTORIALS/hybrids/e04_*/HSE06
cp ../WAVECAR .
mpirun -np 2 vasp_std

While this is running, compare the HSE06/band/KPOINTS file with the IBZKPT file and the reciprocal coordinates in the band/OUTCAR file of the PBE calculation! What $\mathbf{k}$ points are specified? How are the weights determined.

Click to see the answer!

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

Next, compute the band structure of the HSE06 calculation by entering the following:

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

This calculation will take some time. In the meantime, do you know why we set NELMIN in the INCAR file?

Click to see the answer!

Since we restart the calculation from an existing WAVECAR file, the total energy of the system is already converged. Therefore VASP will stop updating the orbitals after NELMIN iterations. To obtain the eigenvalues, the Kohn-Sham Hamiltonian is diagonalized in the basis of these orbitals. The orbitals at the additional $\mathbf{k}$ points may not be sufficiently converged yet, so that the basis is not good enough to get accurate eigenvalues. If that happens you will see jumps in the band structure. To avoid this, we ensure that at least four updates are done.

Let us look at the band structure with the hybrid functional now. You will see all the $\mathbf{k}$ points, also the one from the regular grid. Please use the interactive features to zoom to the area between the two L-point labels.

In [10]:
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e04_CaS-band/HSE06/band")

my_calc.band.plot()

If you look closely at the conduction band, slight kinks can still be seen. This can be improved by running a longer run (NELMIN) or using more bands (NBANDS).

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

In [11]:
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/band")
hse06_calc = py4vasp.Calculation.from_path("./e04_CaS-band/HSE06/band")

for calc in (pbe_calc, hse06_calc):
    print(f"{str(calc.system):15} {gap(calc):6.3f}")
CaS band         2.382
CaS HSE06 band   3.457
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

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

Good job! You have finished Part 1 of the tutorials about hybrid functionals!

Go to Top $\uparrow$