Part 2: Optical absorption of LiF

Content

 4 Preparatory ground-state calculation $G_0W_0$
 5 Optical absorption in the independent particle approximation (BSE-IP)
 6 Optical absorption in BSE-RPA
 7 Optical absorption in BSE-TDA
 8 Optical absorption in BSE beyond TDA

4 Preparatory ground-state calculation $G_0W_0$

At the end of this tutorial you will be able to:

  • Set up a one-shot $G_0W_0$ calculation
  • Find the QP band gap

4.1 Task

Compute the $G_0W_0$ band structure and screened Coulomb potential for LiF unit cell based on PBE electronic structure.

In this part of the tutorial, we study LiF. Unlike diamond, LiF has weakly screened interaction between electrons, which leads to a strongly bound Frenkel exciton.

As in the case of diamond, we will use the $G_0W_0$ approximation as a starting point for the BSE calculation. The $GW$ calculation can be performed in three steps:

  1. DFT ground state
  2. DFT with exact diagonalization to compute unoccupied bands
  3. $G_0W_0$ calculation

4.2 Input

The input files to run this example are prepared at $TUTORIALS/BSE/e04_LiF_GW.

POSCAR

LiF
4.026
  0 0.5 0.5
  0.5 0 0.5
  0.5 0.5 0
Li F
1  1
Direct
  0.00 0.00 0.00 Li
  0.50 0.50 0.50 F

ground-state/INCAR

SYSTEM  = LiF DFT

ISMEAR  = 0      ! Gaussian smearing
SIGMA   = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT   = 500    ! energy cutoff

unoccupied-states/INCAR

SYSTEM  = LiF DFT unoccupied states

ISMEAR  = 0      ! Gaussian smearing
SIGMA   = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT   = 500    ! energy cutoff

NBANDS  = 64     ! 

ALGO    = Exact  ! use exact diagonalization of the Hamiltonian
NELM    = 1      ! since we are already converged stop after one step

LOPTICS = T      ! write WAVEDER

INCAR

SYSTEM    = LiF G0W0

ISMEAR    = 0      ! Gaussian smearing
SIGMA     = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT     = 500    ! energy cutoff

NBANDS    = 64     ! consistent with the ground state

! G0W0
ALGO      = EVGW0  ! use "GW0" for VASP.5.X
NELM      = 1
KPAR      = 4

NOMEGA    = 50     ! default is 100

KPOINTS

k-mesh
 0
Gamma
 6 6 6
 0 0 0

POTCAR

GW pseudopotentials of Li and F

4.3 Calculation

These steps are the same as in Part 1 of the tutorial, so please go to Part 1 for a detailed description.

4.3.1 DFT ground state

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

cd $TUTORIALS/BSE/e04_LiF_GW/ground-state
mpirun -np 4 vasp_std
4.3.2 DFT unoccupied bands

Navigate to this example's directory and run VASP by entering the following:

cd $TUTORIALS/BSE/e04_LiF_GW/unoccupied-states
cp ../ground-state/WAVECAR .
mpirun -np 4 vasp_std
4.3.3 GW calculation

Now we can perform $G_0W_0$ using the WAVECAR and WAVEDER from the DFT unoccupied-bands step.

cd $TUTORIALS/BSE/e04_LiF_GW
cp unoccupied-states/{WAVECAR,WAVEDER} .
mpirun -np 4 vasp_std

Check in the OUTCAR file that the calculated direct band gap for LiF is ~13.44 eV.

4.4 Questions

  1. Is LiF a direct of indirect band gap system?
  2. What is the RPA dielectric constant of LiF?

5 Optical absorption in the independent particle approximation (BSE-IP)

At the end of this tutorial you will be able to:

  • Compute optical absorption in the independent particle approximation

5.1 Task

Compute the IP spectrum of LiF based on $G_0W_0$ quasiparticles.

In the independent particle approximation LHARTREE = .FALSE. and LADDER = .FALSE.. Hence, the BSE Hamiltonian only consists of the diagonal term ($A=D$), i.e., the interaction between the particles is neglected. For many systems where the screening is strong and the electron-hole binding energy is small, such an approximation provides a reasonable description of the optical absorption.

5.2 Input

The input files to run this example are prepared at $TUTORIALS/BSE/e05_LiF_IP.

The direct and the exchange terms of the BSE Hamiltonian are not taken into account, i.e., LHARTREE=.FALSE. and LADDER=.FALSE.

INCAR

SYSTEM    = LiF IP

ISMEAR    = 0      ! Gaussian smearing
SIGMA     = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT     = 500    ! energy cutoff

NBANDS    = 64     ! consistent with the ground state

! BSE
ALGO      = BSE    ! use BSE
ANTIRES   = 0      ! Tamm-Dancoff approximation
NBANDSV   = 5      ! number of conduction bands in BSE
NBANDSO   = 5      ! number of valence bands in BSE
CSHIFT    = 0.2    ! complex shift
OMEGAMAX  = 40     ! max. energy
LHARTREE = .FALSE. ! no RPA diagrams
LADDER   = .FALSE. ! no ladder diagrams

5.3 Calculation

