Requests for technical support from the VASP group should be posted in the VASP-forum.
Bethe Salpeter and Casida equations: a first look
VASP offers a powerful module for solving time-dependent DFT (TD-DFT) and time-dependent Hartree-Fock equations (TDHF) (by solving the Casida equation) or the Bethe-Salpeter BSE equation. The routines determine the optical response function (i.e. frequency dependent dielectric function) including excitonic effects on top of standard DFT calculations, hybrid functionals or the GW. We will first introduce a typical calculation steps for timedependent DFT 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 frequency dependent dielectric function (essentially the Casida equations), however, the electron-hole ladder diagrams are either approximated by the the exchange correlation kernel , exact exchange, or by screened exchange 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, again omit that line). The present implementation can be used for non spin-polarized, 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 . To this end, VASP writes this kernel into the following files during the GW run
W0001.tmp, W0002.tmp, W0003.tmp, ...
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
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
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 OMEGAMAX = desired_maximum_excitation_energy
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 advantageously to restrict the number of bands when solving the Cassida/BSE equations (the computational cost grows very steeply with NBANDSV). 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. The tag OMEGAMAX can be used to reduce the number of excitation pairs further. VASP will include only those valence and conduction band pairs, that have an one-electron energy difference (excitation energy) smaller than the supplied threshold. Make sure to perform tests with respect to NBANDSV and OMEGAMAX, 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 with respect to both parameters. 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.
Calculations beyond Tamm Dancoff approximation
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. Furthermore, it is recommended to perform all preceding calculations using the flag LORBITALREAL =.TRUE.. This forces VASP to make the orbitals 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 . The reciprocal lattice vector is supplied by three integer values with . This feature is only supported in vasp.6 (in vasp.5 the feature can be selected, but the results are erroneous).
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
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
and run the BSE calculations. The exactly same dielectric function can be obtained from the GW code
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 (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, 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 (removal is not required in the BSE calculation). Furthermore, NBANDS must be identical to NBANDSO plus NBANDSV, so that the same number of excitation pairs are included in the calculations, and avoid setting OMEGAMAX in the BSE calculation.
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:
The flag LPEAD=.TRUE. is strictly required and enforces a "numerical" differentiation of the orbitals with respect to . Calculating the derivatives of the orbitals with respect to 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