Part 2: XAS K-edge of LiCl via Bethe-Salpeter equation

In this part of the tutorial, we will use the Many-Body Perturbation Theory (MBPT) to simulate X-ray absorption in rock salt LiCl. Within MBPT, the electron-hole interaction or the excitonic effects in the spectra can be accounted for via the Bethe-Salpeter equation (BSE).

As a first step, we need to calculate the quasiparticles and the dielectric tensor in the $GW$ approximation. Within the $GW$ approximation, the particles are not independent. Instead, electrons interact with each other as they propagate through the medium and induce a charge polarization. The electrons are said to be dressed by the polarization cloud they induce in their immediate environment. Such dressed electrons are called quasiparticles.

These quasiparticles (eigenstates and eigenenergies) are obtained by solving the following one-electron equation: $$ \tag{1} (T+V_{ext}+V_h)\psi_{n{\bf k}}({\bf r})+\int d{\bf r}\Sigma({\bf r},{\bf r}', E_{n{\bf k}})\psi_{n{\bf k}}({\bf r}') = E_{n{\bf k}}\psi_{n{\bf k}}({\bf r}). $$ The many-body effects in this equation are introduced via the self-energy $\Sigma({\bf r},{\bf r}',E_{n{\bf k}})$, which is a non-local and energy-dependent potential. In $GW$ the self-energy $\Sigma$, is approximated as $\Sigma=GW$, the product of a Green's function $G$ and the screened Coulomb interaction $W$, hence the name. The dielectric screening of the Coulomb interaction is commonly calculated within the random-phase approximation (RPA).

The quasiparticle equation should, in principle, be solved self-consistently. However, the most common approach is to stop after the first iteration, known as a one-shot $G_0W_0$ approximation. In the $G_0W_0$ approximation, the Green's function and the screened Coulomb interaction are calculated from KS eigenstates and eigenenergies, and the quasiparticle equation is solved only once, i.e., without any self-consistency.

For more details on the theoretical background of GW, check out the VASP Wiki.


3 Band gap of LiCl within the $G_0W_0$ approximation

By the end of this tutorial, you will be able to:

  • run a one-shot $GW$ calculation ($G_0W_0$)
  • obtain the $G_0W_0$ band gap
  • calculate the dielectric function in the RPA approximation

3.1.1 Task

Compute the quasiparticle energies and find the band gap of LiCl with a $G_0W_0$ calculation.

A G$_0$W$_0$ is performed in three steps:

  1. The first step, is a DFT-ground-state calculation.
  2. In the second step one restarts from the self-consistent solution (WAVECAR) of (1) and computes additional unoccupied KS eigenstates. This is done by a single exact diagonalization of the Hamilton matrix. These additional unoccupied KS states are needed to compute the dielectric screening properties. Furthermore, in this second step, we compute the change of the cell periodic part of the KS orbitals with respect to the Bloch vector $\bf{k}$ (this information is stored in the WAVEDER file).
  3. The third and final step is the actual $G_0W_0$ calculation. This step needs the WAVECAR and WAVEDER from (2) as a starting point.
3.1.2 Input files

The input files to run the DFT groundstate calculation are prepared at $TUTORIALS/XAS/e03_LiCl_G0W0/dft. Check them out!

POSCAR:
LiCl in the cubic structure.


rock salt LiCl
5.106
  0.0 0.5 0.5
  0.5 0.0 0.5
  0.5 0.5 0.0
Li Cl
1  1
Direct
  0.00 0.00 0.00 Li
  0.50 0.50 0.50 Cl

INCAR:
A DFT run using the default generalized-gradient-approximation functional of Perdew, Burke, and Ernzerhof (GGA-PBE).


SYSTEM = LiCl DFT ground state 

ISMEAR = 0      ! Gaussian smearing
SIGMA = 0.05    ! small value for smearing
EDIFF  = 1E-8   ! convergence criterion
ENCUT  = 300    ! energy cutoff

KPOINTS:
An equally spaced (6$\times$6$\times$6) Γ-centered grid of $\bf{k}$ points.


k mesh
 0
Gamma centered
 6 6 6
 0 0 0

POTCAR:
For GW calculations, one should use the POTCAR files with _GW appended to their name. These potentials were constructed to reproduce the atomic scattering properties over a larger energy range than the conventional potentials. This is needed when one needs to compute a large number of unoccupied states. You can check the TITEL field in your POTCAR file.


 PAW_PBE Li_GW 23Jun2021
  1.00000000000000
parameters from PSCTR are:
  VRHFIN =Li: s1p0
  LEXCH  = PE
  EATOM  =     5.3001 eV,    0.3895 Ry

  TITEL  = PAW_PBE Li_GW 23Jun2021
  LULTRA =        F    use ultrasoft PP ?
  IUNSCR =        1    unscreen: 0-lin 1-nonlin 2-no
  RPACOR =    1.500    partial core radius
  POMASS =    7.010; ZVAL   =    1.000    mass and valenz
  RCORE  =    2.100    outmost cutoff radius
  RWIGS  =    2.600; RWIGS  =    1.376    wigner-seitz radius (au A)
  ENMAX  =  112.104; ENMIN  =   84.078 eV

3.1.3 Run the calculation

Open a terminal, navigate to this example's directory and run VASP by entering the following:

cd $TUTORIALS/XAS/e03_*/dft
vasp_std

This calculation computes the DFT ground state. What is the size of the band gap? Is this a direct or indirect band gap system?
Hint: find the line fundamental gap: in the OUTCAR file.

Click here to reveal the answer!

The fundamental band gap is 6.4 eV

3.2.1 Task

In step two of this exercise, we compute the unoccupied KS orbitals and the frequency-dependent dielectric function in the independent particle picture. It is necessary to compute roughly 3 to 4 times as many KS orbitals as VASP would by default. In the previous DFT calculation, check the value of NBANDS written to the OUTCAR file! How many of these bands are occupied? Also, compare this to the value of NBANDS in the unoccupied-states/INCAR below!

3.2.2 Input files

unoccupied-states/INCAR:


Increase NBANDS with respect to the default and perform a single exact diagonalization of the Hamiltonian (ALGO=Exact and NELM=1). The LOPTICS tag is set so that the derivative of the Bloch states with respect to the Bloch wavevector $\bf{k}$ is calculated and written to the WAVEDER file. The WAVECAR and WAVEDER files of this calculation will be the starting point for the actual GW calculation.


SYSTEM  = LiCl  DFT unoccupied states

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

NBANDS  = 64     ! maybe even more bands 

ALGO    = Exact  ! use exact diagonalization of the Hamiltonian
NELM    = 1      ! since we are already converged stop after one step

LOPTICS = T      ! write WAVEDER

3.2.3 Run the calculation
cd $TUTORIALS/XAS/e03_*/unoccupied-states
cp ../dft/WAVECAR .
mpirun -np 2 vasp_std

As a byproduct, you are now able to plot the electronic dielectric function using py4vasp!

In [2]:
import py4vasp
optics_calc = py4vasp.Calculation.from_path("./e03_LiCl_G0W0/unoccupied-states")

optics_calc.dielectric_function.plot()

Does this dielectric function include the excitonic effects?

3.3.1 Task

In the third and final step, we perform the actual G$_0$W$_0$ calculation!

3.3.1 Input files

gw/INCAR:
To perform a single-shot GW calculation, specify ALGO=EVGW0 ("eigenvalue GW") and request a single electronic update (NELMGW=1). Do not forget to set NBANDS to the same value that was used to generate the WAVECAR and WAVEDER files in the previous step.


SYSTEM    = LiCl G0W0

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

NBANDS    = 64     ! maybe even more bands
NBANDSGW  = 32     ! calculate QP for 32 bands

! G0W0
ALGO      = EVGW0  ! use "GW0" for VASP.5.X
NELMGW    = 1      ! use "NELM" prior to VASP.6.3

NOMEGA    = 50     ! number of frequencies

3.3.2 Run the calculation
cd $TUTORIALS/XAS/e03_*/gw
cp ../unoccupied-states/{WAVECAR,WAVEDER} .
mpirun -np 4 vasp_std

You can find the quasiparticle energies QP-energies and the renormalization factor Z for each ${\bf k}$-point and band in the OUTCAR file. Search for the corresponding lines!

How large is the quasiparticle band gap in the G$_0$W$_0$ approximation? How does it compare with the experiment?

Hint: determine the band gap by inspecting the QP energies:

QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
for sc-GW calculations column KS-energies equals QP-energies in previous step
and V_xc(KS)=  KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1  + <V_x>^1)

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     -14.3214     -15.7668     -16.5974     -14.7351     -24.1387       0.7761       2.0000       1.6048
     2      -0.8608      -2.4378     -16.7439     -14.7633     -18.9286       0.7962       2.0000       0.0651
     3      -0.8608      -2.4378     -16.7439     -14.7633     -18.9286       0.7962       2.0000       0.0651
     4      -0.8608      -2.4378     -16.7439     -14.7633     -18.9286       0.7962       2.0000       0.0651
     5       5.5388       6.4198      -7.5780      -8.6292      -4.1817       0.8381       0.0000      -0.0690
     6      10.9956      12.3371      -7.9055      -9.5202      -4.2661       0.8308       0.0000      -0.0875
     7      10.9956      12.3371      -7.9055      -9.5202      -4.2661       0.8308       0.0000      -0.0875
     8      10.9956      12.3371      -7.9055      -9.5202      -4.2661       0.8308       0.0000      -0.0875
