View source for Bethe-Salpeter-equations calculations

__TOC__

== Bethe Salpeter and Casida equations: a first look ==


VASP offers a powerful module for solving the Bethe-Salpeter<ref name="albrecht1998"/><ref name="rohlfing1998"/> 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 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 {{TAG|ALGO}} needs to be set to {{TAG|ALGO}}=''TDHF'' or {{TAG|ALGO}}=''BSE''. Internally {{TAG|ALGO}}=''TDHF'' and {{TAG|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 <math>f_{xc}</math>, screened exchange or by  <math>W(\omega\to 0)</math> calculated in a preceding GW calculations.

The first case, {{TAG|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:

 {{TAGBL|System}}  = Si
 {{TAGBL|NBANDS}} = 16 ! or any larger desired value
 {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
 {{TAGBL|ALGO}} = D   ! Damped algorithm often required for HF type calculations, ALGO = Normal might work as well
 {{TAGBL|LHFCALC}} = .TRUE. ; {{TAGBL|AEXX}} = 0.3 ; {{TAGBL|HFSCREEN}} = 0.2
 {{TAGBL|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 {{TAG|LHFCALC}} (in that case,  of course {{TAGBL|ALGO}} = NORMAL can be used).
In the second step, the dielectric function is evaluated by solving the Casida equation

 {{TAGBL|System}}  = Si
 {{TAGBL|NBANDS}} = 16
 {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
 {{TAGBL|ALGO}} = TDHF
 {{TAGBL|LHFCALC}} = .TRUE. ; {{TAGBL|AEXX}} = 0.3 ; {{TAGBL|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, {{TAGBL|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 {{TAG|LOPTICS}}=.TRUE., these calculations go beyond the independent particle 
approximation and include exchange correlation effects.


The second case, {{TAG|ALGO}}=''BSE'' is in general more accurate, but requires a preceding GW step (GW+BSE) to determine the screened Coulomb kernel <math>W_{GG'}(q,\omega \to 0 )</math>. 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 {{TAG|W000}}?.tmp store only the diagonal elements of the kernel and are fairly small, whereas the files {{TAG|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.

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

After the {{TAG|WAVECAR}} and {{TAG|WAVEDER}} files are written the actual GW calculation can be executed with 

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

Make certain that in the GW calculation, the flag  {{TAG|LWAVE}}=''.TRUE.'' is set, so that the {{TAG|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 

 {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE.

should be added in order to update the {{TAG|WAVEDER}} using finite differences ({{TAG|LPEAD}}=.TRUE.). A specific GW type is selected with the {{TAG|ALGO}} tag in the fifth line. For more information see [[GW calculations]].

The actual BSE calculation is initiated using:

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

It is strongly recommended to set the flags {{TAG|NBANDSO}} and {{TAG|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. {{TAG|NBANDSO}} controls how many occupied bands (below the Fermi-level) are included,
whereas {{TAG|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 {{TAG|OUTCAR}} file whether these agree with what the user desires. 
During the BSE calculations, VASP will first try to read the {{TAG|WFULL000}}?.tmp files and then, if these are missing, the {{TAG|W000}}?.tmp files. For small isotropic bulk systems, results with the more approximate files {{TAG|W000}}?.tmp are usually very similar to the results obtained using {{TAG|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, {{TAG|ALGO}}=''TDHF'' and {{TAG|ALGO}}=''BSE'', the dielectric function, as well as the calculated pair-excitation energies can be found in
the file {{TAG|vasprun.xml}}.

== Calculations beyond Tamm Dancoff approximation ==
The VASP BSE code can perform calculations beyond the Tamm Dancoff approximation (TDA) <ref name="sander2015"/>, by setting the 
{{TAG|ANTIRES}} =2 in the {{TAG|INCAR}} file e.g.

 {{TAGBL|System}}  = Si
 {{TAGBL|NBANDS}} = same as in GW calculation
 {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
 {{TAGBL|ALGO}} = BSE  
 {{TAGBL|ANTIRES}} = 2   ! go beyond TDA
 {{TAGBL|LORBITALREAL}} = .TRUE.
 {{TAGBL|NBANDSO}} = 4 ; {{TAGBL|NBANDSV}} = 8

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

== Calculations at finite wave vectors ==
For {{TAG|ANTIRES}} =2, VASP can also calculate the dielectric function at a q-vector compatible
with the k-point grid (finite momentum excitons). 

 {{TAGBL|System}}  = Si
 {{TAGBL|NBANDS}} = same as in GW calculation
 {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
 {{TAGBL|ALGO}} = BSE  
 {{TAGBL|ANTIRES}} = 2 
 {{TAGBL|KPOINT_BSE}} = 3 -1 0 0  # k-point index,  three integers
 {{TAGBL|LORBITALREAL}} = .TRUE.
 {{TAGBL|NBANDSO}} = 4 ; {{TAGBL|NBANDSV}} = 8

In the {{TAG|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 <math> \bf G</math>. The reciprocal lattice vector
is supplied by three integer values  <math> n_i</math> with <math> {\bf G}= n_1 {\bf G}_1+n_2 {\bf G}_2+n_3 {\bf G}_3</math>.
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

 {{TAGBL|LADDER}} = .FALSE. ;  {{TAGBL|LHARTREE}} = .FALSE.

This should yield exactly the same dielectric function as the preceding calculation using the {{TAG| LOPTICS}} = .TRUE.
We recommend to set the complex shift manually in the BSE and preceding optics calculations e.g. {{TAG | 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

 {{TAGBL|ANTIRES}} = 2
 {{TAGBL|LADDER}} = .FALSE. ;  {{TAGBL|LHARTREE}} = .TRUE.
 {{TAGBL|CSHIFT}} = 0.4

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

 {{TAGBL|ALGO}} = CHI ; {{TAGBL|NOMEGA}} = 200
 {{TAGBL|CSHIFT}} = 0.4

Just make sure that a large {{TAG|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 [[ACFDT/RPA calculations|RPA contributions]] to the correlation energies for each q-point (see {{TAG|OUTCAR}} from
{{TAG|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 

 {{TAGBL|ANTIRES}} = 2
 {{TAGBL|LADDER}} = .FALSE. ;  {{TAGBL|LHARTREE}} = .TRUE.
 {{TAGBL|KPOINT_BSE}} = 2 0 0 0

exactly the same correlation energy should be found in the {{TAG|OUTCAR}} file of the BSE calculation:

 plasmon correlation energy        -0.5716670828

For exact compatibility, {{TAG | ENCUT}} and {{TAG | 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 {{TAG|WAVEDER}} prior to the RPA calculation.

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

 {{TAGBL|ALGO}} = Nothing
 {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE.

The flag {{TAG|LPEAD}}=''.TRUE.'' is strictly required and enforces a "numerical" differentiation of the orbitals with respect to <math>k</math>. Calculating the derivatives of the orbitals with respect to <math>k</math> 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 ==
{{TAG|ALGO}},
{{TAG|LOPTICS}},
{{TAG|LHFCALC}},
{{TAG|LRPA}},
{{TAG|LADDER}},
{{TAG|LHARTREE}},
{{TAG|NBANDSV}},
{{TAG|NBANDSO}}


{{sc|BSE|Examples|Examples that use this tag}}

== References ==
<references>
<ref name="albrecht1998">[http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.80.4510 S. Albrecht, L. Reining, R. Del Sole, and G. Onida, Phys. Rev. Lett. 80, 4510 (1998).]</ref>
<ref name="rohlfing1998">[http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.81.2312 M. Rohlfing, S.G. Louie, Phys. Rev. Lett. 81, 2312 (1998).]</ref>

<ref name="sander2015">[https://doi.org/10.1103/PhysRevB.92.045209 S. Tobias, E. Maggio, and G. Kresse, Phys. Rev. B 92.4 (2015).]</ref>
</references>



----
[[The_VASP_Manual|Contents]]

[[Category:INCAR]][[Category:BSE]][[Category:Howto]]