**Part 2: Optical absorption of LiF**¶

**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:

- DFT ground state
- DFT with exact diagonalization to compute unoccupied bands
- $G_0W_0$ calculation

**4.2 Input**¶

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

.

```
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
```

```
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
```

```
k-mesh
0
Gamma
6 6 6
0 0 0
```

*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**¶

- Is LiF a direct of indirect band gap system?
- 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.

```
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

```
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**¶

- Are excitonic effects accounted for in the IP approximation?
- 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`

.

```
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

```
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**¶

- Are excitonic effects accounted for in the RPA approximation?
- 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`

.

```
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

```
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**¶

- Does the Tamm-Dancoff approximation account for the excitonic effects?
- 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**¶

```
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

```
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**¶

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