Click here to reveal the answer!

The fundamental QP band gap is 8.9 eV

Compare the dielectric function obtained within the independent-particle approximation (IPA) to the dielectric function in the random-phase approximation (RPA). The latter accounts for the fact that an electron that travels through a medium interacts (polarizes) its surroundings.

In [3]:
from py4vasp import Calculation

calc = Calculation.from_path("./e03_LiCl_G0W0/gw")
calc.dielectric_function.plot("Re(IPA, RPA)")

How does the dielectric function change from the IPA to the RPA? Look, for instance, to the change in the real part of the dielectric function at $\omega=0$ eV. Dielectric constant found in the experiment is 2.7. Which of the approximations is closer to the experimental value?

Similarly to the DFT calculations, the $GW$ calculations require a thorough convergence analysis. In addition to the usual DFT parameters, the $GW$ calculations require a convergence with the number of empty states NBANDS and the number of frequencies NOMEGA. In this tutorial, we have limited ourselves to a relatively low convergence.

3.4 Questions

  1. What quantity is approximated by $GW$?
  2. Where does VASP write the quasiparticle energies?
  3. At what level of approximation does VASP compute the screened Coulomb interaction in $GW$ per default?

4 X-ray absorption spectra in LiCl via BSE

By the end of this tutorial, you will be able to:

  • simulate X-ray absorption spectra with and without the excitonic effects
  • calculate and visualize the exciton charge density
  • distinguish dark and bright excitons

