Bethe-Salpeter-equations calculations

VASP offers a powerful module for solving time-dependent DFT (TD-DFT) and time-dependent Hartree-Fock equations (TDHF) (the Casida equation) or the Bethe-Salpeter (BSE) equation. These approaches are used for obtaining the frequency-dependent dielectric function with the excitonic effects and can be based on the ground-state electronic structure in the DFT, hybrid-functional or GW approximations.

Solving Bethe-Salpeter and Casida equations

To take into account the excitonic effects or the electron-hole interaction, one has to use approximations beyond the independent-particle (IP) and the random-phase approximations (RPA). In VASP, it is done via the algorithm selected by ALGO=BSE or ALGO=TDHF, which essentially solves the same equations (Casida/Bethe-Salpeter) but mainly differ in the way the screening of the Coulomb potential is treated. The TDHF approach uses the exact-correlation kernel $f_{\rm {xc}}$, whereas BSE requires the $W(\omega \to 0)$ from a preceding GW calculation. Thus, in order to perform TDHF or BSE calculations, one has to provide the ground-state wavefunctions (WAVECAR) and the derivatives of the orbitals with respect to $k$ (WAVEDER). In addition, the BSE calculation requires files storing the screened Coulomb kernel produced in a GW calculation, i.e. WXXXX.tmp.

Both TDHF and BSE approaches require a preceding ground-state calculation, however, the TDHF does not need the preceding GW and can be performed with the DFT or hybrid-functional orbitals and energies.

Time-dependent Hartree-Fock

The TDHF calculations can be performed in two steps: the ground-state calculation and the optical absorption calculation. For example, optical absorption of bulk Si can be performed with a hybrid-functional electronic structure, where the number of bands is increased to include the relevant conduction bands:

System    = Si
ISMEAR    = 0 
SIGMA     = 0.05
NBANDS    = 16      ! or any larger desired value
ALGO      = D       ! Damped algorithm often required for HF type calculations, ALGO = Normal might work as well
LHFCALC   = .TRUE. 
AEXX      = 0.3 
HFSCREEN  = 0.2
LOPTICS   = .TRUE.  ! can also be done in an additional intermediate step

In the second step, the dielectric function is evaluated by solving the Casida equation

System    = Si
ISMEAR    = 0 
SIGMA     = 0.05
NBANDS    = 16     
ALGO      = TDHF
AEXX      = 0.3 
HFSCREEN  = 0.2

THDF/BSE calculations can be performed for non-spin-polarized, spin-polarized, as well as noncollinear (spin-orbit) cases. There is, however, one caveat. The local exchange-correlation kernel is approximated by the density-density part only. This makes predictions for spin-polarized systems less accurate than for non-spin-polarized systems.

Time-dependent DFT

Within the TD-DFT approximation, the Fock exchange is not included in the exchange-correlation kernel and the ladder diagrams are not taken into account. Hence, only the local contributions in $f_{\rm {xc}}$ are present

System    = Si
ISMEAR    = 0 
SIGMA     = 0.05
NBANDS    = 16     
ALGO      = TDHF
LFXC      = .TRUE.
AEXX      = 0.0 

Since the ladder diagrams are not included in the TD-DFT calculation, the resulting dielectric function lacks the excitonic effects.

Warning: Currently, in the local $f_{\rm {xc}}$, the correlation contribution is not taken into account and only the exchange terms are included in the LDA approximation.

Bethe-Salpeter equation

