Part 1: Optical absorption of diamond carbon¶
0 Bethe-Salpeter-equation formalism ¶
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 ¶
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
1.1 Task¶
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:
- DFT ground state
- DFT with exact diagonalization to compute unoccupied bands
- $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.
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.
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"
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.
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)))
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
import py4vasp
ipa = py4vasp.Calculation.from_path("./e01_C_GW/unoccupied-states")
ipa.dielectric_function.plot("Im")
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
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¶
- What is the difference between direct and fundamental band gaps?
- 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 ¶
At the end of this tutorial you will be able to:
- Perform calculations of the optical absorption including the excitonic effects
2.1 Task¶
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).
2.2 Input¶
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.
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
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
k-mesh
0
Gamma
6 6 6
0 0 0
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
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¶
- How does the e-h interaction affect the dielectric constant?
3 Optical absorption of diamond through TDDDH ¶
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
3.1 Task¶
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:
- DDH ground state
- 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$.
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))
G, ε = np.loadtxt(os.path.expanduser("./e03_C_TDHF/eps_diag.dat")).T
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}$")
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
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
k-mesh
0
Gamma
6 6 6
0 0 0
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
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¶
- How can one determine the parameters for the model dielectric function in Eq. (14)?
- What is the main advantage of using TDDDH over BSE?