Short introduction

The Bethe-Salpeter equation (BSE) formalism can be used to account for the electron-hole interaction or the excitonic effects in the dielectric function. The Bethe-Salpeter equation can be formulated as a Dyson equation for the two-particle polarizability $L$ $$ \tag{2} L=L_0+L_0(v-W)L $$ Two terms describe the interaction: the bare Coulomb $V$ and the screened interaction $W$. The excitation energies $\omega_\lambda$ can be found by solving this equation as an eigenvalue problem. In the Tamm-Dancoff approximation (TDA) it reads

$$\tag{3} AX_\lambda=\omega_\lambda X_\lambda~, $$ where A is the two-particle Hamiltonian

$$\tag{4} A_{vc}^{v'c'} = (\varepsilon_v-\varepsilon_c)\delta_{vv'}\delta_{cc'} + \langle cv'|V|vc'\rangle - \langle cv'|W|c'v\rangle, $$ where $v,v'$ and $c,c'$ are indices of occupied and unoccupied states, correspondingly. In the case of the optical region, the occupied states $v,v'$ only include the valence states. However, in the case of X-ray absorption spectra (XAS), $v,v'$ have to include the excited core states as well. The energies and orbitals are typically taken from $GW$ quasiparticles.

In reciprocal space, the matrix $A$ can be written as

$$ \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{q})_\mathbf{G}\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 v\mathbf{k}|e^{i(\mathbf{q+G})}|v'\mathbf{k}'\rangle \langle v'\mathbf{k}'|e^{-i(\mathbf{q+G})}|v\mathbf{k}\rangle, \end{split} \tag{5} \end{equation} $$ where indices $v$ and $v'$ run over all occupied states, i.e., core and valence, $\Omega$ is the cell volume, $\bar{V}$ is the bare Coulomb potential without the long-range part

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

and the screened Coulomb potential