The BSE calculations require a preceding GW step to determine the screened Coulomb kernel $W_{GG'}(q,\omega \to 0)$. The details on GW calculations can be found in the practical guide to GW calculations. Here, we note that during the GW calculation, VASP writes this kernel into the following files

W0001.tmp, W0002.tmp, ..., W{NKPTS}.tmp

and

WFULL0001.tmp, WFULL0002.tmp, ..., WFULL{NKPTS}.tmp.

The files W000?.tmp store only the diagonal terms of the kernel and are fairly small, whereas the files WFULL000?.tmp store the full matrix. It is important to make sure in the GW step that the flag LWAVE = .TRUE. is set, so that the WAVECAR stores the one-electron GW energies and the one-electron orbitals, if the GW calculation is self-consistent.

For the self-consistent GW calculations the following flags should be added

LOPTICS   = .TRUE. 
LPEAD     = .TRUE.

in order to update the WAVEDER using finite differences (LPEAD = .TRUE.). The type of GW calculation is selected with the ALGO tag, which is discussed in great detail in the practical guide to GW calculations.

Once the GW step is completed, the BSE calculation can be performed using the following setup

System    = Si
NBANDS    = same as in GW calculation
ISMEAR    = 0
SIGMA     = 0.05
ALGO      = BSE
NBANDSO   = 4       ! determines how many occupied bands are used
NBANDSV   = 8       ! determines how many unoccupied (virtual) bands are used
OMEGAMAX  = maximum excitation energy 

Considering that quasiparticle energies in GW converge very slowly with the number of unoccupied bands and require large NBANDS, the number of bands included in the BSE calculation should be restricted explicitly by setting the occupied and unoccupied bands (NBANDSO and NBANDSV) included in the BSE Hamiltonian.

VASP tries to use sensible defaults, but it is highly recommended to check the OUTCAR file and make sure that the right bands are included. The tag OMEGAMAX specifies the maximum excitation energy of included electron-hole pairs and the pairs with the one-electron energy difference beyond this limit are not included in the BSE Hamiltonian. Hint: The convergence with respect to NBANDSV and OMEGAMAX should be thoroughly checked as the real part of the dielectric function, as well as the correlation energy, is usually very sensitive to these values, whereas the imaginary part of the dielectric function converges quickly.

At the beginning of the BSE calculation, VASP will try to read the WFULL000?.tmp files and if these files are not found, VASP will read the W000?.tmp files. For small isotropic bulk systems, the diagonal approximation of the dielectric screening may be sufficient and yields results very similar to the calculation with the full dielectric tensor WFULL000?.tmp. Nevertheless, for molecules and atoms as well as surfaces, the full screened Coulomb kernel is strictly required.

Both TDHF and BSE approaches write the calculated frequency-dependent dielectric function as well as the excitonic energies in the vasprun.xml file.

Calculations beyond Tamm-Dancoff approximation

The TDHF and BSE calculations beyond the Tamm-Dancoff approximation (TDA) can be performed by setting the ANTIRES = 2 in the INCAR file

System    = Si
NBANDS    = same as in GW calculation
ISMEAR    = 0
SIGMA     = 0.05
ALGO      = BSE  
ANTIRES   = 2   ! beyond Tamm-Dancoff 
NBANDSO   = 4 
NBANDSV   = 8

The flag LORBITALREAL = .TRUE. forces VASP to make the orbitals $\phi ({\bf {r}})$ real valued at the Gamma point as well as k-points at the edges of the Brillouin zone. This can improve the performance of BSE/TDHF calculations but it has to be used consistently with the ground-state calculation.

Calculations at finite wavevectors

VASP can also calculate the dielectric function at a ${\bf {q}}$-vector compatible with the k-point grid (finite-momentum excitons).

System       = Si
NBANDS       = same as in GW calculation
ISMEAR       = 0 
SIGMA        = 0.05
ALGO         = BSE  
ANTIRES      = 2 
KPOINT_BSE   = 3 -1 0 0  ! q-point index,  three integers
LORBITALREAL = .TRUE.
NBANDSO      = 4 
NBANDSV      = 8

The tag KPOINT_BSE sets the ${\bf {q}}$-point and the shift at which the dielectric function is calculated. The first integer specifies the index of the ${\bf {q}}$-point and the other three values shift the provided ${\bf {q}}$-point by an arbitrary reciprocal vector ${\bf {G}}$. The reciprocal lattice vector is supplied by three integer values $n_{i}$ with ${\bf {G}}=n_{1}{\bf {G}}_{1}+n_{2}{\bf {G}}_{2}+n_{3}{\bf {G}}_{3}$. This feature is only supported as of VASP.6 (in VASP.5 the feature can be enabled, but the results are erroneous).

Scaling of the BSE equation

The scaling of the BSE/Casida equation strongly limits its application to large systems. The main limiting factor is the diagonalization of the BSE/TDHF Hamiltonian. The rank of the Hamiltonian is

$N_{\rm {rank}}=N_{k}\times N_{c}\times N_{v}$,

where $N_{k}$ is the number of k-points in the Brillouin zone and $N_{c}$ and $N_{v}$ are the number of conduction and valence bands correspondingly. The diagonalization of the matrix scales cubically with the matrix rank, i.e. $N_{\rm {rank}}^{3}$. Despite the fact that this matrix diagonalization is usually the bottleneck for bigger systems, the construction of the BSE Hamiltonian also scales unfavorably and can play a dominant role in big systems, i.e.

$N_{k}\times N_{q}\times (N_{v}\times N_{v}\times N_{G}\times N_{c}\times N_{c})$,

where $N_{q}$ is the number of q-points and $N_{G}$ number of G-vectors.

Consistency tests

In order to verify the results obtained in the BSE calculation, one can perform a number of consistency tests.

First test: IP dielectric function

The BSE code can be used to reproduce the independent particle spectrum if the RPA and the ladder diagrams are switched off

LADDER   = .FALSE. 
LHARTREE = .FALSE.

This should yield exactly the same dielectric function as the preceding calculation with LOPTICS = .TRUE. We recommend to set the complex shift manually in the BSE as well as the preceding optics calculations, e.g. CSHIFT = 0.4. The dielectric functions produced in these calculations should be identical.

Second test: RPA dielectric function

The RPA/GW dielectric function can be used to verify the correctness of the RPA dielectric function calculated via the BSE algorithm. The RPA dielectric function in the BSE code can be calculated by switching off the ladder diagrams while keeping the RPA terms, i.e., the BSE calculation should be performed with the following tags

ANTIRES   = 2
LHARTREE  = .TRUE.
LADDER    = .FALSE.
CSHIFT    = 0.4

The same dielectric function should be obtained via the GW code by setting these flags

ALGO      = CHI 
NOMEGA    = 200
CSHIFT    = 0.4

Make sure that a large CSHIFT is selected as the GW code calculates the polarizability at very few frequency points. Note that the GW code does not use the TDA, so ANTIRES=2 is required for the TDHF/BSE calculation. In our experience, the agreement can be made practically perfect provided sufficient frequency points are used and all available occupied and virtual orbitals are included in the BSE step.

Third test: RPA correlation energy

The BSE code can be used to calculate the correlation energy via the plasmon equation. This correlation energy can be compared with the [[ACFDT/RPA calculations|RPA contributions]] to the correlation energies for each ${\bf {q}}$-point, which can be found in the OUTCAR file of the ACFDT/RPA calculation performed with ALGO = RPA:

q-point correlation energy      -0.232563      0.000000
q-point correlation energy      -0.571667      0.000000
q-point correlation energy      -0.176976      0.000000

For instance, if the BSE calculation is performed at the second ${\bf {q}}$-point

ANTIRES    = 2
LADDER     = .FALSE.
LHARTREE   = .TRUE.
KPOINT_BSE = 2 0 0 0

the same correlation energy should be found in the corresponding OUTCAR file:

plasmon correlation energy        -0.5716670828

For exact compatibility, ENCUT and ENCUTGW should be set to the same values in all calculations, while the head and wings of the dielectric matrix should not be included in the ACFDT/RPA calculations, i.e. remove the WAVEDER file prior to the ACFDT/RPA calculation. In the BSE/RPA calculation removing the WAVEDER file is not required. Furthermore, NBANDS in the ACFDT/RPA calculation must be identical to the number of included bands NBANDSO plus NBANDSV in the BSE/RPA, so that the same number of excitation pairs are included in both calculations. Also, the OMEGAMAX tag in the BSE calculation should not be set.

Common issues

If the dielectric matrix contains only zeros in the vasprun.xml file, the WAVEDER file was not read or is incompatible to the WAVEDER file. This requires a recalculation of the WAVEDER file. This can be achieved even after GW calculations using the following intermediate step:

ALGO      = Nothing
LOPTICS   = .TRUE.
LPEAD     = .TRUE.

The flag LPEAD = .TRUE. is strictly required and enforces a "numerical" differentiation of the orbitals with respect to $k$. Calculating the derivatives of the orbitals with respect to $k$ analytically is not possible at this point, since the Hamiltonian that was used to determine the orbitals is unknown to VASP.

Related tags and articles

ALGO, LOPTICS, LHFCALC, LRPA, LADDER, LHARTREE, NBANDSV, NBANDSO, OMEGAMAX, LFXC, ANTIRES

---