Part 3: Bethe-Salpeter equation calculations


9 G0W0 electronic structure

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

  • Run a G0W0 calculation starting with the PBE+U wavefunctions.
  • Calculate the unoccupied states Kohn-Sham orbitals and the optical absorption in the independent particle approximation (IPA).
  • Determine the G0W0 gap of NiO.
  • Calculate optical absorption in the random-phase approximation (RPA).

9.1 Task

Calculation of G0W0 quasiparticle energies and screened interaction using a 3x3x3 k-mesh for antiferromagnetic NiO using PBE+U (U=6.181, J=1.134)

An alternative approach to modeling strongly correlated systems, is to use Many-Body Perturbation Theory (MBPT) to correct the band gap and simulate optical properties of AFM NiO.

As a first step, we need to find the solutions of the Hedin's equation in the $GW$ approximation, i.e., quasiparticles that include the polarization of charge induced around them as they propagate through the medium. The screening in the $GW$ self-energy is obtained by the random-phase approximation (RPA) for the polarizability function $\chi$:

$$\tag{1} \chi(1,2)=-i\int d(3,4) G(1,3) G(4,1) \Gamma(3,4 ; 2), $$ where the vertex function is $\Gamma(3,4;2)=\delta(3-2)\delta(4-2)$ and $G$ are Green's functions.

The RPA approximation that uses the PBE electronic structure as a starting point for the description of screening relies on partial error cancellation. Namely, the band gap underestimated in the PBE approach leads to an overestimation of the screening, while the lack of vertex corrections causes an underestimation of the screening. Thus, in semiconductors where the bands are dominated by delocalized $s$ and $p$ states, this error cancellation has been shown to work well and yields a good description of the screening [Phys. Rev. Lett. 99 (2007) 246403].

However, in systems where the bands are dominated by localized $d$ states, this error cancellation is often violated. In the case of antiferromagnetic (AFM) NiO, where the PBE band gap is much more severely underestimated compared to experiment, so the RPA approximation leads to a strong overestimation of the screening. In particular, the dielectric constant is reported to be about 15, while the experimental value is approximately 5.7 [Phys. Rev. Materials 2 (2018) 073803]. This issue becomes even more severe in Mott-insulating materials, where PBE yields a metallic state.

There are several approaches that can remedy the error in the description of screening. The more sophisticated approach requires solving Hedin’s equations self-consistently and including vertex corrections in the polarizability function, as in the $QSG\tilde{W}$ approximation [Phys. Rev. B 108 (2023) 165104, Phys. Rev. Lett. 99 (2007) 246403]. However, this method is computationally very demanding and lies beyond the scope of this tutorial. A much cheaper alternative is to start from an improved electronic structure compared to PBE. In this tutorial, we use the PBE+U approximation, where the Hubbard $U$ is determined from preceding cRPA calculations (see part 1 of these strongly correlated tutorials).

In this exercise, we calculate the G0W0 quasiparticle energies and the screened interaction W, starting from the PBE+U wavefunctions. In general, the convergence of the PBE+U calculation can be improved by restarting from a preceding standard PBE calculation. In this case we omit this preliminary step and start with the PBE+U calculation.

This calculations consists of three steps:

  • Calculate the PBE+U ground-state electronic structure of AFM NiO.
  • Include a large number of empty states and compute the orbital derivatives.
  • Perform the G0W0 calculation, starting from the electronic structure obtained using the PBE+U approach.

In this calculation, you will calculate the G0W0 quasiparticle energies for antiferromagnetic (AFM) NiO.

9.2 PBE+U ground state calculation

9.2.1 Input files

The input files to run the DFT ground state calculation are prepared at $TUTORIALS/strongly_correlated/e09_G0W0/dft. Check them out!


POSCAR:
The unit cell of AFM NiO has been increased to 4 atoms to be able to describe the antiferromagnetic ordering of the moments.


NiO
4.171
   1.0 0.5 0.5
   0.5 1.0 0.5
   0.5 0.5 1.0
   Ni    O
    2    2
Direct
     0.000000000         0.000000000         0.000000000
     0.500000000         0.500000000         0.500000000
     0.250000000         0.250000000         0.250000000
     0.750000000         0.750000000         0.750000000

