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.
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:
- The first step, is a DFT-ground-state calculation.
- 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).
- 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!
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
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.
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¶
- What quantity is approximated by $GW$?
- Where does VASP write the quasiparticle energies?
- 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.
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.
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¶
- What is the difference between bright and dark excitons?
- 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..
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.
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.
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¶
- What is the difference between the RPA and BSE dielectric functions?
- How can the exciton binding energy be estimated?
- 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.
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.
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¶
- At what atom is the exciton localized?
- How localized is the first bright exciton?