# Part 3: Efficient Brillouin zone sampling and analysis of the excitons¶

### 9 Efficient Brillouin zone sampling$\uparrow$¶

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

• Use shifted grids to speed up converged optical absorption calculations

BSE/TDHF calculations become very demanding if many k points are included because we need to diagonalize a matrix with the rank $N_k\times N_v\times N_c$. In this tutorial, you will learn to approximate a dense k-point mesh by separate calculations of coarser shifted meshes. In particular, you will calculate the TDDDH absorption spectrum of LiF with a $16\times16\times16$ k-points mesh using 8 calculations with a coarse $4\times4\times4$ grid.

Perform a TDDDH calculation averaged over multiple shifted in the following steps:

• Determine the irreducible k points $\mathbf{k}_n^{ir}$ with the corresponding weights $W_n$ from a DFT calculation with $4\times4\times4$ mesh.
• Perform a TDDDH calculation shifted to each of $\mathbf{k}_n^{ir}$ and extract the dielectric function
• Calculate the resulting dielectric function by adding the results of the calculations with the shifted grids weighted by $W_n$.

Perform a standard PBE calculation with Γ centered $4\times4\times4$ mesh for LiF.

##### 9.1.1 Input¶

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


INCAR.DFT

SYSTEM  = LiF DFT

ISYM     = -1    ! no symmetry in BZ
ISMEAR  = 0      ! Gaussian smearing
SIGMA   = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT   = 450    ! energy cutoff
NBANDS  = 12
KPAR    = 4


INCAR.DDH

SYSTEM   = LiF 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.47   ! the range-separation parameter
AEXX     = 0.53   ! the fraction of exact exchange
ISYM     = -1     ! no symmetry in BZ
NBANDS   = 12
LOPTICS  = .TRUE.
KPAR     = 4
CSHIFT   = 0.2


INCAR.TDHF

SYSTEM   = LiF TDHF exact exchange
LMODELHF = .TRUE. ! range-separated hybrid functional
ISMEAR   = 0      ! Gaussian smearing parameter
SIGMA    = 0.01   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT    = 450    ! energy cutoff BZ
LHFCALC  = .TRUE. ! exact exchange
LMODELHF = .TRUE. ! range-separated hybrid functional
HFSCREEN = 1.47   ! the range-separation parameter
AEXX     = 0.53   ! the fraction of exact exchange
LFXC     = .TRUE. ! local fxc
ISYM     = -1     ! no symmetry in BZ
NBANDS   = 12     ! total number of considered bands
ALGO     = TDHF   ! use TDHF algorithm
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 = 25     ! max. energy
NEDOS    = 2001   ! number of points in dielectric function


KPOINTS

k-mesh
0
Gamma
4 4 4
0 0 0


POTCAR

GW pseudopotentials of Li and F

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

cd $TUTORIALS/BSE/e09_LiF_TDHF_shifted_grid/DFT mpirun -np 4 vasp_std  Collect the irreducible k points$\mathbf{k}_n^{ir}$and weights$W_n$in the POINTS file grep -A11 "irreducible k-points:" OUTCAR | tail -8 > POINTS  #### 9.2 Task: Computing the dielectric function¶ Compute the optical absorption of LiF in the TDDDH approximation using multiple shifted grids. The next step will require to perform 8 TDDDH calculations, which might take around 10 min. Thus, we recommend to navigate to this example's directory and start the script. Then continue to read the description and to inspect the script. cd$TUTORIALS/BSE/e09_LiF_TDHF_shifted_grid
cp $TUTORIALS/BSE/e09_LiF_TDHF_shifted_grid/DFT/POINTS . ./run.sh  The exciton binding energy converges linearly with respect to the distance between two adjacent k points$\mathbf{q} = \mathbf{k} - \mathbf{k}'$. Thus, although the averaging over multiple grids can be a reasonable approximation for the spectra, the exciton binding energy is overestimated and has to be done with the explicit Brillouin-zone sampling. Another common approach to improving the convergence of the spectra with respect to the k points is based on a small asymmetric shift. This shift breaks the symmetry of the regular k-point grid and includes larger number of nonequivalent points as compared to the Γ-centered grid. A commonly used shift$\mathbf{s}=\{1/12,3/12,5/12\}$breaks most symmetry relations and includes the maximum number of nonequivalent points. run.sh #!/bin/bash VASP="mpirun -np 4 vasp_std" for count in {1..8}; do echo "Calculating k-point:"$count
mkdir -p "$count" cp POSCAR POTCAR$count/
sed 4q KPOINTS > $count/KPOINTS shift=$(sed -n "$count"p POINTS) echo$shift >> $count/KPOINTS weight=$(echo $shift | awk '{print$NF}')

