Part 1: An overview of available functionals


1 Band gap of Si with Perdew-Burke-Ernzerhof (PBE), and PBE0

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!

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

Click to see the answer!

The POTCAR file already defines PBE as the density functional. So it is not necessary to explicitly set the GGA tag. AEXX, AGGAC and ALDAC correspond to the prefactors in 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.

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

Can you compare the total energy of PBE and PBE0?

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.

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

In [3]:
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))
Si bandgap
PBE: 0.6162992491537098
PBE0: 1.842131901640176

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
In [4]:
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.

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 Ar band gap with PBE, B3LYP and Hartree-Fock (HF) method

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!

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, 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?

Click to see the answer!

GGA=PE, AGGAX=1-AEXX=0. And due to AEXX=1 also AGGAC and ALDAC are automatically set to zero. See AEXX on the VASP Wiki!

The B3LYP INCAR file also switches the evaluation of the Fock exchange on, but additionally the density functional is changed using the GGA tag and a special fraction of each exchange and correlation energy is set. Compare the fractions to Equation (2)!

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!

In [5]:
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-part1/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 [6]:
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))
e02_Ar-gap 8.503551608603647
HF 17.61200373040635
B3LYP 10.312770386932943
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

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

3 Band-gap optimization for MgO with Heyd-Scuseria-Ernzerhof (HSE) functional

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!

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

HSE06/INCAR


SYSTEM   = MgO HSE06

ENCUT    = 500
ISMEAR   = 0
SIGMA    = 0.01

LHFCALC  = T
HFSCREEN = 0.2
AEXX     = 

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-part1/e03_*
vasp_std

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

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

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 [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 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))
AEXX:  0.45 , Band gap (eV):  7.067342854440238
AEXX:  0.50 , Band gap (eV):  7.371794000924975
AEXX:  0.55 , Band gap (eV):  7.679841295771869
AEXX:  0.60 , Band gap (eV):  7.989659257522117
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?

In [9]:
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))
HFSCREEN:  0.2 , Band gap (eV):  8.313243887149676
HFSCREEN:  0.4 , Band gap (eV):  7.155342084436393
HFSCREEN:  0.6 , Band gap (eV):  6.326655843105992
Click to see the answer!

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

3.4 Questions

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

4 Band structure of CaS using PBE and HSE06

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!

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

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

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, 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!

In [10]:
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.

In [11]:
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!

In [12]:
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?

In [13]:
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))
PBE   band gap (eV):  2.3718690482822926
HSE06 band gap (eV):  3.228142587546561
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 $\vec{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$