# Part 1: Optical absorption of diamond carbon¶

### 0 Bethe-Salpeter-equation formalism$\uparrow$¶

The formalism of the Bethe-Salpeter equation (BSE) allows to include the excitonic effects in the calculation of the dielectric function. The excitation energies are the eigenvalues $\omega_\lambda$ of the following linear problem

$$\tag{1} \left(\begin{array}{cc} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & \mathbf{A}^* \end{array}\right)\left(\begin{array}{l} \mathbf{X}_\lambda \\ \mathbf{Y}_\lambda \end{array}\right)=\omega_\lambda\left(\begin{array}{cc} -\mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{1} \end{array}\right)\left(\begin{array}{l} \mathbf{X}_\lambda \\ \mathbf{Y}_\lambda \end{array}\right)~.$$

The matrices $A$

$$\tag{2} A_{vc}^{v'c'} = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'} + \langle cv'|V|vc'\rangle - \langle cv'|W|c'v\rangle$$

and $B$

$$\tag{3} B_{vc}^{v'c'} = \langle vv'|V|cc'\rangle - \langle vv'|W|c'c\rangle.$$

describe the transitions between the occupied $v,v'$ and unoccupied $c,c'$ states.

Two terms describe the interaction in the BSE Hamiltonian: the bare Coulomb interaction $V$ and and the screened one $W$. The energies and orbitals are usually obtained in a $GW$ calculation.

In reciprocal space, the matrix $A$ is written as

$$\tag{4} \begin{equation} \begin{split} A_{vc\mathbf{k}}^{v'c'\mathbf{k}'} & = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'}\delta_{\mathbf{kk}'}+ \frac{2}{\Omega}\sum_{\mathbf{G}\neq0}\bar{V}_\mathbf{G}(\mathbf{q})\langle c\mathbf{k}|e^{i\mathbf{Gr}}|v\mathbf{k}\rangle\langle v'\mathbf{k}'|e^{-i\mathbf{Gr}}|c'\mathbf{k}'\rangle \\ & -\frac{2}{\Omega}\sum_{\mathbf{G,G}'}W_{\mathbf{G,G}'}(\mathbf{q},\omega)\delta_{\mathbf{q,k-k}'} \langle c\mathbf{k}|e^{i(\mathbf{q+G})}|c'\mathbf{k}'\rangle \langle v'\mathbf{k}'|e^{-i(\mathbf{q+G})}|v\mathbf{k}\rangle, \end{split} \end{equation}$$

where $\Omega$ is the cell volume, $\bar{V}$ is the bare Coulomb potential without the long-range part

$$\tag{5} \bar{V}_{\mathbf{G}}(\mathbf{q})=\begin{cases} 0 & \text { if } G=0 \\ V_{\mathbf{G}}(\mathbf{q})=\frac{4 \pi}{|\mathbf{q}+\mathbf{G}|^2} & \text { else } \end{cases}~,$$

and the screened Coulomb potential

$$\tag{6} W_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega)=\frac{4 \pi \epsilon_{\mathbf{G}, \mathbf{G}^{\prime}}^{-1}(\mathbf{q}, \omega)}{|\mathbf{q}+\mathbf{G}|\left|\mathbf{q}+\mathbf{G}^{\prime}\right|}.$$