dft/INCAR:
A DFT run using the PBE+U approach. Below we emphasize some of the tags specific for this type of calculations:

  • ISPIN turns on spin-polarized calculations.
  • LASPH should be considered a standard setting, especially for magnetic calculations. It includes aspherical contribution to the gradient corrections inside the PAW spheres.
  • MAGMOM sets the magnetic moments for each specific cite, i.e., AFM in our case.
  • As we wish to restart the calculation, it is important to set LMAXMIX to 4 for $d$ electrons (and 6 for $f$ electrons). Otherwise the charge density written to the CHGCAR file is incomplete. But in any case, it might speed up convergence as the charge-density mixer consideres states up to LMAXMIX.
  • LDAU switch the DFT+U method.
  • LDAUTYPE select Dudarev's approach.
  • LDAUL specifies the l-quantum number for which the on-site interaction is added.
  • LDAUU sets the effective on-site Coulomb interactions.
  • LDAUJ Sets the effective on-site exchange interactions.

SYSTEM   = NiO
ISMEAR   = 0
SIGMA    = 0.05
ISPIN    = 2
MAGMOM   = 2.0 -2.0 2*0
LDAU      = .TRUE.
LDAUTYPE  = 2
LDAUL     = 2 -1
LDAUU     = 6.181 0.00
LDAUJ     = 1.134 0.00
LMAXMIX   = 4 
LASPH    = T
EDIFF    = 1E-8

KPOINTS:
An equally spaced (3$\times$3$\times$3) Γ-centered grid of $\bf{k}$ points. This Brillouin zone (BZ) sampling allows us to obtain the results for the purpose of this tutorial. A proper convergence analysis w.r.t. the k-point sampling is required for a fully converged calculation.


Automatically generated mesh
  0
Gamma
 3 3 3

POTCAR:
For G0W0 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. For improved results, potentials with the semicore s- and p- states treated in the valence (_sv and _pv), can be required to produce accurate results.


   PAW_PBE Ni_GW 31Mar2010
   PAW_PBE O_GW_new 19Mar2012             

9.2.2 Run the calculation

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

cd $TUTORIALS/strongly_correlated/e09_G0W0/dft
mpirun -np 4 vasp_std

Find the band gap in the OUTCAR file by searching the line fundamental gap or use py4vasp to determine the gap.

In [2]:
import py4vasp
dft = py4vasp.Calculation.from_path("./e09_G0W0/dft")
dft.bandgap.print()
Band structure
--------------
                       spin independent             spin component 1             spin component 2
val. band max:               3.826558                     3.826558                     3.826558
cond. band min:              7.392539                     7.392539                     7.392539
fundamental gap:             3.565981                     3.565981                     3.565981
VBM @ kpoint:       0.3333   0.3333   0.3333     0.3333   0.3333   0.3333     0.3333   0.3333   0.3333
CBM @ kpoint:      -0.3333   0.3333   0.3333    -0.3333   0.3333   0.3333    -0.3333   0.3333   0.3333

lower band:                  3.826558                     3.826558                     3.826558
upper band:                  7.688477                     7.688477                     7.688477
direct gap:                  3.861919                     3.861919                     3.861919
@ kpoint:           0.3333   0.3333   0.3333     0.3333   0.3333   0.3333     0.3333   0.3333   0.3333

Fermi energy:                4.057135

The fundamental bandgap of AFM NiO determined experimentally is 4.3 eV [Phys. Rev. Lett. 53 (1984) 2339]. The PBE exchange–correlation functional strongly underestimates the bandgap, yielding a value of about 0.95 eV (try turning off the U to check or see Ref. [J. Phys. Chem. A 121 (2017) 3318]). Adding the Hubbard correction significantly improves the bandgap to ~3.6 eV, although it still strongly underestimates it.

With this calculation, you have successfully obtained the WAVECAR file necessary for generating the virtual orbitals used in the GW approximation.

9.3 Compute additional unoccupied Kohn-Sham orbitals

In step two of this exercise, we obtain a large number of unoccupied states and compute the derivatives of the orbitals w.r.t. the $\textbf{k}$-points (stored in WAVEDER) that are required for the G0W0 calculation.

9.3.1 Input files

