Bethe-Salpeter-equations calculations: Difference between revisions

From VASP Wiki
VisualWikitext
No edit summary
 
(32 intermediate revisions by 4 users not shown)
Line 1: Line 1:
VASP offers a powerful module for solving the Bethe-Salpeter (BSE) equation{{cite|albrecht:prl:98}}{{cite|rohlfing:prl:98}}. The BSE can be 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. 
__TOC__
__TOC__


== Bethe Salpeter and Casida equations: a first look ==
== Solving Bethe-Salpeter equation ==
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/ACFDT: Correlation energy in the Random Phase Approximation|RPA]]). In VASP, it is done via the algorithm selected by {{TAG|ALGO}} = BSE. These essentially solves the same equations (Casida/Bethe-Salpeter) but differ in the way the screening of the Coulomb potential is treated. The TDHF approach uses the exact-correlation kernel <math>f_{\rm xc}</math>, whereas BSE requires the <math>W(\omega \to 0)</math> from a preceding ''GW '' calculation. Thus, in order to perform  TDHF or BSE calculations, one has to provide the ground-state orbitals ({{FILE|WAVECAR}}) and the derivatives of the orbitals with respect to <math>k</math> ({{FILE|WAVEDER}}). In addition, the BSE calculation requires files storing the screened Coulomb kernel produced in a ''GW'' calculation, i.e., {{FILE|Wxxxx.tmp}}.