Navigate to this example's directory and run VASP by entering the following:

cd $TUTORIALS/BSE/e05_LiF_IP
cp $TUTORIALS/BSE/e04_LiF_GW/{WAVECAR,WAVEDER,W*.tmp} .
mpirun -np 4 vasp_std

Plot the calculated IP dielectric function with py4vasp

In [5]:
import py4vasp
import numpy as np
import os

E, ε = np.loadtxt(os.path.expanduser("./e09_LiF_TDHF_shifted_grid/expt.dat")).T

ipa = py4vasp.Calculation.from_path("./e05_LiF_IP")
(
    py4vasp.plot(E, ε, "Expt") +
    ipa.dielectric_function.plot("BSE(Im)").label("BSE-IPA")
)

The agreement with the experimental spectrum is not good. Most importantly, the strong excitonic peak at the absorption edge is not present. The absorption edge is at 13.4 eV, which means that electrons and holes do not form a bound state which would reduce the transition energy.

This spectrum should be equivalent to the one obtained if LOPTICS = .TRUE. with the same orbitals and energies.

5.4 Questions

  1. Are excitonic effects accounted for in the IP approximation?
  2. For what systems can the IP be a good approximation?

6 Optical absorption in BSE-RPA

At the end of this tutorial you will be able to:

  • Compute the RPA spectrum of LiF based on $G_0W_0$ quasiparticles

6.1 Task

Compute the RPA spectrum of LiF based on $G_0W_0$ quasiparticles.

In the RPA approximation, the ladder diagrams are not calculated, i.e., LADDER = .FALSE., and only the exchange interaction is included LHARTREE = .TRUE. in the Hamiltonian $A = D + K^x$.

6.2 Input

The input files to run this example are prepared at $TUTORIALS/BSE/e06_LiF_RPA.

INCAR

SYSTEM    = LiF RPA

ISMEAR    = 0      ! Gaussian smearing
SIGMA     = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT     = 500    ! energy cutoff

NBANDS    = 64     ! consistent with the ground state

! BSE
ALGO      = BSE    ! use BSE
ANTIRES   = 0      ! Tamm-Dancoff approximation
NBANDSV   = 5      ! number of conduction bands in BSE
NBANDSO   = 5      ! number of valence bands in BSE
CSHIFT    = 0.2    ! complex shift
OMEGAMAX  = 40     ! max. energy
LHARTREE = .TRUE.  ! RPA diagrams
LADDER   = .FALSE. ! no ladder diagrams

6.3 Calculation

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

cd $TUTORIALS/BSE/e06_LiF_RPA
cp $TUTORIALS/BSE/e04_LiF_GW/{WAVECAR,WAVEDER,W*.tmp} .
mpirun -np 4 vasp_std

Plot the calculated RPA dielectric function with py4vasp

In [7]:
import py4vasp
import numpy as np
import os

E, ε = np.loadtxt(os.path.expanduser("./e09_LiF_TDHF_shifted_grid/expt.dat")).T

ipa = py4vasp.Calculation.from_path("./e05_LiF_IP")
rpa = py4vasp.Calculation.from_path("./e06_LiF_RPA")
(
    py4vasp.plot(E, ε, "Expt") +
    ipa.dielectric_function.plot("BSE(Im)").label("BSE-IPA") +
    rpa.dielectric_function.plot("BSE(Im)").label("BSE-RPA")
)

The agreement with the experiment has not improved and that the spectrum is slightly blueshifted. This shift is due to the repulsive exchange interaction.

Here, the RPA approximation only includes the bubble diagrams in the Tamm-Dancoff approximation, i.e., the terms B and B* in Eq. (1) are neglected and the calculated spectrum is not equivalent to ALGO = CHI due to the lack of the coupling between excitation and de-excitation "bubbles".

The RPA approximation has similar shortcoming as the IP, i.e., a lack of excitons. However, it accounts for the interaction with the plasmon excitations and often provides a good description of the EELS spectrum.

6.4 Questions

  1. Are excitonic effects accounted for in the RPA approximation?
  2. For what systems can the RPA be a good approximation?

7 Optical absorption in BSE-TDA

At the end of this tutorial you will be able to:

  • Perform the optical absorption in the Tamm-Dancoff approximation

7.1 Task

Compute the optical absorption of LiF in the Tamm-Dancoff approximation based on the $G_0W_0$ quasiparticles.

LiF is a wide-gap semiconductor with a small dielectric constant $\varepsilon_\infty=1.9$ and hence weak screening. Due to this weak screening, the attractive interaction between electrons and holes is large, which leads to the strong localization of excitons and large exciton binding energy 0.1-1 eV. Strong electron-hole interaction is typical for alkali halides and organic molecular systems.

In the Tamm-Dancoff approximation, the bubble and the ladder diagrams are taken into account, i.e., LHARTREE = .TRUE. and LADDER = .TRUE. However, the matrices $B$ and $B^*$ are neglected leading to the Eq. (8), where the Hamiltonian consists of three terms: $A=D+K^x+K^D$.

7.2 Input

The input files to run this example are prepared at $TUTORIALS/BSE/e07_LiF_BSE.