$$\tag{7} 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{8} \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.

4.1.1 Task

Compute the Li K-edge XAS spectrum of LiCl.

After completing the $G_0W_0$ calculation, we can solve the Bethe-Salpeter equation.

To perform a BSE calculation for the core states apart from the usual POTCAR, KPOINTS, and POSCAR files, which should be the same as in the preceding $G_0W_0$ step, we need the following files:

  • WAVECAR produced in the $G_0W_0$ step which contain the quasiparticle energies and orbitals
  • WFULLxxxx.tmp containing the RPA dielectric tensor, which was produced in the $G_0W_0$ calculation
  • WAVEDER* the orbitals derivatives are only required if valence orbitals are included

In the BSE calculation, we first construct the Hamiltonian, i.e., we compute all the matrix elements using the quasiparticles and the screened potential $W$. Then, the eigenvalues and the eigenvectors of Eq. (3) are found. Finally, the BSE dielectric function is calculated.

4.1.2 Input

To run a BSE calculation, we need to select the BSE algorithm ALGO = BSE. In contrast to the BSE calculation for the optical region, where we need to specify the number of valence and conduction bands, in the BSE calculation for the core states, only the conduction states need to be specified, while the number of valence states should be set to zero, i.e., NBANDSO = 0. To specify the core state of interest, we use the same tags as in the ΔSCF approach, i.e., ICORELEVEL, CLNT, CLN, CLL.

INCAR


SYSTEM     = LiCl BSE

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

! BSE
ALGO       = BSE   ! BSE algorithm
ANTIRES    = 0     ! Tamm-Dancoff approximation
NBANDS     = 64    ! same as in GW
NBANDSV    = 16    ! number of conduction bands in BSE
NBANDSO    = 0     ! no valence states for the XAS calculation
ICORELEVEL = 2     ! required for treating the core states in BSE
CLNT       = 1     ! species of the core levels (Li)
CLN        = 1     ! principal quantum number (n=1)
CLL        = 0     ! orbital quantum number (l=1 (s))
CSHIFT     = 0.25   ! complex shift/broadening

4.1.3 Calculation

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

cd $TUTORIALS/XAS/e04_LiCl_BSE/bse
cp $TUTORIALS/XAS/e03*/gw/{WAVECAR,WAVEDER,W*.tmp} .
mpirun -np 4 vasp_std

It can also be helpful to check the OUTCAR file to ensure that the calculations are performed as expected. A good point to check is the rank of the BSE Hamiltonian rank= 3456 in the OUTCAR file, which should correspond to the product of three numbers $N_c\times N_v\times N_k$, i.e., the number of cores states, number of conduction states, and k-points in the full Brillouin zone.

When the calculation is finished, we can analyze the results.

We can find the oscillator strengths in the vasprun.xml file.

<varray name="opticaltransitions" >
 <v>    56.9319     0.0000 </v>
 <v>    57.3676     1.3507 </v>
 <v>    57.3676     1.3507 </v>
 <v>    57.3677     1.3505 </v>
 <v>    57.9221     0.0000 </v>
 <v>    58.2199     0.0000 </v>
 <v>    58.3970     0.0887 </v>
 <v>    58.3970     0.0887 </v>
 <v>    58.3970     0.0887 </v>
 <v>    58.4844     0.0000 </v>

In this example, we see that the first transition is a dark exciton while the excitons at higher energies are bright, i.e., have non-zero oscillator strengths. The oscillator strengths are calculated in the dipole approximation with the corresponding selection rules. The transitions that violate these selection rules are not optically active (i.e., dark), and their strength is zero. Transitions that follow the selection rules are optically active (i.e., bright). Transitions 2, 3, and 4 are degenerate and have the same energy. The BSE dielectric function is also stored in vasprun.xml and can be found between lines:

 <dielectricfunction>
       ....
 </dielectricfunction>

Tip: If you get only zeros in the imaginary dielectric function, make sure that WAVECAR was read successfully, i.e., find these lines in stdout

the WAVECAR file was read successfully

Now, we can plot the BSE dielectric function and compare it to the experimental spectrum.

In [4]:
import py4vasp
import numpy as np
import os