Here, the dielectric function $\epsilon_\mathbf{G,G'}(\mathbf{q})$ describes the screening in $W$ within the random-phase approximation (RPA)

$$\tag{7} \epsilon_{\mathbf{G}, \mathbf{G}^{\prime}}^{-1}(\mathbf{q}, \omega)=\delta_{\mathbf{G}, \mathbf{G}^{\prime}}+\frac{4 \pi}{|\mathbf{q}+\mathbf{G}|^2} \chi_{\mathbf{G}, \mathbf{G}^{\prime}}^{\mathrm{RPA}}(\mathbf{q}, \omega).$$

Although the dielectric function is frequency-dependent, the static approximation $W_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega=0)$ is considered a standard for practical BSE calculations. In a simplified notation, we refer to the first term of $A$ as diagonal $D$, the second term as exchange $K^x$, and the third term as direct $K^D$ contribution. The expression for $B_{vc\mathbf{k}}^{v'c'\mathbf{k}'}$ is analogous to Eq. (4) replacing the corresponding indices.

The Tamm-Dancoff approximation (TDA) neglects the coupling between excitations and de-excitations, i.e., the off-diagonal elements $B$ and $B^*$. Hence, the TDA reduces Eq. (1) to

$$\tag{8} AX_\lambda=\omega_\lambda X_\lambda~.$$

The macroscopic dielectric function results from the eigenvalues $\omega_\lambda$ and eigenvectors $X_\lambda$ found via the BSE

$$\tag{9} \epsilon_M(\mathbf{q},\omega)= 1+\lim_{\mathbf{q}\rightarrow 0}v(q)\sum_{\lambda} \left|\sum_{c,v,k}\langle c\mathbf{k}|e^{i\mathbf{qr}}|v\mathbf{k}\rangle X_\lambda^{cv\mathbf{k}}\right|^2 \times \left(\frac{1}{\omega_\lambda - \omega - i\delta} + \frac{1}{\omega_\lambda+\omega + i\delta}\right)~.$$

In the optical limit ($q \rightarrow 0$), the Taylor expansion $e^{i\mathbf{qr}}\approx 1 + i\mathbf{qr} + O(\mathbf{q}^2)$ yields $\langle c\mathbf{k}|e^{i\mathbf{qr}}|v\mathbf{k}\rangle \approx i\mathbf{q}\langle c\mathbf{k}|\mathbf{r}|v\mathbf{k}\rangle \propto \mathbf{q}\langle c\mathbf{k}|\nabla|v\mathbf{k}\rangle$.

Time-dependent density functional theory (TDDFT) provides an alternative approach to obtain the dielectric function. In the Casida formulation, the TDDFT theory leads to the same linear problem as Eq. (1). However, the interaction between electron and holes is described by the exchange-correlation kernel $f_{xc}$

$$\tag{10} A_{vc}^{v'c'} = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'} + \langle cv'|V|vc'\rangle - \langle cv'|f_{\rm xc}|c'v\rangle~.$$$$\tag{11} B_{vc}^{v'c'} = \langle vv'|V|cc'\rangle - \langle vv'|f_{\rm xc}|c'c\rangle~.$$

The kernel $f_{\rm xc}$ is defined from the exchange-correlation functional $$\tag{12} f_{\rm xc}(\mathbf{r,r'})=\frac{\delta v_{\rm xc}(\mathbf{r})}{\delta \rho(\mathbf{r'})}~.$$ Hence, the electron-hole interaction in Eq. (10) reads $$\tag{13} \langle cv'|f_{\rm xc}|c'v\rangle=\gamma\langle cv'|V|c'v\rangle+(1-\gamma)\langle cv'|f^{\rm ALDA}_{\rm x}|vc'\rangle + \langle cv'|f^{\rm ALDA}_{\rm c}|vc'\rangle~,$$ where $\gamma$ describes the screening of the attractive interaction and is determined by the fraction of Fock exchange of an exchange-correlation hybrid functional. This term can be a constant or a function, depending on the approximation. $f_{\rm x}^{ALDA}$ and $f_{\rm c}^{ALDA}$ are the local adiabatic exchange and correlation kernels.