The input files are prepared in $TUTORIALS/strongly_correlated/e09_G0W0/unoccupied/. POSCAR, POTCAR, and KPOINTS are the same as in the dft directory. Check them out!


unoccupied/INCAR:
We use the INCAR file from the previous step and add the tags to compute a larger number of empty states and compute the orbital derivatives:

  • ALGO sets the electronic minimization algorithm; Exact performs an exact diagonalization of the Kohn-Sham (KS) Hamiltonian, which is a preferred way to compute a large number of empty bands in small cells.
  • LOPTICS calculates the frequency-dependent dielectric matrix after the electronic ground state has been determined and computes the orbital derivatives (WAVEDER)
  • NBANDS defines the number of bands, i.e., all the occupied and conduction bands of interest. To obtain accurate QP energies in GW, the number of bands needs to be increased.

SYSTEM    = NiO
ISMEAR    = 0
SIGMA     = 0.05
ISPIN     = 2
MAGMOM    = 2.0 -2.0 2*0
LDAU      = .TRUE.
LDAUTYPE  = 2
LDAUL     = 2 -1
LDAUU     = 6.181 0.00
LDAUJ     = 1.134 0.00
LMAXMIX   = 4 
LASPH     = T

EDIFF    = 1E-8

ALGO     = Exact
NBANDS   = 256
LOPTICS  = .TRUE.

9.3.2 Run the calculation

To calculate the unoccupied states, go to the unoccupied directory and run the calculation:

cd $TUTORIALS/strongly_correlated/e09_G0W0/unoccupied
cp ../dft/WAVECAR .
mpirun -np 4 vasp_std

At the end of this calculation, you should find both the WAVECAR and WAVEDER files in the unoccupied directory.

LOPTICS=.TRUE. computes the orbital derivatives and the dielectric function in the independent particle approximation (IPA). You can analyze the dielectric function by running the next cell.

In [3]:
import py4vasp
unoccupied=py4vasp.Calculation.from_path("./e09_G0W0/unoccupied/")
unoccupied.dielectric_function.plot()

By hovering the cursor over the plot you can find the value of the real part of $\varepsilon$ at $\omega=0$, i.e., the dielectric constant. You can see that it is overestimated compared to the experimental value of 5.7. This is the consequence of the underestimated band gap. In cases where the gap is well reproduced, the dielectric constant in IP approximation is underestimated due to the lack of contributions from the electron-hole interaction.

9.4 Compute the G0W0 quasiparticles and the screened potential

In step three of this exercise, we will calculate the G0W0 quasiparticle energies of AFM NiO as well as the screened potential $W$ using the WAVECAR and WAVEDER from the previous step.

9.4.1 Input files

The input files are prepared in $TUTORIALS/strongly_correlated/e09_G0W0/g0w0/. POSCAR, POTCAR, and KPOINTS are the same as in the dft directory.


g0w0/INCAR:
Use the same number of bands as in the previous step where the unoccupied orbitals and their derivatives were calculated. Notice here that the Hubbard correction should not be applied in the G0W0 calculation.

  • NOMEGA specify the number of real frequencies in the polarizability and the self-energy.
  • ALGO set the G0W0 algorithm (EVGW0) and NELM=1.
  • PRECFOCK and ENCUTGW the values are reduced relative to the default settings for the purpose of this tutorial. A proper convergence analysis is necessary for all convergence parameters.

SYSTEM   = NiO
ISPIN    = 2
MAGMOM   = 2.0 -2.0 2*0

EDIFF    = 1E-8

ISMEAR   = 0
SIGMA    = 0.05

ALGO     = EVGW0
NELM     = 1
NBANDS   = 128
NOMEGA   = 50
PRECFOCK = low
ENCUTGW = 200

9.4.2 Run the calculation

To perform the GW approximation, go to the g0w0 directory and run the calculation:

cd $TUTORIALS/strongly_correlated/e09_G0W0/g0w0
cp ../unoccupied/{WAVECAR,WAVEDER} .
mpirun -np 4 vasp_std

The quasiparticle energies are written to the OUTCAR file. Find their values in the column QP-energies after the line spin component 1 and determine the G0W0 bandgap.