E, eps = np.loadtxt(os.path.expanduser("./e04_LiCl_BSE/bse/expt.dat")).T
bse = py4vasp.Calculation.from_path("./e04_LiCl_BSE/bse/")
fig = bse.dielectric_function.to_plotly("Im(BSE)")
fig.add_scatter(x=E, y=eps, line_width=3, name="Expt.")
fig.data[0]['x']-=55.3
fig.data[0]['name']="BSE"
fig.update_xaxes(range=[-10, 30])
fig.update_layout(height=600,width=600)
fig.show()

The comparison of the calculated BSE spectrum shows a good agreement with the experimental spectrum taken from Shirley (2004) [J. Electron Spectrosc. Relat. Phenom. 579 (2004) 137-140]. The most prominent features are reproduced, and their relative positions are found to be in good agreement with the experiment. The absolute position of the XAS spectrum is not reproduced well in the BSE approach because the core states are frozen in the PAW approach and are taken directly from the POTCAR file. Furthermore, the $GW$ QP shifts are only applied to the valence and conduction states, while the core states are fixed at the DFT level.

It is important to note that the BSE spectrum is not well converged with the grid of $6\times6\times6$ and requires much denser sampling. The converged results can be found in Unzog et al. (2022) and Olovsson et al. (2009) [Phys. Rev. B (2022) 106 155133][Phys. Rev. B 79 (2009) 041102(R)].

Furthermore, the dielectric function has to be thoroughly converged with respect to the number of conduction states. The excitons are strongly localized in XAS, and a large number of conduction states might be required to achieve a properly converged spectrum.

The comparison of the spectra calculated with the $\Delta$SCF and BSE approaches can be found in Unzog et al. (2022) [Phys. Rev. B (2022) 106 155133]. As shown in the paper, the two approaches agree well, and the BSE spectrum is overall closer to the experiment.

4.1.4 Questions

  1. What is the difference between bright and dark excitons?
  2. Why do we need to shift the calculated spectrum in the absolute scale?

4.2 Task

Compute the Li K-edge XAS spectrum of LiCl in the RPA approximation.

The BSE algorithm in VASP allows us to calculate the dielectric function in various approximations. To emphasize the importance of the excitonic effects in XAS, we will calculate the dielectric function in RPA and compare the results to the BSE dielectric function.

4.2.2 Input

To select the RPA approximation we need to disable the ladder diagrams or the screened Coulomb interaction in the two-particle Hamiltonian LADDER=.FALSE., but we will keep the bare Coulomb interaction LHARTEE=.TRUE..

INCAR


SYSTEM     = LiCl BSE

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

! BSE
ALGO       = BSE     ! BSE algorithm
LADDER     = .FALSE. ! RPA approximation
LHARTREE   = .TRUE.  ! 
ANTIRES    = 0       ! Tamm-Dancoff approximation
NBANDS     = 64      ! same as in GW
NBANDSV    = 16      ! number of conduction bands in BSE
NBANDSO    = 0       ! no valence states for the XAS calculation
ICORELEVEL = 2       ! required for treating the core states in BSE
CLNT       = 1       ! species of the core levels (Li)
CLN        = 1       ! principal quantum number (n=1)
CLL        = 0       ! orbital quantum number (l=1 (s))
CSHIFT     = 0.25    ! complex shift/broadening

4.2.3 Calculation

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

cd $TUTORIALS/XAS/e04_LiCl_BSE/rpa
cp $TUTORIALS/XAS/e03*/gw/{WAVECAR,WAVEDER,W*.tmp} .
mpirun -np 4 vasp_std

When the calculation is complete, we can plot the dielectric functions and compare the BSE and RPA dielectric functions for the core-state excitations.

In [5]:
import py4vasp

rpa = py4vasp.Calculation.from_path("./e04_LiCl_BSE/rpa")
bse = py4vasp.Calculation.from_path("./e04_LiCl_BSE/bse")
(
    rpa.dielectric_function.plot("BSE(Im)").label("RPA") +
    bse.dielectric_function.plot("BSE(Im)").label("BSE")
)

The only difference between these dielectric functions is the excitonic effects. Hence, this comparison clearly shows the importance of the interaction between the core-hole and electrons for XAS. The binding energy of the first bright exciton can be estimated as the energy difference between the eigenvalue found in the RPA and BSE approximations.

Click here to reveal the answer!

The binding energy of the first bright exciton is 59.08-57.37=1.71 eV

Plot the real part of the calculated dielectric function.