To perform BSE/TDHF calculations, one needs quasiparticle orbitals $\psi_{n,k}$, orbital derivatives $\nabla_k \psi_{n,k}$, and energies $\varepsilon_{n,k}$ from a preceding ground-state calculation. In VASP these quantities are stored in WAVECAR and WAVEDER. BSE requires in addition the screened Coulomb potential $W_{G,G'}(\omega=0)$ produced in $GW$ calculations. VASP write $W$ to WFULLxxxx.tmp and Wxxxx.tmp files.

### 1 Preparatory ground-state $G_0W_0$ calculation$\uparrow$¶

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
• Plot the RPA dielectric function obtained in the GW calculation

Perform a $G_0W_0$ calculation based on the PBE ground state for diamond.

$GW$ calculations typically require a large number of unoccupied states for converged quasiparticle energies. The convergence with respect to the number of states, k points, and frequency points has to be thoroughly checked. For relatively small cells, such as our case here, it can be faster to perform an exact diagonalization of the Hamiltonian.

In this task, we use the $G_0W_0$ approximation as a starting point for our BSE calculations. 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

#### 1.2 Input¶

The input files to run this example are prepared at $TUTORIALS/BSE/e01_C_GW. We take the unit cell with the experimental lattice parameter. POSCAR Diamond 3.567 0.5 0.5 0.0 0.0 0.5 0.5 0.5 0.0 0.5 C 2 Cartesian 0.00 0.00 0.00 C 0.25 0.25 0.25 C  We use a relatively small number of k-points, but it should be sufficient to reach the convergence of the band gap within 0.1 eV. KPOINTS k-mesh 0 Gamma 6 6 6 0 0 0  ground-state/INCAR for the PBE ground state calculation SYSTEM = C DFT ISMEAR = 0 ! Gaussian smearing SIGMA = 0.01 ! smearing in eV ENCUT = 450 ! energy cutoff  We include 60 empty orbitals, which should be sufficient to get the band gap within 0.1 eV. unoccupied-states/INCAR for the PBE calculation of the unoccupied orbitals SYSTEM = C DFT unoccupied states ISMEAR = 0 ! Gaussian smearing SIGMA = 0.01 ! smearing in eV, small sigma is required to avoid partial occupancies ENCUT = 450 ! energy cutoff NBANDS = 64 ! include 60 unoccupied bands ALGO = Exact ! use exact diagonalization of the Hamiltonian NELM = 1 ! since we are already converged, stop after one step LOPTICS = T ! write orbitals derivatives to WAVEDER CSHIFT = 0.6 ! broadening of the dielectric function  We use 50 frequency points which is smaller than the default but satisfies our convergence criteria of 0.1 eV for the band gap. INCAR for the$G_0W_0$calculation SYSTEM = C G0W0 ISMEAR = 0 ! Gaussian smearing SIGMA = 0.01 ! smearing in eV, small sigma is required to avoid partial occupancies ENCUT = 450 ! energy cutoff NBANDS = 64 ! the same as in the previous step ALGO = G0W0 ! use the evGW algorithm NELM = 1 ! only one iteration for G0W0 KPAR = 4 ! number of k-points treated in parallel NOMEGA = 50 ! the default is 100  The potentials recommended for the GW calculation are labeled "_GW" POTCAR GW pseudopotentials of C. #### 1.3 Calculation¶ ##### 1.3.1 DFT unoccupied bands¶ Open a terminal, navigate to this example's directory by entering following command cd$TUTORIALS/BSE/e01_C_GW/ground-state


and run VASP inside this directory on four CPU cores by executing the command

mpirun -np 4 vasp_std


We can use py4vasp to obtain the direct and the fundamental band gaps. The direct band gap is found as the energy difference between the highest occupied and the lowest unoccupied bands at the same k-point, while the fundamental gap is found as the energy difference between the valence band maximum (VBM) and the conduction band minimum (CBM), i.e., energies at different k-points.

The fundamental band gap of diamond from experiment is known to be 5.85 eV (Phys. Rev. Materials 2, 073803 (2018)) supporting the well-known fact that the PBE functional strongly underestimates the band gap of diamond.

The provided experimental value is corrected for the zero-phonon renormalization (ZPR), which effectively "opens" the experimental band gap. These effects are naturally present in the experiment but not taken into account in our calculations. Hence, we correct the experimental band gap to be able to make a meaningful comparison.

In :
import py4vasp

def fgap(calc):
homo = gs.band.to_frame().query('occupations > 0.5').max()['bands']
lumo = gs.band.to_frame().query('occupations < 0.5').min()['bands']
return  lumo - homo

def dgap(calc):
energies = gs.band.to_dict()
dir_gaps = [min(b[o < 0.5]) - max(b[o > 0.5]) for b,o in zip(energies["bands"],energies["occupations"])]
return min(dir_gaps)

gs = py4vasp.Calculation.from_path("./e01_C_GW/ground-state")

print("PBE fundamental band gap for C: {0:.2f}".format(fgap(gs)))
print("PBE direct band gap for C: {0:.2f}".format(dgap(gs)))


PBE fundamental band gap for C: 4.16
PBE direct band gap for C: 5.60

##### 1.3.2 DFT unoccupied bands¶

Now, navigate to the directory with the input for calculating the unoccupied states, copy the WAVECAR file from the previous step, and perform a calculation by entering the following commands:

cd $TUTORIALS/BSE/e01_C_GW/unoccupied-states cp ../ground-state/WAVECAR . mpirun -np 4 vasp_std  In the OUTCAR file find a block of lines starting with the following lines and make sure that the right number of bands was calculated.  k-point 1 : 0.0000 0.0000 0.0000 band No. band energies occupation  Make sure that the WAVEDER file is present in the directory. Since LOPTICS is enabled in the INCAR, VASP calculates the dielectric function in the independent particle approximation, which can be extracted using py4vasp In : import py4vasp ipa = py4vasp.Calculation.from_path("./e01_C_GW/unoccupied-states") ipa.dielectric_function.plot("Im")  ##### 1.3.3$G_0W_0$calculation¶ Now we can perform$G_0W_0$using the WAVECAR and WAVEDER from the previous step. Navigate to the folder with the$GW$input files and run the code cd$TUTORIALS/BSE/e01_C_GW
cp unoccupied-states/{WAVECAR,WAVEDER} .
mpirun -np 4 vasp_std


The QP eigenvalues are written into the OUTCAR file.

k-point   1 :       0.0000    0.0000    0.0000
band No.  KS-energies  QP-energies   sigma(KS)   V_xc(KS)     V^pw_x(r,r')   Z       occupation Imag(sigma)

1      -8.7939     -10.5231     -17.0439     -14.7585     -25.4072       0.7566    2.0000       1.4530
2      12.6743      11.8907     -17.7523     -16.8014     -19.2036       0.8241    2.0000       0.0666
3      12.6743      11.8907     -17.7523     -16.8014     -19.2036       0.8241    2.0000       0.0666
4      12.6743      11.8907     -17.7523     -16.8014     -19.2036       0.8241    2.0000       0.0666
5      18.2756      19.2915     -14.2721     -15.5097      -9.3175       0.8208    0.0000      -0.1051
6      18.2756      19.2915     -14.2721     -15.5097      -9.3175       0.8208    0.0000      -0.1051


We find the direct band gap of 7.40 eV and the fundamental band gap of 5.58 eV, which is quite close to the experimental value of 5.85 eV.

Now, we can plot the RPA dielectric function

In :
import py4vasp

ipa = py4vasp.Calculation.from_path("./e01_C_GW/unoccupied-states")
rpa = py4vasp.Calculation.from_path("./e01_C_GW")
(
ipa.dielectric_function.plot("Im").label("IPA") +
rpa.dielectric_function.plot("RPA(Im)").label("RPA")
)


#### 1.4 Questions¶

1. What is the difference between direct and fundamental band gaps?
2. What is the difference between the dielectric function calculated with the LOPTICS tag and the one obtained in the $GW$ calculation?

### 2 Optical absorption in BSE$\uparrow$¶

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

• Perform calculations of the optical absorption including the excitonic effects

Compute the optical absorption spectrum using the QP and the screened potential calculated in the previous step.

Now that we have WAVECAR, WAVEDER, and W*.tmp files, we can perform the BSE calculation in the Tamm-Dancoff approximation. In the BSE calculation, first the Hamiltonian is constructed using the QP orbitals and the screened potential in the RPA approximation. Then, the eigenvalues and the eigenvectors of Eq. (8) are found. Finally, the BSE dielectric function is calculated using Eq. (9).

The input files to run this example are prepared at $TUTORIALS/BSE/e02_C_BSE. The BSE calculation is performed by setting ALGO = BSE. We construct the BSE Hamiltonian for all four occupied and four empty states, which corresponds to the highest transition energy of ~37 eV. The NBANDS tag has to be set to the same value as in the preceding GW calculation and the orbitals derivatives, i.e., WAVECAR and WAVEDER. INCAR SYSTEM = C BSE ISMEAR = 0 ! Gaussian smearing SIGMA = 0.01 ! smearing in eV, small sigma is required to avoid partial occupancies ENCUT = 450 ! energy cutoff NBANDS = 64 ! should be the same as in the GS ! BSE ALGO = BSE ! BSE calculation ANTIRES = 0 ! Tamm-Dancoff approximation NBANDSV = 4 ! number of conduction bands in BSE NBANDSO = 4 ! number of valence bands in BSE CSHIFT = 0.6 ! complex shift/broadening  POSCAR Diamond 3.567 0.5 0.5 0.0 0.0 0.5 0.5 0.5 0.0 0.5 C 2 Cartesian 0.00 0.00 0.00 C 0.25 0.25 0.25 C  KPOINTS k-mesh 0 Gamma 6 6 6 0 0 0  POTCAR GW pseudopotentials of C. #### 2.3 Calculation¶ Navigate to this example's directory and run VASP by entering the following: cd$TUTORIALS/BSE/e02_C_BSE
cp $TUTORIALS/BSE/e01_C_GW/{WAVECAR,WAVEDER,W*.tmp} . mpirun -np 4 vasp_std  The calculated optical transitions with the corresponding oscillator strengths can be found in vasprun.xml  <varray name="opticaltransitions" > <v> 6.8588 0.0000 </v> <v> 6.8588 0.0000 </v> <v> 6.8588 0.0000 </v> <v> 6.8638 8.4561 </v> <v> 6.8638 8.4561 </v> <v> 6.8638 8.4561 </v> <v> 6.8732 0.0000 </v> <v> 6.8732 0.0000 </v> <v> 6.9028 0.0000 </v> <v> 8.4940 0.0000 </v> <v> 8.4940 0.0000 </v> <v> 8.5027 0.0000 </v> <v> 8.5027 0.0000 </v> <v> 8.5027 0.0000 </v> <v> 8.5475 0.0000 </v> <v> 8.7058 15.9682 </v>  the BSE dielectric function is also stored in vasprun.xml  <dielectricfunction> <imag> <array> <dimension dim="1">gridpoints</dimension> <field>energy</field> <field>xx</field> <field>yy</field> <field>zz</field> <field>xy</field> <field>yz</field> <field>zx</field> <set> <r> 0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 </r> <r> 0.0374 0.0020 0.0020 0.0020 0.0000 0.0000 0.0000 </r> <r> 0.0749 0.0040 0.0040 0.0040 0.0000 0.0000 0.0000 </r> <r> 0.1123 0.0061 0.0061 0.0061 0.0000 -0.0000 0.0000 </r> <r> 0.1498 0.0081 0.0081 0.0081 0.0000 -0.0000 0.0000 </r>  The average dielectric function over xx, yy, zz can be extracted using the following AWK script awk 'BEGIN{i=0} /<dielectricfunction>/,\ /<\/dielectricfunction>/ \ {if ($1=="<r>") {a[i]=$2 ; b[i]=($3+$4+$5)/3 ; c[i]=$4 ; d[i]=$5 ; i=i+1}} \
END{for (j=0;j<i/2;j++) print a[j],b[j],b[j+i/2]}' vasprun.xml > optics.dat


It can also be useful to check the OUTCAR file to make sure that the calculations are performed as expected. The bands included in the BSE Hamiltonian can be found in the line

Bands included in the BSE spin= 1 VB(min)= 1 VB(max)= 4 CB(min)= 5 CB(max)= 8

To get the rank of the BSE matrix find the following line

BSE (scaLAPACK) attempting allocation of 0.048 Gbyte rank= 3456

Tip: If you get all zeros in the BSE dielectric function, make sure that WAVECAR and WAVEDER files were read successfully, i.e., find these lines in stdout

the WAVECAR file was read successfully

Now we can compare the BSE dielectric function

In :
import py4vasp

ipa = py4vasp.Calculation.from_path("./e01_C_GW/unoccupied-states")
rpa = py4vasp.Calculation.from_path("./e01_C_GW/")
bse = py4vasp.Calculation.from_path("./e02_C_BSE/")
(
ipa.dielectric_function.plot("Im").label("IPA") +
rpa.dielectric_function.plot("RPA(Im)").label("RPA") +
bse.dielectric_function.plot("BSE(Im)").label("BSE")
)


#### 2.4 Questions¶

1. How does the e-h interaction affect the dielectric constant?

### 3 Optical absorption of diamond through TDDDH$\uparrow$¶

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

• Find parameters for the model dielectric function
• Perform calculations of the optical absorption via TDDFT including the excitonic effects

Find parameters for the model dielectric function and calculate the optical absorption spectrum via TDDFT

The $GW$ step is computationally demanding, however it is possible to achieve similar accuracy by using the hybrid-functional approach instead and circumvent the $GW$ calculation.

In the dielectric-dependent hybrid functional (DDH), the screening of Fock exchange is described by the function

$$\tag{14} \epsilon^{-1}(|\mathbf{q+G}|)=1-(1-\epsilon_\infty^{-1})e^{-(|\mathbf{q+G}|)^2/4\mu^2}.$$

The same function can be used for screening the attractive Coulomb interaction between electrons and holes. Such approximation is called TDDDH (see detail in Phys. Rev. Research 2, 032019 (2020)).

In this tutorial we will use the dielectric-dependent hybrid functional for the ground state calculation.

The DDH calculation can be performed in two steps:

1. DDH ground state
2. DDH with exact diagonalization to compute unoccupied bands

#### 3.2 Input¶

The input files to run this example are prepared at $TUTORIALS/BSE/e03_C_TDHF. The parameters can be found by fitting the model to the RPA dielectric function in the vasprun.xml of the GW calculation. In order to do that, you can run the following commands cd$TUTORIALS/BSE/e03_C_TDHF
awk ' /<varray name="epsilon_diag" >/,/<\/varray>/ \
{if ($1=="<v>") {print$2, $3}} '$TUTORIALS/BSE/e01_C_GW/vasprun.xml  > eps_diag.dat


The following script plots the RPA dielectric function with the parameters: $\varepsilon=5.7$ and $\mu=1.79$.

In :
import numpy as np
import os
from py4vasp import plot

def eps_model(g,mu,eps_inf):
return 1. - (1.-1./eps_inf)*np.exp(-g**2/(4.*mu**2))

mesh = np.linspace(0, G.max(), num=200)
model = eps_model(mesh, mu=1.70, eps_inf=5.7)

plot(G, ε, "RPA", marker="o") + plot(mesh, model, "Model", xlabel="$G (Å^{-1})$", ylabel="$ε^{-1}$")


POSCAR

Diamond
3.567
0.5 0.5 0.0
0.0 0.5 0.5
0.5 0.0 0.5
C
2
Cartesian
0.00 0.00 0.00 C
0.25 0.25 0.25 C


ground-state/INCAR

SYSTEM   = C Range-separated Dielectric-Dependent Hybrid functional

ISMEAR   = 0      ! Gaussian smearing
SIGMA    = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT    = 450    ! energy cutoff
LHFCALC  = .TRUE. ! exact exchange
LMODELHF = .TRUE. ! range-separated hybrid functional
HFSCREEN = 1.70   ! the range-separation parameter
AEXX     = 0.175  ! the fraction of exact exchange



unoccupied-states/INCAR

SYSTEM   = C RS-DDH unoccupied states

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

LHFCALC  = .TRUE. ! exact exchange
LMODELHF = .TRUE. ! range-separated hybrid functional
HFSCREEN = 1.70   ! the range-separation parameter
AEXX     = 0.175  ! the fraction of exact exchange

NBANDS   = 8      ! include 4 unoccupied bands

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   = C TDHF

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

LMODELHF = .TRUE. ! range-separated hybrid functional
HFSCREEN = 1.70   ! the range-separation parameter
AEXX     = 0.175  ! the fraction of exact exchange
LFXC     = .TRUE. ! local xc kernel

NBANDS   = 8     ! consistent with the ground state

ALGO     = TDHF   ! use TDHF algorithm
ANTIRES   = 0     ! Tamm-Dancoff approximation
NBANDSV  = 4      ! number of conduction bands in BSE
NBANDSO  = 4      ! number of valence bands in BSE
CSHIFT   = 0.6    ! complex shift/broadening


KPOINTS

k-mesh
0
Gamma
6 6 6
0 0 0


POTCAR

GW pseudopotentials of C

#### 3.3 Calculation¶

##### 3.3.1 DDH ground state¶

Navigate to this example's directory

cd $TUTORIALS/BSE/e03_C_TDHF/ground-state  copy the wave functions from the DFT calculation to speed up the convergence for the hybrid functional calculation cp$TUTORIALS/BSE/e01_C_GW/ground-state/WAVECAR .


and run VASP

mpirun -np 4 vasp_std


This yields the DDH ground state.

You can find in the OUTCAR file that the DDH direct gap is 7.5 eV and the fundamental gap is 5.6 eV.

##### 3.3.2 DDH unoccupied bands¶

As in the case with the PBE calculation, we will use the exact diagonalization algorithm for calculating unoccupied orbitals.

Navigate to this example's directory

cd $TUTORIALS/BSE/e03_C_TDHF/unoccupied-states  copy the wave functions from the previous DDH calculation cp$TUTORIALS/BSE/e03_C_TDHF/ground-state/WAVECAR .


and run VASP

mpirun -np 4 vasp_std


Now we have the WAVECAR and WAVEDER files and we can proceed with the TDDDH calculation.

##### 3.3.3 Optical absorption in TDDDH¶

Navigate to this example's directory

cd $TUTORIALS/BSE/e03_C_TDHF/  copy the wave functions and their derivatives calculated via DDH cp$TUTORIALS/BSE/e03_C_TDHF/unoccupied-states/{WAVECAR,WAVEDER} .


run VASP

mpirun -np 4 vasp_std


Make sure that the rank of the BSE Hamiltonian in OUTCAR is the same as in the BSE calculation (2.3).

Plot the dielectric function

In :
import py4vasp

ipa = py4vasp.Calculation.from_path("./e01_C_GW/unoccupied-states")
rpa = py4vasp.Calculation.from_path("./e01_C_GW/")
bse = py4vasp.Calculation.from_path("./e02_C_BSE")
tdhf = py4vasp.Calculation.from_path("./e03_C_TDHF")
(
ipa.dielectric_function.plot("Im").label("IPA") +
rpa.dielectric_function.plot("RPA(Im)").label("RPA") +
bse.dielectric_function.plot("BSE(Im)").label("BSE") +
tdhf.dielectric_function.plot("TDHF(Im)").label("TDDDH")
)


Compare the TDDDH and BSE spectra. You can see the excellent agreement between BSE and TDDDH. The sparse k-point grid does not allow us to see the difference between the RPA and the spectra with e-h interaction. However, in Fig. 2 of Phys. Rev. Lett. 107, 186401 (2011) you can see how important the excitonic effects are.

The dielectric screening in diamond is quite strong (the dielectric constant is ~5.7), hence the excitons are only weakly bound. Such excitons are called Wannier-Mott excitons. In part 2, we will study a system with a weak screening and strongly bound Frenkel excitons - LiF.

#### 3.4 Questions¶

1. How can one determine the parameters for the model dielectric function in Eq. (14)?
2. What is the main advantage of using TDDDH over BSE?