You should find a value of 4.6 eV, which is reasonably close to the experimental value of 4.3 eV. However, it is important to keep in mind that our calculation is not fully converged with respect to the $\textbf{k}$-point mesh, the number of frequencies, and the number of bands. Therefore, a proper convergence analysis is necessary.

Next, plot the dielectric function and determine the value of the RPA dielectric constant.

In [4]:
import py4vasp
gw = py4vasp.Calculation.from_path("./e09_G0W0/g0w0")
gw.dielectric_function.plot("RPA(Re,Im)")

The RPA dielectric constant is 5.76, which in this case is higher then the experimental value.

However, if we perform the same calculations with a more stringent convergence criteria, i.e., $k$-points, NBANDS, NOMEGA, we will find the RPA dielectric constant of 5.2 and the band gap of 4.25 eV. The dielectric constant is underestimated in RPA compared to experiment mostly due to the lack of the excitonic effects. Thus, we conclude that by using the Hubbard correction in the PBE+U approach, we are able to improve the electronic structure compared to standard PBE. G0W0 starting from PBE+U yields both a good description of the electronic structure, i.e., the fundamental band gap, and the screened interaction W, which is required for the subsequent BSE calculations.

9.5 Questions

  1. What important effects are missing in the RPA dielectric function?
  2. What approximation for the dielectric function is used for LOPTICS=.TRUE.?

10 Optical properties of NiO

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

  • Calculated optical absorption of NiO using BSE.
  • Calculate electron-energy loss spectra (EELS) of NiO.
  • Find the spin-flip excitations in NiO.

10.1 Optical absorption spectrum via BSE

Calculate optical absorption spectra via the Bethe-Salpeter equation (BSE) for antiferromagnetic NiO.

Within MBPT, the electron-hole interaction or excitonic effects in the spectra can be accounted for via the Bethe-Salpeter equation (BSE): $$\tag{2} \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{3} A_{vc}^{v'c'} = (\varepsilon_v^{QP}-\varepsilon_c^{QP})\delta_{vv'}\delta_{cc'} + \langle cv'|V|vc'\rangle - \langle cv'|W|c'v\rangle $$

and $B$

$$\tag{4} B_{vc}^{v'c'} = \langle vv'|V|cc'\rangle - \langle vv'|W|c'c\rangle. $$

describe the transitions between occupied $v,v'$ and unoccupied $c,c'$ states coupled via bare Coulomb interaction $V$ and the screened potential $W$. The solutions of this equation $\omega_\lambda$ and $\left(X_\lambda,Y_\lambda\right)$ are the excitonic states.