cd "$count" cp ../INCAR.DFT ./INCAR$VASP

cp ../INCAR.DDH ./INCAR
$VASP cp ../INCAR.TDHF ./INCAR$VASP
../extract_optics.sh vasprun.xml
awk -v w=$weight '{print$1,$2*w}' optics.dat > optics_weighted.dat cd .. done paste {1..8}/optics_weighted.dat | awk '{print$1,($2+$4+$6+$8+$10+$12+$14+$16)/24};' > optics.dat


Plot the resulting dielectric function

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

# account for phonon renormalization
E_bse -= 1.15

plot((E_bse, ε_bse, "BSE"), (E_expt, ε_expt, "Expt"))


LiF has a large zero-phonon renormalization of 1.15 eV, which is not included in the band gap calculation. To simplify the comparison, we shifted the BSE spectrum by 1.15 eV. The overall agreement with experiment is good.

#### 9.3 Questions¶

• Can the approximation of the multiple shifted grids be used for finding accurate exciton binding energies?

### 10 Exciton analysis$\uparrow$¶

It can be useful to inspect properties of a particular excitonic state. In VASP it is possible to write the lowest NBSEEIG eigenvectors into a separate file BSEFATBAND. A convenient way to visually represent an eigenvector is a band structure plot where the coupling coefficients are depicted as circles at the corresponding bands and k points (fatband).

In order to write the BSEFATBAND file, the number of eigenvectors has to be provided NBSEEIG.

Here we are going to plot and investigate the fatband structure for diamond.

#### Task 10.1: Ground state calculation with DDH¶

In this tutorial we will use the range-separated 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

The input files to run this example are prepared at $TUTORIALS/BSE/e10_C_FATBAND/DFT. 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 8 8 8 0 0 0  DFT/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 4 empty orbitals DDH/INCAR for the DDH calculation of the unoccupied orbitals 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.26 ! the range-separation parameter AEXX = 0.175 ! the fraction of exact exchange NBANDS = 8 ! maybe even more bands LOPTICS = T ! write WAVEDER  POTCAR GW pseudopotentials of C. ##### 10.1.2 DFT¶ Navigate to the DFT directory and run VASP by entering the following: cd$TUTORIALS/BSE/e10_C_FATBAND/DFT
mpirun -np 4 vasp_std

##### 10.1.3 DDH¶

Now copy WAVECAR to the DDH directory and perform a DDH calculation with 4 empty bands and write WAVEDER file

cd $TUTORIALS/BSE/e10_C_FATBAND/DDH cp ../DFT/WAVECAR . mpirun -np 4 vasp_std  Find the direct and fundamental band gaps of diamond. #### Task 10.2: TDDDH calculation¶ Compute the TDDDH spectrum of diamond and write the eigenvectors for the first 10 excitons. Now that the DDH ground state is obtained, we can perform the TDDDH calculation. ##### 10.2.1 Input¶ The input files to run this example are prepared at $TUTORIALS/BSE/e10_C_FATBAND/.

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.26   ! the range-separation parameter
AEXX     = 0.175  ! the fraction of exact exchange
LFXC     = .TRUE. ! local fxc

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
NBSEEIG  = 10     ! number of BSE/TDHF eigenvectors


Navigate to the working directory and run VASP by entering the following:

cd $TUTORIALS/BSE/e10_C_FATBAND/ cp$TUTORIALS/BSE/e10_C_FATBAND/DDH/{WAVECAR,WAVEDER} .
mpirun -np 4 vasp_std


The BSEFATBAND file has to following structure

  rank of the BSE matrix                NBSEEIG
1BSE eigenvalue    E_BSE      IP-eigenvalue:    E_IP
Kx Ky Kz E_h E_e Abs(X_BSE)/W_k NB_h NB_h Re(X_BSE)+ i*Im(X_BSE)


where E_BSE and E_IP are the BSE and IP transition energies, Kx Ky Kz the k point coordinates, E_h and E_e the hole and electron eigenvalues, X_BSE the component of the eigenvector, W_k the weight of the k point, and NB_h NB_e the hole and electron orbital numbers.

For example, if we would like to make a fatband plot of the first bright exciton, we should first find it in the vasprun.xml file. The first transition with non-zero oscillator strength is the 4th state.

  <varray name="opticaltransitions" >
