Bethe-Salpeter-equations calculations

VASP offers a powerful and highly predictive module for Bethe-Salpeter^[1]^[2] BSE and time-dependent Hartree-Fock calculations (TDHF). It can be used to calculate the response function including excitonic effects on top of GW or hybrid functional calculations. We will first introduce a typical calculation and then report on the required flags in more detail.

To calculate spectra beyond the independent particle approximation and beyond the random phase approximation (RPA) the flag ALGO needs to be set to ALGO=TDHF or ALGO=BSE. Internally ALGO=TDHF and ALGO=BSE use identical routines to calculate the dielectric function, however, the electron-hole ladder diagrams are either approximated by the screened exchange or by $W(\omega \to 0)$ calculated in a preceding GW calculations.

The first case, ALGO=TDHF is simpler. The calculations need to be done in two steps. The first step is a standard DFT or hybrid functional calculation, where the number of bands is increased to include the relevant conduction bands in the calculation:

System  = Si
NBANDS = 16 ! or any larger desired value
ISMEAR = 0 ; SIGMA = 0.05
ALGO = D
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 Cassida equation

System  = Si
NBANDS = 16
ISMEAR = 0 ; SIGMA = 0.05
ALGO = TDHF
LHFCALC = .TRUE. ; AEXX = 0.3 ; HFSCREEN = 0.2

In this case the exchange kernel, as selected in the fifth line, should be identical to the previous groundstate calculations. The present implementation can be used for spin-polarized as well as non-collinear (spin-orbit) cases. There is, however, one caveat. The local exchange-correlation kernel is not exactly included and approximated by the density-density part only. This makes predictions for spin polarized systems less accurate then for non-spin polarized systems.

The second case, ALGO=BSE requires a preceding GW step (GW+BSE). Make certain that in the GW calculation, the flag LWAVE=.TRUE. is set, so that the WAVECAR file is updated to store the one-electron GW energies and possibly the one-electron orbitals as determined in scGW calculations. See GW Calculations for details about GW calculations.

The actual BSE calculation is initiated using:

System  = Si
NBANDS = same as in GW calculation
ISMEAR = 0 ; SIGMA = 0.05
ALGO = BSE
NBANDSO = 4 ; NBANDSV = 8 # determines how many occupied and virtual bands are used

It is strongly recommended to set the flags NBANDSO and NBANDSV: in GW calculations usually a very large number of bands is included. For the BSE, it is however desirable to restrict the number of bands when solving the Cassida/BSE equations. NBANDSO controls how many occupied bands (below the Fermi-level) are included, whereas NBANDSV determines the number of virtual (unoccupied) bands included above the Fermi-level. VASP tries to use sensible defaults, but it is always wise to check the OUTCAR file whether these agree with what the user desires.

For ALGO=BSE the particle-hole ladders are approximated by the $W(\omega \to 0)$ calculated in the preceding GW calculations. To this end, vasp writes the following files during the GW step

W0001.tmp  W0002.tmp  W0003.tmp

and

WFULL0001.tmp  WFULL0002.tmp  WFULL0003.tmp

The files W000?.tmp store only the diagonal elements of the screened exchange $W$ , and are therefore fairly small, whereas the files WFULL000?.tmp store the full matrix (the integer corresponds to the k-point index). During the BSE calculations, VASP will first try to read the WFULL000?.tmp files and then, if these are missing, the W000?.tmp files. For small isotropic bulk systems, results with the more approximate files W000?.tmp are usually very similar to the results obtained using WFULL000?.tmp, however, for molecules and atoms as well as surfaces it is strictly required to use the full screened Coulomb kernel $W$ .

In both cases, ALGO=TDHF and ALGO=BSE, the dielectric function, as well as the calculated pair-excitation energies can be found in the file vasprun.xml.

Calculations beyond Tamm Dancoff (TDA) approximation

The VASP BSE code can perform calculations beyond the TDA ^[3], by setting the ANTIRES =2 in the INCAR file e.g.

System  = Si
NBANDS = same as in GW calculation
ISMEAR = 0 ; SIGMA = 0.05
ALGO = BSE  
ANTIRES = 2 
LORBITALREAL = .TRUE.
NBANDSO = 4 ; NBANDSV = 8

Of course this also works for TDHF calculations. Furthermore, it is recommended to perform all preceding calculations using the flag LORBITALREAL =.TRUE.. This forces VASP to make the orbitals $\phi ({\bf {r}})$ real valued at the Gamma point as well as k-points at the edge of the Brillouin zone.

Calculations at finite wave vectors

For ANTIRES =2, VASP can also calculate the dielectric function at a 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  # k-point index,  three integers
LORBITALREAL = .TRUE.
NBANDSO = 4 ; NBANDSV = 8

In the KPOINT_BSE line, four integer values must be supplied. The first one specifies the index of the k-point at which the polarizability is supposed to be calculated. The other three values allow to shift the supplied k-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 in vasp.6 (in vasp.5 the feature can be selected, but the results are erroneous).

Always double check

There are many ways to double check the BSE code, and here we highlight some possible cross checks. First the BSE code can be executed to reproduce the independent particle spectrum by adding the lines

LADDER = .FALSE. ;  LHARTREE = .FALSE.

This should yield exactly the same dielectric function as the preceding calculation using the LOPTICS = .TRUE. We recommend to set the complex shift manually in the BSE and preceding optics calculations e.g. CSHIFT = 0.4, to obtain exactly the same results. This should yield absolutely identical results.

Second, one can use the RPA/GW routines to compare the dielectric function at the level of the RPA. In the BSE calculation simply add

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

and run the BSE calculations. The exactly same dielectric function can be obtained from the GW code

ALGO = CHI ; NOMEGA = 200
CSHIFT = 0.4

Just make sure that a large CSHIFT is selected (the GW code calculates the polarizability only at few frequency points). Note that the GW code does not use the Tamm-Dancoff approximation. In our experience the agreement can be made practically perfect (if sufficient frequency points are used, and if all occupied and virtual orbitals are included in the BSE calculation).

Furthermore, the TDHF/ BSE code calculates the correlation energy via the plasmon equation. This can be compared with the RPA contributions to the correlation energies for each q-point (see OUTCAR from 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 second k-point is selected in the BSE calculations

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

exactly the same correlation energy should be found in the OUTCAR file of the BSE calculation:

plasmon correlation energy        -0.5716670828

(for exact compatability ENCUT and ENCUTGW should be set to the same values in the calculations and the head and wings should not be included in the RPA calculations, e.g. rm WAVEDER prior to the RPA calculations).

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 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) at this point.

Related Tags and Sections

ALGO, LOPTICS, LHFCALC, LRPA, LADDER, LHARTREE, NBANDSV, NBANDSO

Examples that use this tag

References

Contents

[albrecht1998-1] S. Albrecht, L. Reining, R. Del Sole, and G. Onida, Phys. Rev. Lett. 80, 4510 (1998).

[rohlfing1998-2] M. Rohlfing, S.G. Louie, Phys. Rev. Lett. 81, 2312 (1998).

[sander2015-3] S. Tobias, E. Maggio, and G. Kresse, Phys. Rev. B 92.4 (2015).

[1]

[2]

[3]