In the Tamm-Dancoff approximation, the coupling between resonant and anti-resonant parts is neglected, i.e., $B=0$. Hence, the BSE equation reduces to $$\tag{5} A(\mathbf{q})|X_\lambda\rangle=E_\lambda|X_\lambda\rangle, $$ where in the general case of finite momentum $\mathbf{q}$ the Hamiltonian can be written as $$\tag{6} A_{\mathbf{k}cv}^{\mathbf{k}'c'v'}(\mathbf{q})= \delta_{\mathbf{k}\mathbf{k'}} \delta_{v v'} \delta_{c c'}\left(\varepsilon_{\mathbf{k}+\mathbf{q} v}^{QP}-\varepsilon_{\mathbf{k} c}^{QP}\right) -\left(f_{\mathbf{k}+\mathbf{q} v}-f_{\mathbf{k} c}\right) K_{\mathbf{k} c, \mathbf{k}+\mathbf{q} v, \mathbf{k'}+\mathbf{q} v', \mathbf{k'} c'}. $$ The kernel $K$ contains the bare Coulomb potental $V$ and the static screened potential $W$. In the collinear case, the spin-structure of the kernel can be written in the basis $|\uparrow\uparrow\rangle$,$|\uparrow\downarrow\rangle$,$|\downarrow\uparrow\rangle$,$|\downarrow\downarrow\rangle$ $$\tag{7} K=\left[\begin{array}{cccc} V-W & 0 & 0 & V \\ 0 & 0 & -W & 0 \\ 0 & -W & 0 & 0 \\ V & 0 & & V-W \end{array}\right]. $$ Thus, the BSE Hamiltonian decouples into two classes of solutions: one containing the bare Coulomb potential $V$, corresponding to spin-conserving transitions, and another without it, corresponding to spin-flip transitions. The spin-conserving excitations can be optically active (bright) or forbidden by the dipole selection rules (dark). The spin-flip transitions in the collinear case are always dark. When spin–orbit coupling is included and the spins become noncollinear, the decomposition above is no longer possible, and the full kernel $K$ must be solved.

In this exercise, you will solve the Bethe-Salpeter equation (BSE) in the Tamm-Dancoff approximation for the spin-conserving solutions [Phys. Rev. B 92 (2015) 045209].

10.1.1 Input files

The input files to run this example are prepared at $TUTORIALS/strongly_correlated/e10_BSE/absorption.

We use the same POSCAR, KPOINTS, and POTCAR files as in the preceding G0W0 step (e09). In the INCAR file, we select the BSE algorithm and specify the number of valence and conduction band states that will be included in the BSE kernel.

There are three different algorithms for solving the Bethe–Salpeter equation. The most accurate approach is exact diagonalization (IBSE = 2), which we use in this calculation. This method allows us to directly obtain both the eigenvalues and eigenvectors, which are required for computing the oscillator strengths. However, this algorithm scales as N3, where N is the size of the BSE Hamiltonian matrix.


SYSTEM   = NiO
ISPIN    = 2
MAGMOM   = 2.0 -2.0 2*0

EDIFF    = 1E-8

ISMEAR   = 0
SIGMA    = 0.05

NBANDS   = 128
ALGO     = BSE
IBSE     = 2
NBANDSO  = 16 ! number of occupied states treated in BSE
NBANDSV  = 16 ! number of unoccupied states treated in BSE
OMEGAMAX = 16
PRECFOCK = low
ENCUTGW = 200
CSHIFT = 0.3

10.1.2 Run the calculation

Navigate to the e10_BSE/absorption directory and run the calculation:

cd $TUTORIALS/strongly_correlated/e10_BSE/absorption
cp ../../e09_G0W0/g0w0/W* .
mpirun -np 4 vasp_std
In [5]:
import py4vasp
import numpy as np
import os
bse = py4vasp.Calculation.from_path("./e10_BSE/absorption")
E_c, ε_im_c, e_im_c  = np.loadtxt(os.path.expanduser("./e10_BSE/absorption/optics_converged.dat")).T
E, ε = np.loadtxt(os.path.expanduser("./e10_BSE/absorption/expt.dat")).T
(
    py4vasp.plot(E, ε, "Expt") +
    py4vasp.plot(E_c, ε_im_c, "converged") +
    bse.dielectric_function.plot("BSE(Im)").label("BSE")
)

The calculated spectrum is in reasonably good agreement with experiment, i.e., the position and the overall shape of the spectrum are reproduced. The spectrum is truncated at 16 eV, which is the result of the OMEGAMAX limit that we set in our calculation. For reference, we provide the results for a well-converged BSE calculation to show that by sampling $k$-points in 3x3x3 grid we are able to correctly reproduce the position of the main peak and overall the shape of the spectrum.

The solutions of the BSE equation in the static approximation are real-valued and have infinite lifetime. However, in reality, the excitons have a finite lifetime due to coupling to other excitons, with phonons being the most dominant contribution. These effects are completely neglected in our approach, and hence to account for the lifetime broadening, we introduce a small complex shift in the eigenvalues, which gives a Lorentzian broadening. This shift is set by CSHIFT and is usually adjusted to match the width of the prominent features to experiment. In this case, we used the value of 0.3 eV.

To find the oscillator strengths, open the vasprun.xml file and find the block after the line <varray name="opticaltransitions" >. You can do this by entering the following line into the terminal:

grep -A 16 '<varray name="opticaltransitions" >' vasprun.xml
In [6]:
%%bash
cd e10_BSE/absorption
grep -A 16 '<varray name="opticaltransitions" >' vasprun.xml
cd ..
  <varray name="opticaltransitions" >
   <v>     1.6086     0.0000 </v>
   <v>     1.6260     0.0000 </v>
   <v>     1.6281     0.0000 </v>
   <v>     1.6281     0.0000 </v>
   <v>     1.6355     0.0000 </v>
   <v>     1.6355     0.0000 </v>
   <v>     2.6415     0.0000 </v>
   <v>     2.6526     0.0000 </v>
   <v>     2.6641     0.0000 </v>
   <v>     2.6641     0.0000 </v>
   <v>     2.6658     0.0000 </v>
   <v>     2.6658     0.0000 </v>
   <v>     3.9272    18.5672 </v>
   <v>     4.0086    11.6043 </v>
   <v>     4.0086    11.6043 </v>
   <v>     4.0461     0.0000 </v>

The first column is the energy of the excitonic state, and the second column is the oscillator strength. The first 12 transitions have the strength of zero. These are the transitions that violate the selection rules and are called "dark" excitons. The first "bright" exciton appears at around 4 eV, which is also visible in the absorption spectrum.

10.1.3 Questions

  1. Why do some of the transitions have zero strength?
  2. Can spin-flip excitations be bright in the case of collinear spins?
  3. How do we set the Lorentzian broadening?

10.2 EELS via BSE

Calculate electron energy loss spectrum (EELS) via the Bethe-Salpeter equation (BSE) for antiferromagnetic NiO.

In electron-energy loss spectroscopy (EELS), the incident electrons carry a finite momentum and can therefore excite states not only vertically, as in the case of photons, but also at a finite momentum transfer $q$. The value of this momentum transfer has to be provided via the KPOINT_BSE tag. The spectrum can only be calculated for the specific vector $q$ included in the calculation. The index of the corresponding $k$-point can be found in the OUTCAR file after the line irreducible k-points.

The EELS signal can be obtained from the macroscopic dielectric function calculated within the BSE: $$\tag{8} \operatorname{EELS}(\mathbf{q}, \omega)=-\operatorname{Im}\left[\epsilon^{-1}(\mathbf{q}, \omega)\right] $$

For optical spectra, the Tamm–Dancoff approximation has been shown to work well. However, finite-momentum calculations must include the coupling between resonant and anti-resonant transitions (ANTIRES=2) [Phys. Rev. B 91, 045209 (2015)].

To observe plasmon excitations in the EEL spectrum, we need to compute the spectrum over a wide energy range (OMEGAMAX=45). This significantly increases the size of the BSE matrix, making exact diagonalization too expensive for the purposes of this tutorial. Instead, we can use one of the iterative algorithms: IBSE=1 (time-evolution algorithm) or IBSE=3 (Lanczos). In principle, the Lanczos algorithm is usually the fastest way to compute the dielectric function. However, for finite-momentum calculations, we need to solve the full BSE (ANTIRES=2). However, the current implementation of the Lanczos algorithm does not support calculations for full BSE and we will use the time-evolution algorithm instead. Notice, that in the case of IBSE=1, only the dielectric function is written out to vasprun.xml, while the <varray name="opticaltransitions" > block is missing.

10.2.1 Input files

The input files to run this example are prepared at $TUTORIALS/strongly_correlated/e10_BSE/eels.


By analyzing the OUTCAR from the previous step we can find the list of $k$-points in the irreducible BZ:

Found     13 irreducible k-points:

Following reciprocal coordinates:
           Coordinates               Weight
    0.000000  0.000000  0.000000      1.000000
--> 0.250000  0.000000  0.000000      6.000000 <--
    0.500000  0.000000  0.000000      3.000000
...

In experiments, spectra are typically measured at small momentum transfer, and the calculation may therefore require a dense $k$-point sampling in order to include the corresponding $q$-point. In this case, we select the second $k$-point (0.250000 0.000000 0.000000) with KPOINT_BSE = 2.


INCAR


SYSTEM   = NiO
ISPIN    = 2
MAGMOM   = 2.0 -2.0 2*0

EDIFF    = 1E-8

ISMEAR   = 0
SIGMA    = 0.05

NBANDS   = 128
ALGO     = BSE
ANTIRES  = 2
IBSE     = 1
NBANDSO  = 12
NBANDSV  = 12
PRECFOCK = low
ENCUTGW = 200
OMEGAMAX = 45
CSHIFT = 0.3
KPOINT_BSE = 2

10.2.2 Run the calculation

Navigate to the e10_BSE/eels directory and run the calculation:

cd $TUTORIALS/strongly_correlated/e10_BSE/eels
cp ../../e09_G0W0/g0w0/W* .
mpirun -np 4 vasp_std

Run the following command to extract the dielectric function and plot the EEL spectrum.

In [7]:
import py4vasp
import numpy as np
import os
eels = py4vasp.Calculation.from_path("./e10_BSE/eels")
data = eels.dielectric_function.read("BSE")
E, ε = np.loadtxt(os.path.expanduser("./e10_BSE/eels/expt.dat")).T
E_c, ε_im_c, ε_re_c  = np.loadtxt(os.path.expanduser("./e10_BSE/eels/optics_converged.dat")).T

E_eps=data['energies']
eps=data['dielectric_function']
(
    py4vasp.plot(E, ε, "Expt", xlabel="Energy (eV)", ylabel="dielectric function ε") +
    py4vasp.plot(E_c, ε_im_c/(ε_im_c**2+ε_re_c**2), "converged") +
    py4vasp.plot(E_eps, eps.imag/((eps.imag)**2+(eps.real)**2), "EELS")
)

The calculated EELS spectrum is in a reasonable agreement with experiment. The most prominent feature in the spectrum starts with a steep rise of the loss function starting at 20 eV, which corresponds to longitudinal plasma frequency or the collective charge oscillations of the valence-electron.

10.2.3 Questions

  1. Is Tamm-Dancoff a good approximation for finite $q$ calculations?
  2. What type of excitations is the peak at 20 eV associated with?

10.3 Spin-flip excitations

Calculate the spin-flip excitations via the Bethe-Salpeter equation (BSE) for antiferromagnetic NiO.

Perform a BSE calculation with LHARTREE=.FALSE. to find the spin-flip solutions.

10.3.1 Input files

The input files to run this example are prepared at $TUTORIALS/strongly_correlated/e10_BSE/spin_flip.


INCAR


SYSTEM   = NiO
ISPIN    = 2
MAGMOM   = 2.0 -2.0 2*0

EDIFF    = 1E-8

ISMEAR   = 0
SIGMA    = 0.05

NBANDS   = 128
ALGO     = BSE
IBSE     = 2
NBANDSO  = 16
NBANDSV  = 16
PRECFOCK = low
ENCUTGW = 200
LHARTREE = .FALSE. ! for spin-flip transitions

10.3.2 Run the calculation

Navigate to the e10_BSE/spin_flip directory and run the calculation:

cd $TUTORIALS/strongly_correlated/e10_BSE/spin_flip
cp ../../e09_G0W0/g0w0/W* .
mpirun -np 4 vasp_std

Open the vasprun.xml file and find the block after the line :

grep -A 16 '<varray name="opticaltransitions" >' $TUTORIALS/strongly_correlated/e10_BSE/spin_flip/vasprun.xml
In [8]:
%%bash
cd e10_BSE/spin_flip
grep -A 16 '<varray name="opticaltransitions" >' vasprun.xml
cd ..
  <varray name="opticaltransitions" >
   <v>     1.5514     0.0000 </v>
   <v>     1.5514     0.0000 </v>
   <v>     1.5714     0.0000 </v>
   <v>     1.5714     0.0000 </v>
   <v>     1.5714     0.0000 </v>
   <v>     1.5714     0.0000 </v>
   <v>     2.2156     0.0000 </v>
   <v>     2.2156     0.0000 </v>
   <v>     2.2427     0.0000 </v>
   <v>     2.2427     0.0000 </v>
   <v>     2.2427     0.0000 </v>
   <v>     2.2427     0.0000 </v>
   <v>     3.8905    10.1435 </v>
   <v>     3.8905    10.1435 </v>
   <v>     3.9148     0.0000 </v>
   <v>     3.9148     0.0000 </v>

Note that the energy of the first transition is 50 meV lower than in the absorption calculation. This shift arises from the absence of the repulsive bare Coulomb contribution (LHARTREE=.FALSE.). The resulting transitions correspond to spin-flip excitations and therefore violate spin conservation for optical dipole transitions, making them optically inactive. These excitations include the collective spin excitations, i.e., magnons, which can be analyzed within BSE.

10.3.3 Questions

  1. Are spin-flip transitions optically active?
  2. What term in the BSE Hamiltonian is turning to zero for collinear spin-flip transitions?

Well done! You have finished Part 3 of the strongly correlated tutorial!