Bethe-Salpeter-equations calculations

VASP offers a powerful and highly predictive module for solving the Bethe-Salpeter BSE equation, time-dependent DFT (TD-DFT) as well as time-dependent Hartree-Fock equations (TDHF) (the latter two by solving the Casida equation). It can be used to calculate the optical response function (i.e. frequency dependent dielectric function) including excitonic effects on top of GW, standard DFT calculations, or even hybrid functionals. We will first introduce a typical calculational steps 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 (essentially the Casida equations), however, the electron-hole ladder diagrams are either approximated by the the exchange correlation kernel $f_{xc}$, 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   ! 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

If course, it is also possible to start from a standard DFT calculations, by omitting the entire line starting with LHFCALC (in that case, of course ALGO = NORMAL can be used). In the second step, the dielectric function is evaluated by solving the Casida 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 ground state calculations (if the ground state calculation was not including exact exchange, simply again omit that line). 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.

To summarize, ALGO = TDHF allows to calculate the frequency dependent response function after a corresponding ground state calculation. This is done by solving the Casida equation. Compared to LOPTICS=.TRUE., these calculations go beyond the independent particle approximation and include exchange correlation effects.

The second case, ALGO=BSE is in general more accurate, but requires a preceding GW step (GW+BSE) to determine the screened Coulomb kernel $W_{GG'}(q,\omega \to 0 )$. To this end, VASP writes this kernel into the following files during the GW run

W0001.tmp, W0002.tmp, W0003.tmp, ...

and

WFULL0001.tmp, WFULL0002.tmp, WFULL0003.tmp, ...

The files W000?.tmp store only the diagonal elements of the kernel and are fairly small, whereas the files WFULL000?.tmp store the full matrix (the integer corresponds to the q-point index). A GW calculation in VASP consists of two steps. The first one is a self-consistent DFT step similar to above but with significantly more bands, e.g.

System  = Si
NBANDS = 64  ! large number of bands 
ISMEAR = 0 ; SIGMA = 0.05
ALGO = Exact ! exact diagonalization
LOPTICS = .TRUE. ! write derivative of wavefunction

After the WAVECAR and WAVEDER files are written the actual GW calculation can be executed with

System  = Si
NBANDS = 64  ! large number of bands 
ISMEAR = 0 ; SIGMA = 0.05
LWAVE = .TRUE. ! update wavefunction
ALGO = select GW flavor

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. For scGW runs

LOPTICS = .TRUE. ; LPEAD = .TRUE.

should be added in order to update the WAVEDER using finite differences (LPEAD=.TRUE.). A specific GW type is selected with the ALGO tag in the fifth line. For more information see 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. 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.

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 approximation

The VASP BSE code can perform calculations beyond the Tamm Dancoff approximation (TDA) , 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   ! go beyond TDA
LORBITALREAL = .TRUE.
NBANDSO = 4 ; NBANDSV = 8

This also works for TDHF calculations as well. 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 q-point at which the polarizability is supposed to be calculated. The other three values allow to shift the supplied 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 in vasp.6 (in vasp.5 the feature can be selected, but the results are erroneous).

Consistency checks

In the following we highlight three ways to check BSE results for consistency.

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. 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 implementation does not use the TDA. In our experience, however, the agreement can be made practically perfect provided sufficient frequency points are used and all occupied and virtual orbitals are included in the BSE step.

Third, 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 q-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 compatibility, ENCUT and ENCUTGW should be set to the same values in all calculations and the head and wings should not be included in the RPA calculations, e.g. rm WAVEDER prior to the RPA calculation.

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.

Related Tags and Sections

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

Examples that use this tag

References

---

Contents