In summary, 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.


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<ref name="albrecht1998"/><ref name="rohlfing1998"/> 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.
== Bethe-Salpeter equation calculation ==
We will first introduce a typical calculation steps for timedependent DFT
The BSE calculations require a preceding ''GW'' step to determine the screened Coulomb kernel <math>W_{GG'}(q,\omega \to 0 )</math>. The details on ''GW'' calculations can be found in the practical guide to [[Practical guide to GW calculations|''GW'' ]]  calculations. Here, we note that during the ''GW'' calculation, VASP writes this kernel into the following files
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 frequency dependent dielectric function (essentially
W0001.tmp, W0002.tmp, ..., W{NKPTS}.tmp
the Casida equations), however, the electron-hole ladder diagrams are either approximated by the the exchange correlation kernel <math>f_{xc}</math>, exact exchange, or by screened exchange <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:
and
 
{{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
WFULL0001.tmp, WFULL0002.tmp, ..., WFULL{NKPTS}.tmp.
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
The files {{FILE|Wxxxx.tmp}} store only the diagonal terms of the kernel and are fairly small, whereas the files {{TAG|WFULLxxxx.tmp}} store the full matrix. It is important to make sure in the ''GW'' step that the flag {{TAG|LWAVE}} = .TRUE. is set, so that the {{FILE|WAVECAR}} stores the one-electron ''GW'' energies and the one-electron orbitals, if the ''GW'' calculation is self-consistent.
{{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, 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.
For the self-consistent ''GW'' calculations the following flags should be added


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


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


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
Once the ''GW'' step is completed, the BSE calculation can be performed using the following setup


  W0001.tmp, W0002.tmp, W0003.tmp, ...
  {{TAG|SYSTEM}}    = Si
{{TAG|NBANDS}}    = same as in GW calculation
{{TAG|ISMEAR}}    = 0
{{TAG|SIGMA}}    = 0.05
{{TAG|ALGO}}      = BSE
{{TAG|NBANDSO}}  = 4      ! determines how many occupied bands are used
{{TAG|NBANDSV}}  = 8      ! determines how many unoccupied (virtual) bands are used
{{TAG|OMEGAMAX}}  = desired_maximum_excitation_energy


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


  WFULL0001.tmp, WFULL0002.tmp, WFULL0003.tmp, ...
VASP tries to use sensible defaults, but it is highly recommended to check the {{FILE|OUTCAR}} file and make sure that the right bands are included. The tag {{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 {{TAG|NBANDSV}} and {{TAG|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.


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.
At the beginning of the BSE calculation, VASP will try to read the {{FILE|WFULLxxxx.tmp}} files and if these files are not found, VASP will read the {{FILE|Wxxxx.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 {{TAG|WFULLxxxx.tmp}}. Nevertheless, for molecules and atoms as well as surfaces, the full-screened Coulomb kernel is strictly required.
The first one is a self-consistent DFT step similar to above but with significantly more bands, e.g.


{{TAGBL|System}}  = Si
Both TDHF and BSE approaches write the calculated frequency-dependent dielectric function as well as the excitonic energies in the {{TAG|vasprun.xml}} file.
{{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  
== Model BSE (mBSE) ==
BSE calculations can be performed using a model dielectric function{{cite|bokdam:scr:2016}}{{cite|tal:prr:2020}}. In this approach the calculation of the screened Coulomb potential is not required. Instead, the model dielectric function can be used to describe the screening of the Coulomb potential by setting the tag {{TAG|LMODELHF}} with parameters {{TAG|AEXX}} and {{TAG|HFSCREEN}}.


{{TAGBL|System}}  = Si
Model BSE calculation can be performed the following steps:
{{TAGBL|NBANDS}} = 64  ! large number of bands
# ground-state calculation
{{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
# GW calculation (optional in model BSE calculation)
{{TAGBL|LWAVE}} = .TRUE. ! update wavefunction
# optical absorption calculation via model BSE
{{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 example, an optical absorption calculation of bulk Si can be performed using a model dielectric function as described in Ref. {{cite|tal:prr:2020}}.
For scGW runs


  {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE.
  {{TAG|SYSTEM}}    = Si
{{TAG|ISMEAR}}    = 0
{{TAG|SIGMA}}    = 0.05
{{TAG|NBANDS}}    = 16      ! or any larger desired value
{{TAG|ALGO}}      = D      ! Damped algorithm often required for HF type calculations, {{TAG|ALGO}} = Normal might work as well
{{TAG|LHFCALC}}  = .TRUE.
{{TAG|LMODELHF}} = .TRUE.  
{{TAG|AEXX}}      = 0.083
{{TAG|HFSCREEN}}  = 1.22
{{TAG|LOPTICS}}   = .TRUE. ! can also be done in an additional intermediate step


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]].
In the second step, the dielectric function is evaluated by solving the Casida equation


The actual BSE calculation is initiated using:
{{TAG|SYSTEM}}    = Si
{{TAG|ISMEAR}}    = 0
{{TAG|SIGMA}}    = 0.05
{{TAG|NBANDS}}    = 16   
{{TAG|ALGO}}      = TDHF
{{TAG|IBSE}}      = 0
{{TAG|NBANDSO}}  = 4      ! number of occupied bands
{{TAG|NBANDSV}}  = 8      ! number of unoccupied bands
{{TAG|LHARTREE}}  = .TRUE.
{{TAG|LADDER}}    = .TRUE.
{{TAG|LFXC}}      = .FALSE. ! local xc kernel is disabled in mBSE
{{TAG|LMODELHF}}  = .TRUE.
{{TAG|AEXX}}      = 0.083
{{TAG|HFSCREEN}}  = 1.22


{{TAGBL|System}}  = Si
== Calculations beyond Tamm-Dancoff approximation ==
{{TAGBL|NBANDS}} = same as in GW calculation
The TDHF and BSE calculations beyond the Tamm-Dancoff approximation (TDA){{cite|sander:prb:15}} can be performed by setting the {{TAG|ANTIRES}} = 2 in the {{TAG|INCAR}} file
{{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
{{TAGBL|OMEGAMAX}} = desired_maximum_excitation_energy


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 advantageously to restrict the number of bands when solving the Cassida/BSE equations (the computational cost grows very steeply with {{TAG|NBANDSV}}). {{TAG|NBANDSO}} controls how many occupied bands (below the Fermi-level) are included,
{{TAG|SYSTEM}}       = Si
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. The tag {{TAG|OMEGAMAX}} can be used to reduce the number of excitation pairs further. VASP will include only those valence and conduction band pairs,
{{TAG|NBANDS}}       = same as in GW calculation
that have an one-electron energy difference (excitation energy) smaller than the supplied threshold. Make sure to perform tests with respect to {{TAG|NBANDSV}} and {{TAG|OMEGAMAX}}, the real part of
{{TAG|ISMEAR}}       = 0
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
{{TAG|SIGMA}}       = 0.05
to both parameters.
{{TAG|ALGO}}         = BSE 
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.
  {{TAG|ANTIRES}}     = 2      ! beyond Tamm-Dancoff
{{TAG|LORBITALREAL}} = .TRUE.  
{{TAG|NBANDSO}}     = 4
{{TAG|NBANDSV}}     = 8


In both cases, {{TAG|ALGO}}=''TDHF'' and {{TAG|ALGO}}=''BSE'', the frequency dependent dielectric function, as well as the calculated pair-excitation energies can be found in
The flag {{TAG|LORBITALREAL}} = .TRUE. forces VASP to make the orbitals <math> \phi({\bf r}) </math> 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 should be used consistently with the ground-state calculation.
the file {{TAG|vasprun.xml}}.


== Calculations beyond Tamm Dancoff approximation ==
== Calculations at finite wavevectors ==
The VASP BSE code can perform calculations beyond the Tamm Dancoff approximation (TDA) <ref name="sander2015"/>, by setting the
VASP can also calculate the dielectric function at a <math>{\bf q}</math>-vector compatible with the k-point grid (finite-momentum excitons).  
{{TAG|ANTIRES}} =2 in the {{TAG|INCAR}} file e.g.


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


This also works for TDHF calculations. 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.
The tag {{TAG|KPOINT_BSE}} sets the <math>{\bf q}</math>-point and the shift at which the dielectric function is calculated. The first integer specifies the index of the <math>{\bf q}</math>-point and the other three values shift the provided <math>{\bf q}</math>-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 as of VASP.6 (in VASP.5 the feature can be enabled, but the results are erroneous).


== Calculations at finite wave vectors ==
== Consistency tests ==
For {{TAG|ANTIRES}} =2, VASP can also calculate the dielectric function at a q-vector compatible
In order to verify the results obtained in the BSE calculation, one can perform a number of consistency tests.
with the k-point grid (finite momentum excitons).  


{{TAGBL|System}}  = Si
=== First test: IP dielectric function ===
{{TAGBL|NBANDS}} = same as in GW calculation
The BSE code can be used to reproduce the independent particle spectrum if the RPA and the ladder diagrams are switched off
{{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
{{TAG|LADDER}}   = .FALSE.  
at which the polarizability is supposed to be calculated. The other three
  {{TAG|LHARTREE}} = .FALSE.
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 ==
This should yield exactly the same dielectric function as the preceding calculation with {{TAG| LOPTICS}} = .TRUE.  We recommend to set the complex shift manually in the BSE as well as the preceding optics calculations, e.g.  {{TAG | CSHIFT}} = 0.4.  The dielectric functions produced in these calculations should be identical.
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
=== 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


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


This should yield exactly the same dielectric function as the preceding calculation using the {{TAG| LOPTICS}} = .TRUE.
The same dielectric function should be obtained via the ''GW'' code by setting these flags
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
{{TAG|ALGO}}      = CHI
BSE calculation simply add
{{TAG|NOMEGA}}    = 200
{{TAG|CSHIFT}}    = 0.4


{{TAGBL|ANTIRES}} = 2
Make sure that a large {{TAG|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 {{TAG|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.
{{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
=== 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 <math>{\bf q}</math>-point, which can be found in the {{FILE|OUTCAR}} file of the ACFDT/RPA calculation performed with {{TAG|ALGO}} = RPA:
{{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 (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 [[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.232563      0.000000
Line 156: Line 148:
  q-point correlation energy      -0.176976      0.000000
  q-point correlation energy      -0.176976      0.000000


For instance, if the second q-point is selected in the BSE calculations
For instance, if the BSE calculation is performed at the second <math>{\bf q}</math>-point  


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


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


  plasmon correlation energy        -0.5716670828
  plasmon correlation energy        -0.5716670828


For exact compatibility, {{TAG | ENCUT}} and {{TAG | ENCUTGW}} should be set to the same
For exact compatibility, {{TAG|ENCUT}} and {{TAG|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 {{FILE|WAVEDER}} file prior to the ACFDT/RPA calculation. In the BSE/RPA calculation removing the {{FILE|WAVEDER}} file is not required. Furthermore, {{TAG|NBANDS}} in the ACFDT/RPA calculation must be identical to the number of included bands {{TAG|NBANDSO}} plus {{TAG|NBANDSV}} in the BSE/RPA, so that the same number of excitation pairs are included in both calculations. Also, the {{TAG|OMEGAMAX}} tag in the BSE calculation should not be set.
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.
Furthermore, {{TAG|NBANDS}} must be identical to {{TAG|NBANDSO}} plus {{TAG|NBANDSV}}, so that the same number of excitation pairs are included
in the calculations, and avoid setting {{TAG | OMEGAMAX}} in the BSE calculation.


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


  {{TAGBL|ALGO}} = Nothing
  {{TAG|ALGO}}     = Nothing
  {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE.
  {{TAG|LOPTICS}}   = .TRUE.
{{TAG|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.
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 ==
== Related tags and articles ==
{{TAG|ALGO}},
{{TAG|ALGO}},
{{TAG|LOPTICS}},
{{TAG|LOPTICS}},
{{TAG|LPEAD}},
{{TAG|LHFCALC}},
{{TAG|LHFCALC}},
{{TAG|LRPA}},
{{TAG|LRPA}},
Line 188: Line 180:
{{TAG|NBANDSV}},
{{TAG|NBANDSV}},
{{TAG|NBANDSO}},
{{TAG|NBANDSO}},
{{TAG|OMEGAMAX}}
{{TAG|OMEGAMAX}},
 
{{TAG|LFXC}},
 
{{TAG|ANTIRES}},
See also: {{sc|BSE|Examples|Examples that use this tag}}, [[Time Evolution]]
{{TAG|NBSEEIG}},
{{FILE|BSEFATBAND}}


== References ==
== References ==
<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>
 
 


----
----
 
[[Category:Many-body perturbation theory]][[Category:Bethe-Salpeter equations]][[Category:Howto]]
[[Category:Many-Body Perturbation Theory]][[Category:BSE]][[Category:Howto]]

Latest revision as of 17:34, 7 February 2024

VASP offers a powerful module for solving the Bethe-Salpeter (BSE) equation[1][2]. The BSE can be 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 equation

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. These essentially solves the same equations (Casida/Bethe-Salpeter) but differ in the way the screening of the Coulomb potential is treated. The TDHF approach uses the exact-correlation kernel , whereas BSE requires the from a preceding GW calculation. Thus, in order to perform TDHF or BSE calculations, one has to provide the ground-state orbitals (WAVECAR) and the derivatives of the orbitals with respect to (WAVEDER). In addition, the BSE calculation requires files storing the screened Coulomb kernel produced in a GW calculation, i.e., Wxxxx.tmp.

In summary, 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.

Bethe-Salpeter equation calculation

The BSE calculations require a preceding GW step to determine the screened Coulomb kernel . 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 Wxxxx.tmp store only the diagonal terms of the kernel and are fairly small, whereas the files WFULLxxxx.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  = desired_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 WFULLxxxx.tmp files and if these files are not found, VASP will read the Wxxxx.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 WFULLxxxx.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.

Model BSE (mBSE)

BSE calculations can be performed using a model dielectric function[3][4]. In this approach the calculation of the screened Coulomb potential is not required. Instead, the model dielectric function can be used to describe the screening of the Coulomb potential by setting the tag LMODELHF with parameters AEXX and HFSCREEN.

Model BSE calculation can be performed the following steps:

  1. ground-state calculation
  2. GW calculation (optional in model BSE calculation)
  3. optical absorption calculation via model BSE

For example, an optical absorption calculation of bulk Si can be performed using a model dielectric function as described in Ref. [4]. 

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. 
LMODELHF  = .TRUE. 
AEXX      = 0.083
HFSCREEN  = 1.22
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
IBSE      = 0
NBANDSO   = 4       ! number of occupied bands
NBANDSV   = 8       ! number of unoccupied bands
LHARTREE  = .TRUE.
LADDER    = .TRUE.
LFXC      = .FALSE. ! local xc kernel is disabled in mBSE 
LMODELHF  = .TRUE. 
AEXX      = 0.083
HFSCREEN  = 1.22

Calculations beyond Tamm-Dancoff approximation

The TDHF and BSE calculations beyond the Tamm-Dancoff approximation (TDA)[5] 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
LORBITALREAL = .TRUE. 
NBANDSO      = 4 
NBANDSV      = 8

The flag LORBITALREAL = .TRUE. forces VASP to make the orbitals 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 should be used consistently with the ground-state calculation.

Calculations at finite wavevectors

VASP can also calculate the dielectric function at a -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 -point and the shift at which the dielectric function is calculated. The first integer specifies the index of the -point and the other three values shift the provided -point by an arbitrary reciprocal vector . The reciprocal lattice vector is supplied by three integer values with . This feature is only supported as of VASP.6 (in VASP.5 the feature can be enabled, but the results are erroneous).

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 RPA contributions to the correlation energies for each -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 -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 . 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 articles

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

References

  1. ↑ S. Albrecht, L. Reining, R. Del Sole, and G. Onida, Phys. Rev. Lett. 80, 4510-4513 (1998).
  2. ↑ M. Rohlfing and S. G. Louie, Phys. Rev. Lett. 81, 2312-2315 (1998).
  3. ↑ M. Bokdam, T. Sander, A. Stroppa, S. Picozzi, D. D. Sarma, C. Franchini, and G. Kresse, Sci. Rep. 6, 28618 (2016).
  4. ↑ a b A. Tal, P. Liu, G. Kresse, A. Pasquarello, Accurate optical spectra through time-dependent density functional theory based on screening-dependent hybrid functionals, Phys. Rev. Research 2, 032019 (2020)
  5. ↑ T. Sander, E. Maggio, and G. Kresse, Beyond the Tamm-Dancoff approximation for extended systems using exact diagonalization, Phys. Rev. B 92, 045209 (2015).
Retrieved from "https://www.vasp.at/wiki/index.php?title=Bethe-Salpeter-equations_calculations&oldid=23743"