INCAR

SYSTEM    = LiF BSE

ISMEAR    = 0      ! Gaussian smearing
SIGMA     = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT     = 500    ! energy cutoff

NBANDS    = 64     ! consistent with the ground state

! BSE
ALGO      = BSE    ! use BSE
ANTIRES   = 0      ! Tamm-Dancoff approximation
NBANDSV   = 5      ! number of conduction bands in BSE
NBANDSO   = 5      ! number of valence bands in BSE
CSHIFT    = 0.2    ! complex shift
OMEGAMAX  = 40     ! max. energy
LHARTREE  = .TRUE. ! bubble diagrams
LADDER    = .TRUE. ! ladder diagrams

7.3 Calculation

Navigate to this example's directory and run VASP by entering the following:

cd $TUTORIALS/BSE/e07_LiF_BSE
cp $TUTORIALS/BSE/e04_LiF_GW/{WAVECAR,WAVEDER,W*.tmp} .
mpirun -np 4 vasp_std

Plot the calculated TDA dielectric function with py4vasp

In [9]:
import py4vasp
import numpy as np
import os

E, ε = np.loadtxt(os.path.expanduser("./e09_LiF_TDHF_shifted_grid/expt.dat")).T

ipa = py4vasp.Calculation.from_path("./e05_LiF_IP")
rpa = py4vasp.Calculation.from_path("./e06_LiF_RPA")
tda = py4vasp.Calculation.from_path("./e07_LiF_BSE")
(
    py4vasp.plot(E, ε, "Expt") +
    ipa.dielectric_function.plot("BSE(Im)").label("BSE-IPA") +
    rpa.dielectric_function.plot("BSE(Im)").label("BSE-RPA") +
    tda.dielectric_function.plot("BSE(Im)").label("BSE-TDA")
)

Now the spectrum is in a much better agreement with experiment. The first strong excitonic peak is properly captured. Although, the intensities of some of the features do not match the experiment. We will see in Part 3 that these issues are resolved using a denser k-points mesh.

7.4 Questions

  1. Does the Tamm-Dancoff approximation account for the excitonic effects?
  2. What terms of the Eq. (1) are neglected in the TDA?

8 Optical absorption in BSE beyond TDA

At the end of this tutorial you will be able to:

  • Calculate the optical absorption spectrum beyond TDA

8.1 Task

Compute the optical absorption of LiF beyond the Tamm-Dancoff approximation based on the $G_0W_0$ quasiparticles.

Finally, we are going to solve the full BSE equation, i.e., include the resonant - anti-resonant coupling $$H^{BSE}=\begin{pmatrix} A & B\\ B^* & A^*\\ \end{pmatrix} $$

If ANTIRES = 2, the terms $B$ and $B^*$ are included in the calculation, which in principle means that the matrix of rank $2*N_c*N_v*N_k$ must be diagonalized. However, VASP uses the time inversion-symmetry to reduce the problem by two equations of rank $N_c*N_v*N_k$. Thus, the full BSE calculation is computationally only about twice as expensive as the TDA.

8.2 Input

INCAR

SYSTEM    = LiF BSE beyond TDA

ISMEAR    = 0      ! Gaussian smearing
SIGMA     = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT     = 500    ! energy cutoff

NBANDS    = 64     !  consistent with the ground state

! BSE
ALGO      = BSE    ! use BSE
ANTIRES   = 2      ! beyond TDA
NBANDSV   = 5      ! number of conduction bands in BSE
NBANDSO   = 5      ! number of valence bands in BSE
CSHIFT    = 0.2    ! complex shift
OMEGAMAX  = 40     ! max. energy
LHARTREE = .TRUE.  ! RPA diagrams
LADDER   = .TRUE.  ! ladder diagrams

8.3 Calculation

cd $TUTORIALS/BSE/e08_LiF_FBSE
cp $TUTORIALS/BSE/e04_LiF_GW/{WAVECAR,WAVEDER,W*.tmp} .
mpirun -np 4 vasp_std

Plot the calculated full BSE dielectric function with py4vasp

In [11]:
import py4vasp
import numpy as np
import os

E, ε = np.loadtxt(os.path.expanduser("./e09_LiF_TDHF_shifted_grid/expt.dat")).T

ipa = py4vasp.Calculation.from_path("./e05_LiF_IP")
rpa = py4vasp.Calculation.from_path("./e06_LiF_RPA")
tda = py4vasp.Calculation.from_path("./e07_LiF_BSE")
full = py4vasp.Calculation.from_path("./e08_LiF_FBSE")
(
    py4vasp.plot(E, ε, "Expt") +
    ipa.dielectric_function.plot("BSE(Im)").label("BSE-IPA") +
    rpa.dielectric_function.plot("BSE(Im)").label("BSE-RPA") +
    tda.dielectric_function.plot("BSE(Im)").label("BSE-TDA") +
    full.dielectric_function.plot("BSE(Im)").label("BSE-FULL")
)

8.4 Questions

  1. What is the difference between TDA and BSE?
  2. When does one need to use approximations beyond TDA?

Good job! You have finished Part 2!

Go to Top $\uparrow$