Part 3: Efficient Brillouin zone sampling and analysis of the excitons¶
9 Efficient Brillouin zone sampling ¶
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$.
9.1 Task¶
Perform a standard PBE calculation with Γ centered $4\times4\times4$ mesh for LiF.
9.1.1 Input¶
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
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
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
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
k-mesh
0
Gamma
4 4 4
0 0 0
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
from py4vasp import plot
import numpy as np
import os
E_bse, ε_bse = np.loadtxt(os.path.expanduser("./e09_LiF_TDHF_shifted_grid/optics.dat")).T
E_expt, ε_expt = np.loadtxt(os.path.expanduser("./e09_LiF_TDHF_shifted_grid/expt.dat")).T
# 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 ¶
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:
- DDH ground state
- DDH with exact diagonalization to compute unoccupied bands
10.1.1 Input¶
The input files to run this example are prepared at $TUTORIALS/BSE/e10_C_FATBAND/DFT.
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
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
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/.
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
x y circle radius
for the k-points along the Γ-L and Γ-X directions.
import py4vasp
import numpy as np
calc = py4vasp.Calculation.from_path("./e10_C_FATBAND")
data = calc.fatband.read()
# calculate k point distances for band plot
kpoints = calc.kpoint.read()["coordinates"]
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)
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?