In [6]:
import py4vasp

rpa = py4vasp.Calculation.from_path("./e04_LiCl_BSE/rpa")
bse = py4vasp.Calculation.from_path("./e04_LiCl_BSE/bse")
(
    rpa.dielectric_function.plot("BSE(Re)").label("RPA") +
    bse.dielectric_function.plot("BSE(Re)").label("BSE")
)

If we compare the real part of the RPA dielectric function found in the BSE algorithm to the RPA results from the $G_0W_0$ step, we will see that they differ dramatically. This difference is primarily because only the core electrons were considered in the BSE algorithm. The RPA dielectric function from the BSE algorithm is close to 1, which shows that the 1s core states contribute very little to the screening in this system.

Another aspect that distinguishes the RPA dielectric function in the $G_0W_0$ calculation from the one obtained in this example is the Tamm-Dancoff approximation (TDA), which is not numerically significant for this case but important from the conceptual point of view. The RPA dielectric function in the BSE algorithm was calculated in TDA, while this approximation was not used in the $G_0W_0$ calculation.

4.2.4 Questions

  1. What is the difference between the RPA and BSE dielectric functions?
  2. How can the exciton binding energy be estimated?
  3. Why is the RPA dielectric constant found via the BSE algorithm much smaller than the one found in the $G_0W_0$ step?

4.3.1 Task

Compute the exciton wavefunction for the Li K-edge of LiCl.

Visualizing the exciton wavefunction or the corresponding charge density can be useful for analyzing the symmetry and the localization of the exciton in the system.

The exciton wavefunction is written as a function of coordinates of two particles (one hole and one electron) $\psi_\lambda(\mathbf{r}_e,\mathbf{r}_h)=\sum_{vc} X_{vc}^\lambda \phi_v^*(\mathbf{r}_h)\phi_c(\mathbf{r}_e)$. To visualize such a function in 3D space, we need to "fix" one of the coordinates. Considering that in XAS, the core hole is strongly localized, it is common to fix the coordinates of the hole.

4.3.2 Input

To calculate the excitonic wavefunction, we must select the number of eigenvectors we would like to analyze NBSEEIG and fix the core-hole position BSEHOLE to the coordinates of the excited atom.

INCAR


SYSTEM     = LiCl BSE

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

! BSE
ALGO       = BSE     ! BSE algorithm
ANTIRES    = 0       ! Tamm-Dancoff approximation
NBANDSV    = 16      ! number of conduction bands in BSE
NBANDSO    = 0       ! no valence states for the XAS calculation
ICORELEVEL = 2       ! required for treating the core states in BSE
CLNT       = 1       ! species of the core levels (Li)
CLN        = 1       ! principal quantum number (n=1)
CLL        = 0       ! orbital quantum number (l=1 (s))
CSHIFT     = 0.1     ! complex shift/broadening

NBSEEIG    = 3       ! six lowest excitons
BSEHOLE    = 0.0 0.0 0.0 ! fix the core-hole position
LCHARGH5   = .TRUE.  ! to write exciton charge density to vaspout.h5
PRECFOCK   = low     ! smaller FFT grid to reduce the exciton charge density files size

4.3.3 Calculation

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

cd $TUTORIALS/XAS/e04_LiCl_BSE/wf
cp $TUTORIALS/XAS/e03*/gw/{WAVECAR,WAVEDER,W*.tmp} .
mpirun -np 4 vasp_std

When the calculation is complete, we can plot the charge density of the first bright exciton. The first bright exciton is triply-degenerate, and to plot the corresponding wavefunction, we need to sum up all its three components.

In [7]:
import py4vasp

bse = py4vasp.Calculation.from_path("./e04_LiCl_BSE/wf/")

view = bse.exciton.density.to_ngl("2+3+4",isolevel=1,center=True)
c = view.trajectory_0
c.update_ball_and_stick(opacity=0.5)
view

By visualizing the exciton charge density we can see that it is well localized between the two nearest Li atoms and is centered around a Li atom.

By changing the first argument of the method to_ngl we can select the first dark exciton bse.exciton.density.to_ngl("1", isolevel=1, center=True).

4.3.4 Questions

  1. At what atom is the exciton localized?
  2. How localized is the first bright exciton?

Well done! You have finished Part 2 of the X-ray absorption spectra tutorial!