<v>     7.4041     0.0000 </v>
<v>     7.4041     0.0000 </v>
<v>     7.4041     0.0000 </v>
<v>     7.4057     8.0769 </v>
<v>     7.4057     8.0769 </v>
<v>     7.4057     8.0769 </v>
<v>     7.4107     0.0000 </v>
<v>     7.4107     0.0000 </v>
<v>     7.4257     0.0000 </v>
<v>     8.4266     0.0000 </v>
<v>     8.4266     0.0000 </v>
<v>     8.4408     0.0000 </v>
<v>     8.4408     0.0000 </v>
<v>     8.4408     0.0000 </v>
<v>     8.4551     0.0000 </v>
<v>     8.7107    21.6166 </v>
<v>     8.7107    21.6161 </v>
<v>     8.7107    21.6157 </v>
<v>     8.7159     0.0000 </v>
<v>     8.7159     0.0000 </v>


Hence we should plot the eigenstate 4BSE in the BSEFATBAND:

     4BSE eigenvalue    7.40574382      IP-eigenvalue:    7.82309365
0.00000  0.00000  0.00000   -12.9886220  19.6504559    0.0524972  1  5  -0.000103 +i*   0.000000
0.00000  0.00000  0.00000    11.8273623  19.6504559  413.2335502  2  5  -0.788790 +i*   0.170924
0.00000  0.00000  0.00000    11.8273623  19.6504559   27.9545231  3  5  -0.021670 +i*  -0.050114
0.00000  0.00000  0.00000    11.8273623  19.6504559   12.7380207  4  5  -0.009768 +i*  -0.022881
0.00000  0.00000  0.00000   -12.9886220  19.6504559    0.0056932  1  6   0.000008 +i*  -0.000007
0.00000  0.00000  0.00000    11.8273623  19.6504559   27.7696637  2  6  -0.032488 +i*   0.043431
0.00000  0.00000  0.00000    11.8273623  19.6504559  241.1676125  3  6   0.425077 +i*   0.202927
0.00000  0.00000  0.00000    11.8273623  19.6504559   10.3048824  4  6   0.018200 +i*   0.008593
0.00000  0.00000  0.00000   -12.9886220  19.6504559    0.0023290  1  7   0.000003 +i*  -0.000004
...


Knowing the IP energy of this transition and its BSE eigenvalue, we can find the binding energy of the first bright exciton 0.42 eV, which is strongly overestimated due to the sparse k-point grid.

We are going to make the plot with the axis

|k-point|  electron/hole  Abs(X_BSE)
eigenvalue


for the k-points along the Γ-L and Γ-X directions.

In :
import py4vasp
import numpy as np

calc = py4vasp.Calculation.from_path("./e10_C_FATBAND")

# calculate k point distances for band plot
X_Γ = calc.kpoint.path_indices((0,0,0), (0.5,0.5,0.0))[::-1]
Γ_L = calc.kpoint.path_indices((0,0,0), (0.5,0.5,0.5))
X_Γ_L = np.concatenate((X_Γ, Γ_L[1:]))  # remove duplicate Γ point
dists = np.linalg.norm(kpoints[X_Γ_L], axis=1)
dists[:len(X_Γ)] *= -1  # plot Γ-X in negative range

# only one spin channel
spin = 0
bse_index = data["bse_index"][spin,X_Γ_L]
bands = data["bands"][spin,X_Γ_L].T
vbmin = data["first_valence_band"][spin]
cbmin = data["first_conduction_band"][spin]

# select one particular eigenstate
λ = 3
scale = 2
fatband = scale * np.abs(data["fatbands"][λ,bse_index])

# generate plot
def plot_pair(valence, conduction):
return py4vasp.plot(
dists,
np.array((bands[vbmin + valence], bands[cbmin + conduction])),
name=f"valence band {valence + 1} – conduction band {conduction + 1}",
width=fatband[:,conduction,valence],
marker="o",
)

graph = (
py4vasp.plot(
dists, bands, "bands",
xticks={dists.min(): "X", 0.0: "Γ", dists.max(): "L"},
ylabel="$E - E_F (\mathrm{eV})$",
marker="o",
color="black",
) +
plot_pair(1, 0) + plot_pair(2, 1) + plot_pair(3, 2) + plot_pair(1, 1) + plot_pair(2, 0)
)
graph.to_plotly().update_layout(width=600)

NoData: Could not find data in output, please make sure that the provided input
should produce this data and that the VASP calculation already finished.
Also check that VASP did not exit with an error.


In the plot above, the maximum contribution comes from the band maxima at the Γ point near the direct band gap. The fact that the contribution at the adjacent k points is rather small means that this particular excitonic state is well localized in the k space.

#### 10.3 Questions:¶

• What information can be found in the BSEFATBAND file?