Practical guide to GW calculations: Difference between revisions

From VASP Wiki
No edit summary
 
(175 intermediate revisions by 5 users not shown)
Line 1: Line 1:
__TOC__
The [[The GW approximation of Hedin's equations|GW approximation]] is an approximation to the self-energy. GW calculations are available as of VASP.5.X. For details on the implementation and use of the ''GW'' routines, we recommend the papers by Shishkin ''et al.'' {{cite|shishkin:prb:2006}}{{cite|shishkin:prb:2007}}{{cite|shishkin:prl:2007}} and Fuchs ''et al.''{{cite|fuchs:prb:2007}}
Available as of VASP.5.X. For details on the implementation and use of the ''GW'' routines we recommend the papers by Shishkin ''et al.''<ref name="shishkin-PRB74"/><ref name="shishkin-PRB75"/><ref name="shishkin:prl:07"/> and Fuchs ''et al.''<ref name="fuchs-PRB76"/>


= The GW approximation of Hedin's equations =
<span id="GW_in_one_go">
The GW method can be understood in terms of the following eigenvalue equation<ref name="HybertsenLouie"/>
= Single step procedure: GW in one go =


<span id="EVGW">
As of VASP.6.3 all GW approximations can be run in one single run by selecting the corresponding {{TAG|ALGO}} tag and omitting {{TAG|NBANDS}}), for instance like so
<math>
{{TAGBL|System}} = SiC
(T+V_{ext}+V_h)\phi_{n{\bf k}}({\bf r})+\int d{\bf r}\Sigma({\bf r},{\bf r}',\omega=E_{n{\bf k}})\phi_{n{\bf k}}({\bf r}') = E_{n{\bf k}}\phi_{n{\bf k}}({\bf r})
{{TAGBL|ALGO}} = EVGW0, QPGW0, EVGW, QPGW, GW0R or GWR # use an algorithgm described below
</math>
{{TAGBL|NELMGW}} = 1,2,.. # number of self-consistency cycles
</span>
{{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05  ! small sigma is required to avoid partial occupancies
 
{{TAGBL|LOPTICS}} = .TRUE.  # for insulators, omit for metals
Here <math>T</math> is the kinetic energy, <math>V_{ext}</math> the external potential of the nuclei, <math>V_h</math> the Hartree potential and <math>E_{n{\bf k}}</math> the quasiparticle energies with orbitals <math>\phi_{n{\bf k}}</math>. In contrast to DFT, the exchange-correlation potential is replaced by the many-body self-energy <math>\Sigma</math> and should be obtained together with the Green's function <math>G</math>, the irreducible polarizability <math>\chi</math>, the screened Coulomb interaction <math>W</math> and the irreducible vertex function <math>\Gamma</math> in a self-consistent procedure. For completeness, these equations are<ref name="Hedin"/>
Note, {{TAG|NBANDS}} must not be present in the INCAR to select this procedure.
 
{{NB|important|In older versions a two step procedure is required, where the first step is always a DFT calculation and the second step the actual GW calculation.}}
<span id="G">
The two-step procedure is described below.
<math>
G(1,2)=G_0(1,2)+\int d(3,4) G_0(1,3)\Sigma(3,4)G(4,2)
</math>
 
</span>
<span id="chi">
<math>
\chi(1,2)=\int d(3,4) G(1,3)G(4,1)\Gamma(3,4;2)
</math>
 
</span>
<span id="W">
<math>
W(1,2)=V(1,2)+\int d(3,4) V(1,3)\chi(3,4)W(4,2)
</math>
</span>
 
<span>
<span id="Sigma">
<math>
\Sigma(1,2)=\int d(3,4) G(1,3)\Gamma(3,2;4)W(4,1)
</math>
</span>
 
<span>
<span id="I">
<math>
\Gamma(1,2;3)=\delta(1,2)\delta(1,3)+\int d(4,5,6,7)\frac{\delta\Sigma(1,2)}{\delta G(4,5)}G(4,6)G(7,5)\Gamma(6,7;3)
</math>
</span>
 
Here the common notation <math>1=({\bf r}_1,t_1)</math> was adopted and <math>V</math> denotes the bare Coulomb interaction. Note, that these equations are exact and provide an alternative to the Schrödinger equation for the many-body problem. Nevertheless, approximations are necessary for realistic systems. The most popular one is the GW approximation and is obtained by neglecting the equation for the vertex function and using the bare vertex <math>\Gamma(1,2;3)=\delta(1,2)\delta(1,3)</math> instead.  
 
This means that the equations for the polarizability and self-energy reduce to
 
</span>
<span id="chi_GW">
<math>
\chi(1,2)=G(1,2)G(2,1)
</math>
 
<span>
<span id="Sigma_GW">
<math>
\Sigma(1,2)=G(1,2)W(2,1)
</math>
</span>
 
while the equations for the Green's function and the screened potential remain the same. However, in practice, these equations are usually solved in reciprocal space in the frequency domain
 
</span>
<span id="W_reci">
<math>
W_{{\bf G}{\bf G}'}({\bf q},\omega)=\left[\delta_{{\bf G}{\bf G}'}-\chi_{{\bf G}{\bf G}'}({\bf q},\omega)V_{{\bf G}{\bf G}'}({\bf q})\right]^{-1}V_{{\bf G}{\bf G}'}({\bf q})
</math>
 
<span>
<span id="G_GW">
<math>
G_{{\bf G}{\bf G}'}({\bf q},\omega)=\left[\delta_{{\bf G}{\bf G}'}-\Sigma_{{\bf G}{\bf G}'}({\bf q},\omega)G^{0}_{{\bf G}{\bf G}'}({\bf q})\right]^{-1}G^{0}_{{\bf G}{\bf G}'}({\bf q})
</math>
</span>


Here, in the first iteration <math>G^{0}</math> is the non-interacting Green's function
== Caveats ==
The single-step GW procedure performs a DFT step internally with an exact diagonalization of the Kohn-Sham Hamiltonian using the maximum available {{TAGBL|NBANDS}} supported for the chosen {{TAGBL|ENCUT}} value.
Consequently, a large number of unoccupied bands is initialized with random plane-wave coefficients. In rare cases, this yields two linearly dependent column vectors in the Hamiltonian and results in LAPACK errors like "ZPOTRF fails". These errors can be prevented using the two-step GW procedure as described below. Furthermore, one can "ramp up" NBANDS to the maximum value by repeatedly restarting the DFT calculation from a pre-converged {{TAGBL|WAVECAR}} with fewer bands.


<span id="SpectralG">
<span id="DFT_Step_for_GW">
<math>
G^{(0)}({\bf r},{\bf r}',\omega)=\sum_{n{\bf k}}\frac{\phi_{n{\bf k}}^{*(0)}
({\bf r})\phi^{(0)}_{n{\bf k}}
({\bf r}')}{\omega-\epsilon_{n{\bf k}}}</math>
</span>


with <math>\phi_{n{\bf k}}</math> being a set of one-electron orbitals and <math>E_{n{\bf k}}^{(0)}</math> the corresponding energies.That means that GW calculations require a first guess for the one-electron eigensystem, which is usually taken from a preceding DFT step.
= First step: DFT calculation =
 
''GW'' calculations always require a one-electron basis set. Usually this set is obtained from a standard DFT calculation and written into the {{FILE|WAVECAR}} file and can be calculated for instance the following INCAR file:
In the following we discuss those supported by VASP and how the calculations are performed in practice. 
  {{TAGBL|System}} = SiC
 
  {{TAGBL|NBANDS}} = 512
= General outline of a GW calculation =
 
As of VASP.6 all GW approximations can be selected directly via the {{TAG|ALGO}} tag (omitting {{TAG|NBANDS}}) without a preceding DFT calculation. However, for older versions a two step procedure is required, where the first step is always a DFT calculation.
The actual GW calculation is performed in the second step.
 
== First step: DFT calculation ==
''GW'' calculations always require the calculation of a standard DFT {{FILE|WAVECAR}} file in an initial step, using for instance the following INCAR file:
 
  {{TAGBL|System}} = Si
  {{TAGBL|NBANDS}} = 96
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05  ! small sigma is required to avoid partial occupancies
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05  ! small sigma is required to avoid partial occupancies
  {{TAGBL|LOPTICS}} = .TRUE.
  {{TAGBL|LOPTICS}} = .TRUE.
 
Note, that a significant number of empty bands is required for ''GW'' calculations, so that it might be better to perform the calculations in two steps: first a standard ground-state calculation with few unoccupied orbitals only,
Note, that the a significant number of empty bands is required for ''GW'' calculations, so that it might be better to perform the calculations in two steps: first a standard grounstate calculation with few unoccupied orbitals only,
  {{TAGBL|System}} = SiC ground-state occupied orbitals
 
  {{TAGBL|System}} = Si groundstate occupied orbitals
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05  ! small sigma is required to avoid partial occupancies
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05  ! small sigma is required to avoid partial occupancies
  {{TAGBL|EDIFF}} = 1E-8              ! required tight tolerance for groundstate orbitals
  {{TAGBL|EDIFF}} = 1E-8              ! required tight tolerance for ground-state orbitals
 
and, second, a calculation of a large number of unoccupied orbitals
and, second, a calculation of a large number of unoccupied orbitals
 
  {{TAGBL|System}}  = SiC unoccupied orbitals
  {{TAGBL|System}}  = Si unoccupied orbitals
  {{TAGBL|ALGO}} = Exact              ! use exact diagonalization of the Hamiltonian
  {{TAGBL|ALGO}} = Exact              ! use exact diagonalization of the Hamiltonian
  {{TAGBL|NELM}} = 1                  ! since we are already converged stop after one step
  {{TAGBL|NELM}} = 1                  ! since we are already converged stop after one step
  {{TAGBL|NBANDS}} = 96                
  {{TAGBL|NBANDS}} = 512              ! maybe even larger                
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05  ! small sigma is required to avoid partial occupancies
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05  ! small sigma is required to avoid partial occupancies
  {{TAGBL|LOPTICS}} = .TRUE.
  {{TAGBL|LOPTICS}} = .TRUE.
Furthermore, note that the flag {{TAG|LOPTICS}}=.TRUE. is required to write the file {{FILE|WAVEDER}}, which contains the derivative of the orbitals with respect to '''k'''.
This derivative is used to construct the head and wings of the dielectric matrix employing '''k'''&middot;'''p''' perturbation theory and is important to accelerate k-point convergence for insulators and semiconductors. {{NB|warning|For metals, in general, we recommend omitting the {{TAG|LOPTICS}} tag and removing the {{FILE|WAVEDER}} file from the directory.}}


Furthermore note that the flag {{TAG|LOPTICS}}=.TRUE. is required in order to write the file {{FILE|WAVEDER}}, which contains the derivative of the orbitals with respect to '''k'''.
== Optional: Use Hybrid functionals ==
This derivative is used to construct the head and wings of the dielectric matrix employing '''k'''&middot;'''p''' perturbation theory and is important to accelerate k-point convergence for insulators and semiconductors. For metals, in general, we recommend to omit the {{TAG|LOPTICS}} tag and remove the {{FILE|WAVEDER}} file from the directory.
Optionally, one can start a ''GW'' calculation from a [[:Category:Hybrid_functionals|hybrid functional]], such as HSE.  
 
For hybrid functionals, the two step procedure will accordingly involve the following {{FILE|INCAR}} files. In the first step, converged [[List_of_hybrid_functionals#HSE|HSE03]] orbitals are determined (see [[List_of_hybrid_functionals|here]] for a selection of available hybrid functionals):
=== Optional: Use Hybrid functionals ===
  {{TAGBL|System}}  = SiC ground-state occupied orbitals
Optionally, one can start a ''GW'' calculation from a hybrid functional, such as HSE.  
For hybrid functionals, the two step procedure will accordingly involve the following {{FILE|INCAR}} files. In the first step, converged [[Hartree-Fock_and_HF/DFT_hybrid_functionals#HSE|HSE03]] orbitals are determined (usually hybrid functional calculations should be preceeded by standard DFT calculations, we have not documented this step here, see the [[Howto on Hybrid functional calculations]]):
 
  {{TAGBL|System}}  = Si groundstate occupied orbitals
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ALGO}} = Damped ; {{TAGBL|TIME}} = 0.5  ! or {{TAGBL|ALGO}} = Conjugate
  {{TAGBL|ALGO}} = Damped ; {{TAGBL|TIME}} = 0.5  ! or {{TAGBL|ALGO}} = Conjugate
  {{TAGBL|LHFCALC}} = .TRUE. ; {{TAGBL|AEXX}} = 0.25 ; {{TAGBL|HFSCREEN}} = 0.3  
  {{TAGBL|LHFCALC}} = .TRUE. ; {{TAGBL|AEXX}} = 0.25 ; {{TAGBL|HFSCREEN}} = 0.3  
  {{TAGBL|EDIFF}} = 1E-6      ! required tight tolerance for groundstate orbitals
  {{TAGBL|EDIFF}} = 1E-6      ! required tight tolerance for ground-state orbitals
 
Secondly, determine the HSE03 orbitals for unoccupied states:
Secondly, determine the HSE03 orbitals for unoccupied states:
 
  {{TAGBL|System}}  = SiC unoccupied orbitals
  {{TAGBL|System}}  = Si unoccupied orbitals
  {{TAGBL|NBANDS}} = 512      ! maybe even larger
  {{TAGBL|NBANDS}} = 96
  {{TAGBL|ALGO}} = Exact
  {{TAGBL|ALGO}} = Exact
  {{TAGBL|NELM}} = 1          ! since we are already converged stop after one step
  {{TAGBL|NELM}} = 1          ! since we are already converged stop after one step
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|LHFCALC}} = .TRUE. ; {{TAGBL|AEXX}} = 0.25 ; {{TAGBL|HFSCREEN}} = 0.3  
  {{TAGBL|LHFCALC}} = .TRUE. ; {{TAGBL|AEXX}} = 0.25 ; {{TAGBL|HFSCREEN}} = 0.3  
  {{TAGBL|LOPTICS}} = .TRUE.
  {{TAGBL|LOPTICS}} = .TRUE. # for insulators
</span>


== Step 2: quasi-particle GW ==
= Second step: GW calculation=
The actual ''GW'' calculation is done in a second step. Here different ''GW'' flavors are possible and are selected with the {{TAG|ALGO}} tag.  
The actual ''GW'' calculation is done in a second step. Here different ''GW'' flavors are possible and are selected with the {{TAG|ALGO}} tag.  
Note that as of VASP.6 the GW {{TAG|ALGO}} tags have been renamed, see [[ALGO#GW tags in VASP.6|here]] for VASP.5.X tags.


<span id="G0W0">
<span id="G0W0">
=== Single shot quasi-particle energies ===  
== Single shot quasiparticle energies: G<sub>0</sub>W<sub>0</sub> ==
This is the simplest ''GW'' calculation and computationally the most efficient one.  
This is the simplest ''GW'' calculation and computationally the most efficient one.  
A single-shot calculation is often referred to as G<sub>0</sub>W<sub>0</sub> and calculates the quasi-particle energies from a single ''GW'' iteration by neglecting all off-diagonal matrix elements of the self-energy and employing a Taylor expansion of the self-energy around the DFT energies <math>E_{n{\bf k}}^{(0)}</math>. The corresponding equation becomes
A single-shot calculation is often referred to as G<sub>0</sub>W<sub>0</sub> and calculates the quasiparticle energies from a single ''GW'' iteration by neglecting all off-diagonal matrix elements of the self-energy and employing a Taylor expansion of the self-energy around the DFT energies <math>E_{n{\bf q}}^{(0)}</math>. The corresponding equation becomes


<span id="G0W0_eigenvalues">
<span id="G0W0_eigenvalues">
Line 164: Line 83:


In VASP, G<sub>0</sub>W<sub>0</sub> calculations are selected using an {{TAG|INCAR}} file such as
In VASP, G<sub>0</sub>W<sub>0</sub> calculations are selected using an {{TAG|INCAR}} file such as
 
  {{TAGBL|System}}  = SiC
  {{TAGBL|System}}  = Si
  {{TAGBL|NBANDS}} = 512
  {{TAGBL|NBANDS}} = 64
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|NELM}} = 1 ; {{TAGBL|ALGO}} = EVGW0 ! use "GW0" for VASP.5.X
  {{TAGBL|NELMGW}} = 1 ! use NELM for VASP.6.2 and older
{{TAGBL|ALGO}} = EVGW0 ! use "GW0" for VASP.5.X
  {{TAGBL|NOMEGA}} = 50
  {{TAGBL|NOMEGA}} = 50
{{NB|mind|Convergence with respect to the number of empty bands {{TAG|NBANDS}} and with respect to the number of frequencies {{TAG|NOMEGA}} must be checked carefully.}}


Convergence with respect to the number of empty bands {{TAG|NBANDS}} and with respect to the number of frequencies {{TAG|NOMEGA}} must be checked carefully.
To avoid complicated inter-nested tests, we recommend calculating all orbitals that the plane-wave basis set allows to calculate (except for simple tests). For further reading, please consult the section on {{TAG|ENCUTGW}}.


To avoid complicated inter-nested tests, we recommend to calculate all orbitals that the plane wave basis set allows to calculate (except for simple tests). For further reading please consult the section on {{TAG|ENCUTGW}}.
After a successful G<sub>0</sub>W<sub>0</sub> run, VASP will write the quasiparticle energies into the {{TAG|OUTCAR}} file for a set of {{TAG|NBANDSGW}} bands for every k-point in the Brillouin zone. Look for lines similar to  
 
QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
After a successful G<sub>0</sub>W<sub>0</sub> run, VASP will write the quasi-particle energies into the {{TAG|OUTCAR}} file for a set of {{TAG|NBANDSGW}} bands for every k-point in the Brillouin zone. Look for lines similar to  
for sc-GW calculations column KS-energies equals QP-energies in previous step
 
and V_xc(KS)=  KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1  + <V_x>^1)
  k-point  1 :      0.0000    0.0000    0.0000
  k-point  1 :      0.0000    0.0000    0.0000
   band No.  KS-energies  QP-energies  sigma(KS)  V_xc(KS)    V^pw_x(r,r')  Z            occupation Imag(sigma)
   band No.  KS-energies  QP-energies  sigma(KS)  V_xc(KS)    V^pw_x(r,r')  Z            occupation Imag(sigma)
   
   
       1      -6.4074     -6.9316     -11.1153     -10.4469     -17.4566       0.7843       2.0000      1.3376
       1      -7.1627     -8.3040     -14.5626     -12.7276     -21.6682       0.6219       2.0000      1.2037
       2       5.5804      5.1918     -11.9505     -11.4385     -13.4545       0.7589       2.0000      0.0124
       2     -2.0901      -3.4347     -15.7660     -14.2799     -21.7439       0.9048       2.0000      0.6914
       3       5.5804      5.1918     -11.9505     -11.4385     -13.4545       0.7589       2.0000      0.0124
       3     -2.0901      -3.4347     -15.7660     -14.2799     -21.7439       0.9048       2.0000      0.6914
       4       5.5804      5.1918     -11.9505     -11.4385     -13.4545       0.7589       2.0000      0.0124
       4     -2.0901      -3.4347     -15.7660     -14.2799     -21.7439       0.9048       2.0000      0.6914
       5      8.1035      8.3477      -9.7546     -10.0823      -5.3852       0.7449       0.0000     -0.0459
       5      0.4603      -0.4663    -13.7603     -12.5200    -18.1532       0.7471       2.0000       0.2167
       6       8.1035      8.3477      -9.7546    -10.0823      -5.3852       0.7449      0.0000     -0.0459
       6      0.4603     -0.4663    -13.7603     -12.5200    -18.1532       0.7471       2.0000       0.2167
      7      8.1035      8.3477      -9.7546     -10.0823      -5.3852       0.7449       0.0000     -0.0459
The first column is the band index and the third column denotes the quasiparticle energies <math>E_{n{\bf q}}</math>. Column two, four, five and seven refer to the DFT energies <math>E_{n{\bf q}}^{(0)}</math>, diagonal matrix elements of the self-energy <math>\langle \phi^{(0)}_{n{\bf q}}|\Sigma(\omega=E_{n{\bf q}}^{(0)}) |\phi^{(0)}_{n{\bf q}}\rangle</math>, the exchange-correlation potential <math>\langle \phi^{(0)}_{n{\bf q}}|V_{xc} |\phi^{(0)}_{n{\bf q}}\rangle</math> and the renormalization factor <math>Z_{n{\bf q}}</math> defined above, respectively.
      8      8.9494      9.3129    -10.5554    -11.0374      -5.5531      0.7543      0.0000      -0.0829


The first column is the band index and the thrid column denotes the quasi-particle energies <math>E_{n{\bf q}}</math>. Column two, four, five and seven refer to the DFT energies <math>E_{n{\bf q}}^{(0)}</math>, diagonal matrix elements of the self-energy <math>\langle \phi^{(0)}_{n{\bf q}}|\Sigma(\omega=E_{n{\bf q}}^{(0)}) |\phi^{(0)}_{n{\bf q}}\rangle</math>, the exchange-correlation potential <math>\langle \phi^{(0)}_{n{\bf q}}|V_{xc} |\phi^{(0)}_{n{\bf q}}\rangle</math> and the renormalization factor <math>Z_{n{\bf q}}</math> defined above, respectively.
<span id="gw0">
</span>


<span id="gw0">
== Partially self-consistent calculations: EVGW<sub>0</sub> ==
=== Partially self-consistent EVGW<sub>0</sub> calculations ===
In most cases, the ''best'' results (''i.e.'', closest to experiment) are obtained by iterating only <math>G</math> via the spectral representation
In most cases, the ''best'' results (''i.e.'', closest to experiment) are obtained by iterating only <math>G</math> via the spectral representation


Line 203: Line 121:
</span>
</span>


but keeping <math>W</math> fixed to the initial DFT level <math>W_0</math>. This method goes back to Hybertsen and Louie<ref name="HybertsenLouie"/> and can be achieved in two ways.  
but keeping <math>W</math> and the orbitals <math>\phi^{(0)}_{n{\bf q}}</math> fixed to the initial DFT level. This method goes back to Hybertsen and Louie {{cite|hybertsen:prb:1986}} and can be achieved in two ways.  


If the spectral method is not selected ({{TAG|LSPECTRAL}}=.FALSE., requiring much more compute time), the quasi-particle (QP) shifts are iterated automatically four times, and one finds four sets of QP shifts in the {{FILE|OUTCAR}} file. The first one corresponds to the G<sub>0</sub>W<sub>0</sub> case. The {{FILE|INCAR}} file is simply:
If the spectral method is not selected ({{TAG|LSPECTRAL}}=.FALSE., requiring much more compute time), the quasiparticle (QP) shifts are iterated automatically four times, and one finds four sets of QP shifts in the {{FILE|OUTCAR}} file. The first one corresponds to the G<sub>0</sub>W<sub>0</sub> case. The {{FILE|INCAR}} file is simply:
 
  {{TAGBL|System}} = SiC
  {{TAGBL|System}} = Si
  {{TAGBL|NBANDS}} = 512
  {{TAGBL|NBANDS}} = 64
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ALGO}} = EVGW0 ! use "GW0" in VASP.5.X
  {{TAGBL|ALGO}} = EVGW0 ! use "GW0" in VASP.5.X
  {{TAGBL|LSPECTRAL}} =.FALSE.
  {{TAGBL|LSPECTRAL}} =.FALSE.
{{NB|tip|In self-consistent GW calculations, convergence with the number of updated bands {{TAG|NBANDSGW}} must be checked carefully.}}


For technical reasons, it is not possible to iterate <math>G</math> in this manner if {{TAG|LSPECTRAL}}=.TRUE. is set in the {{FILE|INCAR}} file (this is the default). In this case, an iteration number must be supplied in the {{FILE|INCAR}} file using the {{TAG|NELM}}-tag. Usually three to four iterations are sufficient to obtain accurate QP shifts.
For technical reasons, it is not possible to iterate <math>G</math> in this manner if {{TAG|LSPECTRAL}}=.TRUE. is set in the {{FILE|INCAR}} file (this is the default). In this case, an iteration number must be supplied in the {{FILE|INCAR}} file using the {{TAG|NELMGW}} tag. Usually, three to four iterations are sufficient to obtain accurate QP shifts.
 
  {{TAGBL|System}} = SiC
  {{TAGBL|System}} = Si
  {{TAGBL|NBANDS}} = 512
  {{TAGBL|NBANDS}} = 64
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ALGO}} = EVGW0 ! use "GW0" in VASP.5.X
  {{TAGBL|ALGO}} = EVGW0 ! use "GW0" in VASP.5.X
  {{TAGBL|NELM}} = 4
  {{TAGBL|NELMGW}} = 4 ! use NELM in VASP.6.2 and older
</span>
</span>


<span id="scgw0">
The results are found again in the {{TAGBL|OUTCAR}} file
=== Self-consistent QPGW<sub>0</sub> calculations ===  
QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 4
k-point  1 :      0.0000    0.0000    0.0000
  band No. old QP-enery  QP-energies  sigma(KS)  T+V_ion+V_H  V^pw_x(r,r')  Z            occupation Imag(sigma)
      1      -8.6924      -8.7107    -14.2871      5.5647    -21.6681      0.6076      2.0000      1.1648
      2      -3.4692      -3.4806    -15.6742      12.1894    -21.7437      0.7304      2.0000      0.6351
      3      -3.4692      -3.4806    -15.6742      12.1894    -21.7437      0.7304      2.0000      0.6351
      4      -3.4692      -3.4806    -15.6742      12.1894    -21.7437      0.7304      2.0000      0.6351
      5      -0.6957      -0.7006    -13.6827      12.9802    -18.1531      0.7264      2.0000      0.2769
      6      -0.6957      -0.7006    -13.6827      12.9802    -18.1531      0.7264      2.0000      0.2769
      7      -0.6957      -0.7006    -13.6827      12.9802    -18.1531      0.7264      2.0000      0.2769
In contrast to single shot GW calculations, the second column represent now the QP-energies from the previous iteration. 
<span id="qpgw0">
 
== Partially self-consistent quasiparticle calculations: QPGW<sub>0</sub> ==  
If non diagonal components of the self-energy (in the orbital basis) should be included use {{TAG|ALGO}}=QPGW0. The following setting can be used:
If non diagonal components of the self-energy (in the orbital basis) should be included use {{TAG|ALGO}}=QPGW0. The following setting can be used:
 
  {{TAGBL|System}} = SiC
  {{TAGBL|System}} = Si
  {{TAGBL|NBANDS}} = 512
  {{TAGBL|NBANDS}} = 64
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ALGO}} = QPGW0 ! or "scGW0" for VASP.5.2.11 and older  
  {{TAGBL|ALGO}} = QPGW0 ! or "scGW0" for VASP.5.2.11 and older  
  {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE. ! ommit this lines for metals
  {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE. ! ommit this lines for metals
  {{TAGBL|NELM}} = 4
  {{TAGBL|NELMGW}} = 4 ! use NELM for VASP.6.2 and older
 
In this case, the orbitals are updated as well by constructing a hermitian (energy independent) approximation to the self-energy {{cite|shishkin:prl:2007}}. The "static" COHSEX approximation can be selected by setting {{TAG|NOMEGA}} = 1 {{cite|bruneval:prb:06}}. To improve convergence to the ground-state, the charge density (and the charge density only) is mixed using a Kerker type mixing in VASP.5.3.2 and more recent versions (see {{TAG|IMIX}}). The mixing parameters {{TAG|AMIX}}, {{TAG|BMIX}}, {{TAG|AMIX_MAG}}, {{TAG|BMIX_MAG}}, {{TAG|AMIN}} can be adjusted, if convergence problems are encountered.
In this case, the orbitals are updated as well by constructing a hermitian (energy independent) approximation to the self-energy.<ref name="shishkin:prl:07"/> The "static" COHSEX approximation can be selected by setting {{TAG|NOMEGA}} = 1.<ref name="bruneval:prb:06"/> To improve convergence to the groundstate, the charge density (and the charge density only) is mixed using a Kerker type mixing in VASP.5.3.2 and more recent versions (see {{TAG|IMIX}}). The mixing parameters {{TAG|AMIX}}, {{TAG|BMIX}}, {{TAG|AMIX_MAG}}, {{TAG|BMIX_MAG}}, {{TAG|AMIN}} can be adjusted, if convergence problems are encountered.
We strongly urge the user to monitor convergence by inspecting the lines <pre>charge density residual</pre> in the {{FILE|OUTCAR}} files.
We strongly urge the user to monitor convergence by inspecting the lines <pre>charge density residual</pre> in the {{FILE|OUTCAR}} files.


Alternatively the mixing may be switched off by setting {{TAG|IMIX}}=0 and controlling the step width for the orbitals using the parameter {{TAG|TIME}} (which defaults to 0.4). This selects a fairly sophisticated damped MD algorithm, that is also used for DFT methods when {{TAG|ALGO}}=Damped. In general, this method is more reliable for metals and materials with strong charge sloshing.
Alternatively, the mixing may be switched off by setting {{TAG|IMIX}}=0 and controlling the step width for the orbitals using the parameter {{TAG|TIME}} (which defaults to 0.4). This selects a fairly sophisticated damped MD algorithm that is also used for DFT methods when {{TAG|ALGO}}=Damped. This method is generally more reliable for metals and materials with strong charge sloshing.


After every iteration VASP writes following lines into the {{TAG|OUTCAR}} file
After every iteration, VASP writes the following lines into the {{TAG|OUTCAR}} file
QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
    GWSYM:  cpu time  15.8978: real time  15.9528
   
   
  k-point  1 :      0.0000    0.0000    0.0000
  k-point  1 :      0.0000    0.0000    0.0000
   band No. DFT-energies  QP-energies  QP-e(diag)  sigma(DFT)    Z            occupation
   band No. DFT-energies  QP-energies  QP-e(diag)  sigma(DFT)    Z            occupation
   
   
       1      -6.4074     -6.9604     -6.9316     -7.0758       0.7843       2.0000
       1      -7.1626     -8.4217     -8.3038     -8.9978       0.6219       2.0000
       2       5.5804      5.1885      5.1918      5.0683       0.7589       2.0000
       2     -2.0899      -3.4394      -3.4347      -3.5765       0.9047       2.0000
       3       5.5804      5.1885      5.1918      5.0683       0.7589       2.0000
       3     -2.0899      -3.4394      -3.4347      -3.5765       0.9047       2.0000
       4       5.5804      5.1885      5.1918      5.0683       0.7589       2.0000
       4     -2.0899      -3.4394      -3.4347      -3.5765       0.9047       2.0000
       5      8.1035      8.3022      8.3477      8.4313       0.7449       0.0000
       5      0.4604      -0.4787      -0.4663      -0.7800       0.7471       2.0000
       6      8.1035      8.3022      8.3477      8.4313       0.7449       0.0000
       6      0.4604      -0.4787      -0.4663      -0.7800       0.7471       2.0000
       7      8.1035      8.3022      8.3477      8.4313       0.7449       0.0000
       7      0.4604      -0.4787      -0.4663      -0.7800       0.7471       2.0000
       8      8.9494       9.3052       9.3129       9.4313       0.7543       0.0000
       8      5.1013       4.1883       4.2149       3.9518       0.7711       2.0000
For the first iteration, here, the fourth column should be identical to the third column of the G<sub>0</sub>W<sub>0</sub> results discussed above. The third column reports the quasiparticle energies obtained from including the off-diagonal matrix elements in the eigenvalue equation.  


For the first iteration, here the forth column should be identical to the third column of the EVGW0 results discussed above, while the third column reports the quasi-particle energies obtained from including the off-diagonal matrix elements in the eigenvalue equation.  
=== Caveats ===
The ''QPGW0'' (or ''scGW0'' in VASP.5.2.11 and older) must be used with great caution, particularly in combination with symmetry.
Symmetry is handled in a rather sophisticated manner. Specifically, only the minimal number of required combinations of ''q'' and ''k'' points is considered. In this case, symmetry must be applied to restore the full star of ''q''. This is done by determining degenerate eigenvalue/eigenvector pairs and restoring their symmetry according to their irreducible representation. Although the procedure is generally relatively reliable, it fails to work properly if the degenerate states do not possess eigenvalues that are sufficiently close due to insufficient convergence in the preceding DFT calculations. That is because states are treated as degenerate if, and only if, their eigenenergies are within 0.01 eV.


==== Caveats ====
For large supercells with low symmetry, we strongly recommend switching off symmetry.
The ''QPGW0'' (or ''scGW0'' in VASP.5.2.11 and older) must be used with great caution, in particular, in combination with symmetry.
Symmetry is handled in a rather sophisticated manner, specifically, only the minimal number of required combination of ''q'' and ''k'' points is considered. In this case, symmetry must be applied to restore the full star of ''q''. This is done by determining degenerate eigenvalue/eigenvector pairs and restoring their symmetry according to their irreducible representation. Although the procedure is generally rather reliable, it fails to work properly if the degenerate states do not posses eigenvalues that are sufficiently close, due to insufficient convergence in the preceding DFT calculations. States are treated as degenerate if, and only if, their eigenenergies are within 0.01 eV.
 
For large supercells with low symmetry, we strongly recommend to switch off symmetry.
</span>
</span>


<span id="scgw">
<span id="qpgw">


=== Self-consistent EVGW and QPGW calculations ===
== Self-consistent EVGW and QPGW calculations ==
Self-consistent ''QPGW'' calculations are only supported in a QP picture. As for ''QPGW''<sub>0</sub>, it is possible to update the eigenvalues only ({{TAG|ALGO}}=''EVGW'' or ''GW'' for VASP.5.X), or the eigenvalues and one-electron orbitals ({{TAG|ALGO}}=''QPGW'' or ''scGW'' in VASP.5.2.11 and older). In all cases, a QP picture is maintained, ''i.e.'', satellite peaks (shake ups and shake downs) can not be accounted for in the self-consistency cycle. Self-consistent ''QPGW'' calculations can be performed by simply repeatedly calling VASP using:
Self-consistent ''QPGW'' calculations are only supported in a QP picture. As for ''QPGW''<sub>0</sub>, it is possible to update the eigenvalues only ({{TAG|ALGO}}=''EVGW'' or ''GW'' for VASP.5.X), or the eigenvalues and one-electron orbitals ({{TAG|ALGO}}=''QPGW'' or ''scGW'' in VASP.5.2.11 and older). In all cases, a QP picture is maintained, ''i.e.'', satellite peaks (shake ups and shake downs) can not be accounted for in the self-consistency cycle. Self-consistent ''QPGW'' calculations can be performed by simply repeatedly calling VASP using:
 
  {{TAGBL|System}} = SiC
  {{TAGBL|System}} = Si
  {{TAGBL|NBANDS}} = 512
  {{TAGBL|NBANDS}} = 64
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ALGO}} = EVGW    ! "GW" in VASP.5.X, eigenvalues only  or alternatively
  {{TAGBL|ALGO}} = EVGW    ! "GW" in VASP.5.X, eigenvalues only  or alternatively
  {{TAGBL|ALGO}} = QPGW    ! "scGW" in VASP.5.2.11 and older, eigenvalues and one electron orbitals
  {{TAGBL|ALGO}} = QPGW    ! "scGW" in VASP.5.2.11 and older, eigenvalues and one electron orbitals
For QPGW0 or QPGW, nondiagonal terms in the Hamiltonian are accounted for, ''e.g.'' the linearized QP equation is diagonalized, and the one-electron orbitals are updated {{cite|shishkin:prl:2007}}.
Analogously to the other ''GW'' algorithms, the low scaling GW implementation can be selected by using ''EVGWR'' or ''QPGWR'' respectively.
Alternatively (and preferably), the user can specify an electronic iteration counter using {{TAG|NELMGW}} ({{TAG|NELM}} in VASP.6.2 and older):
 
  {{TAGBL|System}} = SiC
For QPGW0 or QPGW non diagonal terms in the Hamiltonian are accounted for, ''e.g.'' the linearized QP equation is diagonalized, and the one electron orbitals are updated.<ref name="shishkin:prl:07"/>
  {{TAGBL|NBANDS}} = 512
Alternatively (and preferably), the user can specify an electronic iteration counter using {{TAG|NELM}}:
 
  {{TAGBL|System}} = Si
  {{TAGBL|NBANDS}} = 64
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|NELM}} = 3
  {{TAGBL|NELMGW}} = 3   ! use NELM in VASP.6.2 and older
  {{TAGBL|ALGO}} = EVGW  ! "GW" in VASP.5.X  
  {{TAGBL|ALGO}} = EVGW  ! "GW" in VASP.5.X  
  # or   
  # or   
  {{TAGBL|ALGO}} = QPGW  ! "scGW" in VASP.5.2.11 and older
  {{TAGBL|ALGO}} = QPGW  ! "scGW" in VASP.5.2.11 and older
In this case, the one-electron energies (=QP energies) are updated 3 times (starting from the DFT eigenvalues) in both G and W. For {{TAG|ALGO}}=''QPGW'' (or {{TAG|ALGO}}=''scGW'' in VASP.5.2.11 and older), the one electron energies and one electron orbitals are updated 3 times {{cite|shishkin:prl:2007}}. As for {{TAG|ALGO}} = ''QPGW0'' (or ''scGW0'' in vasp.5.2.11 and older), the "static" COHSEX approximation can be selected by setting {{TAG|NOMEGA}}=1 {{cite|bruneval:prb:06}}.


In this case, the one-electron energies (=QP energies) are updated 3 times (starting from the DFT eigenvalues) in both G and W. For {{TAG|ALGO}}=''QPGW'' (or {{TAG|ALGO}}=''scGW'' in VASP.5.2.11 and older), the one electron energies and one electron orbitals are updated 3 times.<ref name="shishkin:prl:07"/> As for {{TAG|ALGO}} = ''QPGW0'' (or ''scGW0'' in vasp.5.2.11 and older), the "static" COHSEX approximation can be selected by setting {{TAG|NOMEGA}}=1.<ref name="bruneval:prb:06"/>
To improve convergence to the ground-state, the charge density is mixed using a Kerker type mixing starting with VASP.5.3.2 (see {{TAG|IMIX}}). The mixing parameters {{TAG|AMIX}}, {{TAG|BMIX}}, {{TAG|AMIX_MAG}}, {{TAG|BMIX_MAG}}, {{TAG|AMIN}} can be adjusted, if convergence problems are encountered.
 
Alternatively, the mixing may be switched off by setting {{TAG|IMIX}}=0 and controlling the step width for the orbitals using the parameter {{TAG|TIME}} (which defaults to 0.4). This selects a fairly sophisticated damped MD algorithm that is also used for DFT methods when {{TAG|ALGO}}=Damped. This method is generally more reliable for metals and materials with strong charge sloshing.
To improve convergence to the groundstate, the charge density is mixed using a Kerker type mixing starting with VASP.5.3.2 (see {{TAG|IMIX}}). The mixing parameters {{TAG|AMIX}}, {{TAG|BMIX}}, {{TAG|AMIX_MAG}}, {{TAG|BMIX_MAG}}, {{TAG|AMIN}} can be adjusted, if convergence problems are encountered.
Alternatively the mixing may be switched off by setting {{TAG|IMIX}}=0 and controlling the step width for the orbitals using the parameter {{TAG|TIME}} (which defaults to 0.4). This selects a fairly sophisticated damped MD algorithm, that is also used for DFT methods when {{TAG|ALGO}}=Damped. In general, this method is more reliable for metals and materials with strong charge sloshing.


==== Caveats ====
Additional information about this method is found [[The GW approximation of Hedin's equations#GWLimitations|here]].
Fully self-consistent QPGW calculations with an update of the orbitals in <math>G</math> and <math>W</math><ref name="shishkin:prl:07"/> require significant care and are prone to diverge (QPGW0 calculations are usually less critical). As discussed, above, one can select this mode using:


  {{TAGBL|System}} = Si
=== Caveats ===
  {{TAGBL|NBANDS}} = 64
Fully self-consistent QPGW calculations with an update of the orbitals in <math>G</math> and <math>W</math>{{cite|shishkin:prl:2007}} require significant care and are prone to diverge (QPGW0 calculations are usually less critical). As discussed, above, one can select this mode using:
  {{TAGBL|System}} = SiC
  {{TAGBL|NBANDS}} = 512
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ALGO}} = QPGW  ! or "scGW" in VASP.5.2.11 and older, eigenvalues and one-electron orbitals
  {{TAGBL|ALGO}} = QPGW  ! or "scGW" in VASP.5.2.11 and older, eigenvalues and one-electron orbitals
  {{TAGBL|NELM}} = number of steps
  {{TAGBL|NELMGW}} = number of steps ! use NELM for VASP.6.2 and older
 
However, one ''caveat'' applies to this case: when the orbitals are updated, the derivatives of the orbitals with respect to <math>k</math> (stored in the {{FILE|WAVEDER}} file) will become incompatible with the orbitals.
However, one ''caveat'' applies to this case: when the orbitals are updated, the derivatives of the orbitals with respect to <math>k</math> (stored in the {{FILE|WAVEDER}} file) will become incompatible with the orbitals.
This can cause severe problems and convergence to the incorrect solution. For metals, we recommend to avoid using the {{FILE|WAVEDER}} file alltogether ({{TAG|LOPTICS}}=.TRUE. should not be used in the preparatory DFT runs).
This can cause severe problems and convergence to the incorrect solution.  
{{NB|warning|For metals, in general, we recommend omitting the {{TAG|LOPTICS}} tag and removing the {{FILE|WAVEDER}} file from the directory.}}


For insulators, VASP (version 5.3.2 or higher) can update the {{FILE|WAVEDER}} file in each electronic iteration if the finite difference method is used to calculate the first derivative of the orbitals with respect to <math> k </math>:
For insulators, VASP (version 5.3.2 or higher) can update the {{FILE|WAVEDER}} file in each electronic iteration if the finite difference method is used to calculate the first derivative of the orbitals with respect to <math> k </math>:
 
  {{TAGBL|System}} = SiC
  {{TAGBL|System}} = Si
  {{TAGBL|NBANDS}} = 512
  {{TAGBL|NBANDS}} = 64
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ALGO}} = QPGW ! "scGW" in VASP.5.2.11 and older, eigenvalues and one-electron orbitals
  {{TAGBL|ALGO}} = QPGW ! "scGW" in VASP.5.2.11 and older, eigenvalues and one-electron orbitals
  {{TAGBL|NELM}}  = 10
  {{TAGBL|NELMGW}}  = 10 ! use NELM in VASP.6.2 and older
  {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE.
  {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE.
The combination <tt>{{TAG|LOPTICS}}=.TRUE.; {{TAG|LPEAD}}=.TRUE.</tt> is required since <math>\frac{\delta H-\epsilon_{n{\bf k}}S}{\delta {k}_i}</math> is not available for ''GW'' like methods.
The combination <tt>{{TAG|LOPTICS}}=.TRUE.; {{TAG|LPEAD}}=.TRUE.</tt> is required since <math>\frac{\delta H-\epsilon_{n{\bf k}}S}{\delta {k}_i}</math> is not available for ''GW'' like methods.
{{TAG|LPEAD}}=.TRUE. circumvents this problems by calculating the derivatives of the orbitals using numerical differentiation on the finite k-point grid (this option is presently limited to insulators).
{{TAG|LPEAD}}=.TRUE. circumvents this problem by calculating the derivatives of the orbitals using numerical differentiation on the finite k-point grid (this option is presently limited to insulators).


Vertex corrections are presently not documented. This is a feature still under construction, and we recommend to collaborate with the Vienna group if you are desperately in need of that feature.
Vertex corrections are presently not documented. This is a feature still under construction, and we recommend collaborating with the Vienna group if you desperately need that feature.
 
Note that as of VASP.6 the GW {{TAG|ALGO}} tags have been renamed, see [[ALGO#GW tags in VASP.6|here]] for VASP.5.X tags.
</span>
</span>


<span id="LowGW">
<span id="LowGW">
= Low scaling GW algorithms =
= Low scaling GW algorithms =
</span>
The GW implementations in VASP described in the papers of Shishkin ''et al.''{{cite|shishkin:prb:2006}} {{cite|shishkin:prb:2007}} avoids storage of the Green's function <math>G</math> as well as Fourier transformations between time and frequency domain entirely. That is, all calculations are performed solely on the real frequency axis using Kramers-Kronig transformations for convolutions in the equation of <math>\chi</math> and <math>\Sigma</math> in reciprocal space.


The GW implementations in VASP described in the papers of Shishkin ''et al.''<ref name="shishkin-PRB74"/><ref name="shishkin-PRB75"/> avoids storage of the Green's function <math>G</math> as well as Fourier transformations between time and frequency domain entirely. That is, all calculations are performed solely on the real frequency axis using Kramers-Kronig transformations for convolutions in the equation of <math>\chi</math> and <math>\Sigma</math> in reciprocal space.  
As of VASP.6 a new cubic scaling GW algorithm {{cite|liu:prb:2016}} (called space-time implementation in the following) can be selected. This approach follows the idea of Rojas ''et al.'' {{cite|rojas:prl:1995}} and performs the GW self-consistency cycle on imaginary time <math>t\to i\tau</math> and imaginary frequency axes <math>\omega\to i\omega</math>.
{{NB|tip|Using the low-scaling GW algorithm also calculates the total energy in the Random Phase approximation (RPA), which is described in a [[ACFDT/RPA calculations#ACFDTR/RPAR|separate article]].}}


As of VASP.6 a new cubic scaling GW algorithm<ref name="liu"/> (called space-time implementation in the following) can be selected. This approach follows the idea of Rojas ''et al.''<ref name="rojas"/> and performs the GW self-consistency cycle on imaginary time <math>t\to i\tau</math> and imaginary frequency axes <math>\omega\to i\omega</math>.  
<span id="G0W0R">
== Low scaling, single shot GW calculations: G<sub>0</sub>W<sub>0</sub>R==
The low-scaling analogue of G<sub>0</sub>W<sub>0</sub> is selected with {{TAGBL|ALGO}}=G0W0R.
In contrast to the [[GW calculations#G0W0|single-shot GW calculations on the real-axes]],
here the self-energy <math>\Sigma = G_0 W_0 </math> is determined on the imaginary frequency axis.
To this end, the overall scaling is reduced by one order of magnitude and is cubic with respect to the system size,
because a small value for {{TAG|NOMEGA}} can be used (usually <20).  


Following the [[ACFDT/RPA calculations#ACFDTR/RPAR| low scaling ACFDT/RPA algorithms]] the space-time implementation determines first, the non-interacting Green's function on the imaginary time axis in real space
This algorithm evaluates:
* Single-shot GW quasiparticle energies (from an analytical continuation of the self-energy to the real axis){{cite|liu:prb:2016}}


<math>G({\bf r},{\bf r}',i\tau)=-\sum_{n{\bf k}}\phi_{n{\bf k}}^{*(0)}({\bf r}) \phi_{n{\bf k}}^{(0)}({\bf r}') e^{-(\epsilon_{n{\bf k}}-\mu)\tau}\left[\Theta(\tau)(1-n_{n{\bf k}})-\Theta(-\tau)n_{n{\bf k}}\right]</math>
* Natural orbitals from the first order change of the density matrix (i.e. <math>G_0 \Sigma G_0</math>), see the {{TAG|NATURALO}} tag for more information {{cite|ramberger:jcp:2019}}.
 
{{NB|mind|This selection ignores {{TAG|NELMGW}}.}}
Here <math>\Theta</math> is the step function and <math>n_{n{\bf k}}</math> the occupation number of the state <math>\phi_{n{\bf k}}^{(0)}</math>. Because the Green's function is non-oscillatory on the imaginary time axis it can be represented on a coarse grid <math>\tau_{m}</math>, where the number of time points can be selected in VASP via the {{TAG|NOMEGA}} tag. Usually 16 points are sufficient for insulators and small band gap systems, while metals require finer girds with roughly 24 points.  
Following {{TAGBL|INCAR}} file selects the low-scaling GW algorithm:
 
{{TAGBL|System}} = SiC
Subsequently, the irreducible polarizability is calculated from a contraction of two imaginary time Green's functions
{{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
{{TAGBL|LOPTICS}} = .TRUE. 
{{TAGBL|ALGO}} = G0W0R
{{TAGBL|NOMEGA}} = 12 ! small number of frequencies necessary
Search the {{TAGBL|OUTCAR}} file for the following lines
  QP shifts evaluated in KS or natural orbital/ Bruckner basis
  k-point  1 :      0.0000    0.0000    0.0000
  band No.  KS-energies  sigma(KS)    QP-e(linear)    Z        QP-e(zeros)    Z        occupation    Imag(E_QP)    QP_DIFF TAG
 
        1      -7.1627      -8.6732      -8.2451      0.7166      -8.2346      0.7026      2.0000      -1.3101      0.0000  2
        2      -2.0901      -3.4155      -3.0350      0.7129      -3.0272      0.7011      2.0000      -0.5582      -0.0000  2
        3      -2.0901      -3.4155      -3.0350      0.7129      -3.0272      0.7011      2.0000      -0.5582      0.0000  2
        4      -2.0901      -3.4155      -3.0350      0.7129      -3.0272      0.7011      2.0000      -0.5582      -0.0000  2
        5      0.4603      -0.8219      -0.4904      0.7414      -0.4814      0.7273      2.0000      -0.1902      0.0000  2
        6      0.4603      -0.8219      -0.4904      0.7414      -0.4814      0.7273      2.0000      -0.1902      -0.0000  2
Here column four is obtained by a linearization of the self-energy around the Kohn-Sham energies (second column) and can be compared to the third column of [[GW calculations#G0W0|single-shot GW calculations on the real axis]].
Column six represents another set of QP-energies that is obtained from the roots of the following equation


<math>
<math>
\chi({\bf r},{\bf r}',i\tau_m) = -G({\bf r},{\bf r}',i\tau_m)G({\bf r}',{\bf r},-i\tau_m)
\langle \phi^{(0)}_{n{\bf q}}| T + V_{ext}+V_h+ \Sigma(\omega) | \phi^{(0)}_{n{\bf q}}\rangle -\omega =0
</math>
</math>


Afterwards, the same compressed Fourier transformation as for the [[ACFDT/RPA calculations#ACFDTR/RPAR| low scaling ACFDT/RPA algorithm]] is employed to obtain the irreducible polarizability in reciprocal space on the imaginary frequency axis <math>\chi({\bf r},{\bf r}',i\tau_m) \to \chi_{{\bf G}{\bf G}'}({\bf q},i \omega_n) </math>.
These roots represent the poles of the Green's function in the spectral representation.  
The next step is the computation of the screened potential


<span id=OutputGW>
=== Output description ===
The meaning of each column is explained briefly in the following.
* <code>band No.</code> the band index of KS orbital at given k-point
* <code>KS-energies</code> eigenenergies corresponding to band index
* <code>sigma(KS)</code> diagonal matrix elements of self-energy evaluated at KS energies
* <code>QP-e(linear)</code> quasiparticle energies obtained from linearizing frequency dependence of diagonal self-energy around KS energies
* <code>Z</code> renormalization factor obtained from five-point stencil for derivative of self-energy w.r.t. frequency
* <code>QP-e(zeros)</code> quasiparticle energies obtained from full frequency dependence of self-energy, i.e. real part of complex pole <math>\omega</math> of Green's function
* <code>Z</code> renormalization factor obtained from central difference for derivative of self-energy w.r.t. frequency
* <code>occupation</code> occupation number for band at given k-point
* <code>Imag(E_QP)</code> imaginary part of complex pole <math>\omega</math>, i.e. measure for inverse lifetime of quasi-particle
* <code>QP_DIFF</code> difference of QP energies (of linearized self-energy) obtained from Eq. 77 of Liu et. al.{{cite|liu:prb:2016}} and M. Grumets thesis{{cite|grumet:thesis:2017}}.
</span>
</span>
<span id="W2_reci">
<math>
W_{{\bf G}{\bf G}'}({\bf q},\omega_m)=\left[\delta_{{\bf G}{\bf G}'}-\chi_{{\bf G}{\bf G}'}({\bf q},\omega_m)V_{{\bf G}{\bf G}'}({\bf q})\right]^{-1}V_{{\bf G}{\bf G}'}({\bf q})
</math>
</span>
</span>


followed by the inverse Fourier transform <math>W_{{\bf G}{\bf G}'}({\bf q},i \omega_n) \to \chi({\bf r},{\bf r}',i\tau_m)</math> and the calculation of the self-energy
<span id=RPAFORCES>


<math>
== Optional: RPA Forces ==
\Sigma({\bf r},{\bf r}',i\tau_m) = -G({\bf r},{\bf r}',i\tau_m)W({\bf r}',{\bf r},-i\tau_m)
Optionally, RPA forces can be calculated by adding following line to the {{TAGBL|INCAR}}:
</math>
  {{TAGBL|LRPAFORCE}} = .TRUE.
After the QP-energies, VASP performs a linear-response calculation that is required for the RPA forces.{{cite|ramberger:prl:118}}
Following data block in the {{TAGBL|OUTCAR}} file can be found after a successful run:
POSITION                                      TOTAL RPA FORCE (eV/Angst)
-----------------------------------------------------------------------------------
      0.17542    -0.22348      0.17542        -0.292069      7.581315    -0.292069
      1.12850      1.31044      1.12850        0.304683    -7.605527      0.304683
-----------------------------------------------------------------------------------
    total drift:                                0.012614    -0.024212      0.012614
SUGGESTED UPDATED POSCAR (direct coordinates)  step
-----------------------------------------------------------------------------------
  -0.00958461  -0.00958461  0.13485779        0.04179056  0.04179056  0.00283088
  0.25787833  0.25787833  0.22191754        -0.04337198  -0.04337198  0.00431513
{{NB|warning|Currently RPA forces for metallic systems are not supported.}}
</span>
<span id="EVGW0R">
== Low scaling, partially self-consistent GW calculations: EVGW<sub>0</sub>R==
{{NB|mind|available as of vasp 6.4.0}}
The low-scaling analogue of EVGW<sub>0</sub> is selected with {{TAGBL|ALGO}}=EVGW0R.
Following {{TAGBL|INCAR}} file selects this algorithm:
{{TAGBL|System}} = SiC
{{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
{{TAGBL|LOPTICS}} = .TRUE. 
{{TAGBL|ALGO}} = EVGW0R
{{TAGBL|NELMGW}} = 4  ! number of iterations in G
{{TAGBL|NOMEGA}} = 12 ! small number of frequencies necessary
After each iteration, a similar block of data as for {{TAG|ALGO}}=G0W0R calculations is written to {{FILE|OUTCAR}} showing the {{TAG|NBANDSGW}} updated quasi-particle energies (poles) of the Green's function. 
</span>
</span>
<span id="scGW0R">
== Partially self-consistent GW calculations: GW<sub>0</sub>R ==
The space-time implementation allows for true self-consistent GW calculations. That is, the solution of the Dyson equation for the Green's function can be obtained with a modest computational effort. The main procedure of a self-consistent GW calculation consists of four main steps
* Obtain Green's function <math> G^{(j)}_{\bf G G'}({\bf q},i\omega_n) = \left[ \delta_{\bf GG'} - \Sigma^{(j-1)}_{\bf GG'}({\bf q},i\omega_n)G^{(j-1)}_{\bf G G'}({\bf q},i\omega_n)\right]^{-1}G^{(j-1)}_{\bf G G'}({\bf q},i\omega_n)</math>
* Compute irreducible polarizability <math>  \chi^{(j)}({\bf r},{\bf r}',i\tau_m) = -G^{(j)}({\bf r},{\bf r}',i\tau_m)G^{(j)}({\bf r}',{\bf r},-i\tau_m) </math>
* Determine screened potential <math>W^{(j)}_{{\bf G}{\bf G}'}({\bf q},\omega_m)=\left[\delta_{{\bf G}{\bf G}'}-\chi^{(j)}_{{\bf G}{\bf G}'}({\bf q},\omega_m)V_{{\bf G}{\bf G}'}({\bf q})\right]^{-1}V_{{\bf G}{\bf G}'}({\bf q})</math>
* Calculate GW self-energy  <math>  \Sigma^{(j)}({\bf r},{\bf r}',i\tau_m) =- G^{(j)}({\bf r},{\bf r}',i\tau_m)W^{(j)}({\bf r}',{\bf r},i\tau_m) </math>


<span id="QPR">
This procedure can be selected with the following {{TAG|INCAR}} settings
== Low scaling quasi-particle calculations ==
  {{TAGBL|System}} = SiC
If quasi-particle calculations are performed, the self-energy is transformed to reciprocal space in combination with a Pade approximation and the quasi-particle eigenvalue problem is solved.<ref name="liu"/> Here all previously discussed GW approximations are supported and can be selected by adding an additional ''R'' to the {{TAG|ALGO}} option, ''e.g.'' single-shot GW calculations or partially self-consistent quasi-particle energies can be selected via the following lines
  {{TAGBL|System}} = Si
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
  {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE.   
  {{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE.   
  {{TAGBL|NELM}} = number of iterations wanted  
  {{TAGBL|NELMGW}} = number of iterations wanted ! NELM in 6.2 and older
! only quasi-particle energies 
  {{TAGBL|ALGO}} = GW0R ! ALGO = scGW0R has the same effect here, that is self-consistency in G, no update in W
  {{TAGBL|ALGO}} = EVGW0R
The number of self-consistency steps can be set with the {{TAG|NELMGW}} tag.
! or self-consistent quasi-particle energies and orbitals
{{TAGBL|ALGO}} = QPGW0R


Note, that a preceding DFT calculation is not necessary with this setting. The corresponding output is the same as for {{TAG|ALGO}}=''EVGW0'' and ''QPGW0'', but small deviations in the quasi-particle energies and orbitals can be expected due to the analytical continuation. However, the differences should be smaller than 0.04 eV. 
Due to efficiency, VASP performs each step in the Hartree-Fock basis. This is the reason why there are two sets of QP-energies found after the first iteration (one for the QP-energies in the KS-basis and one for the QP energies in the HF basis)
 
After the second iteration, only the QP energies obtained in the HF basis are printed, and a similar output as follows is found in the {{TAG|OUTCAR}} file  
After a successful run, VASP writes the following lines into the {{TAG|OUTCAR}} file


QP shifts evaluated in HF basis
  k-point  1 :      0.0000    0.0000    0.0000
  k-point  1 :      0.0000    0.0000    0.0000
  band No.  KS-energies  sigma(KS)    QP-e(linear)    Z        QP-e(zeros)    Z        occupation    Imag(E_QP)    QP_DIFF  
  band No.  KS-energies  sigma(KS)    QP-e(linear)    Z        QP-e(zeros)    Z        occupation    Imag(E_QP)    QP_DIFF TAG
   
   
       1      -6.4074     -6.8617     -6.7178       0.6832     -6.7168       0.6789       2.0000      -1.1944       0.0058
       1      -7.1626     -8.6510     -8.2275       0.7154     -8.2173       0.7017       2.0000      -1.3177       0.0000  2
       2       5.5804      5.0902      5.2064       0.7629      5.2073       0.7591       2.0000       0.0001       0.0018
       2     -2.0899      -3.4157      -3.0348       0.7127      -3.0269       0.7008       2.0000     -0.5614       0.0000  2
       3       5.5804      5.0894      5.2058       0.7629      5.2067       0.7591       2.0000       0.0001      0.0018
       3     -2.0899      -3.4157      -3.0348       0.7127      -3.0269       0.7008       2.0000     -0.5614      -0.0000  2
       4       5.5804      5.0893      5.2057       0.7629      5.2066       0.7591       2.0000       0.0001       0.0018
       4     -2.0899      -3.4157      -3.0348       0.7127      -3.0269       0.7008       2.0000     -0.5614       0.0000  2
       5      8.1035      8.4039      8.3336       0.7660      8.3332       0.7632       0.0000      -0.0190       0.0031
       5      0.4604      -0.8170      -0.4857       0.7407      -0.4768       0.7266       2.0000      -0.1945       0.0000  2
       6      8.1035      8.4026      8.3326       0.7660      8.3322       0.7632       0.0000      -0.0192      0.0034
       6      0.4604      -0.8170      -0.4857       0.7407      -0.4768       0.7266       2.0000      -0.1945      -0.0000  2
       7      8.1035      8.4023      8.3324       0.7660      8.3320       0.7632       0.0000      -0.0192       0.0035
       7      0.4604      -0.8170      -0.4857       0.7407      -0.4768       0.7266       2.0000      -0.1945       0.0000  2
       8      8.9494       9.4257       9.3101       0.7574       9.3089       0.7524       0.0000      -0.0410       0.0145
       8      5.1013       4.0069       4.2594       0.7693       4.2645       0.7598       2.0000      -0.0602       0.0000  2
Here the meaning of each column is the same as for the other low-scaling GW algorithms. 
</span>


Here column four is obtained by a linearization of the self-energy around the energies from previous step and has the same meaning as the third column in {{TAG|ALGO}}=''EVGW0'' calculations, while column six is obtained from the roots of
<span id="scGWR">
 
== Fully self-consistent GW caluclations: GWR ==
<math>
If the screened potential should be updated during the self-consistency circle {{cite|grumet:prb:2018}} the following {{TAG|INCAR}} file can be used
\langle \phi^{(0)}_{n{\bf q}}| \Sigma(\omega) | \phi^{(0)}_{n{\bf q}}\rangle -\omega =0
{{TAGBL|System}} = SiC
</math>
{{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
{{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE. 
{{TAGBL|NELMGW}} = number of iterations wanted ! use NELM in VASP.6.2 and older
{{TAGBL|ALGO}} = GWR ! ALGO = scGWR has the same effect here, that is self-consistency in G and W
The output is similar to partially self-consistent GW calculations, with the difference that ''KS-energies'' are replaced by the QP energies from previous iteration.
</span>


If self-consistent quasi-particle calculations are performed VASP uses the latter as new eigenvalues.  
=== Caveats ===
Using this option, similar caveats can be expected as for {{TAG|ALGO}}=''EVGW'' and ''QPGW'' calculations and we recommend to leave out the {{TAG|LOPTICS}} and {{TAG|LPEAD}} line for metals.  


Similarly, low scaling self-consistent quasi-particle calculations can be chosen where the screened potential is updated as well
The cubic scaling space-time GW algorithm requires considerably more memory than the corresponding quartic-scaling implementations, two Green's functions <math>G({\bf r,r'},i\tau_n)</math> have to be stored in real-space. To reduce the memory overhead, VASP exploits Fast Fourier Transformations (FFT) to avoid storage of the matrices on the (larger) real space grid, on the one hand. The precision of the FFT can be selected with {{TAG|PRECFOCK}}, where usually the values ''Fast'' sufficient.


{{TAGBL|System}} = Si
On the other hand, the code avoids storage of redundant information, i.e., both the Green's function and polarizability matrices are distributed as well as the individual imaginary grid points. The distribution of the imaginary grid points can be set by hand with the {{TAG|NTAUPAR}} and {{TAG|NOMEGAPAR}} tags, which splits the imaginary grid points {{TAG|NOMEGA}} into {{TAG|NTAUPAR}} time and {{TAG|NOMEGAPAR}} groups. For this purpose both tags have to be divisors of {{TAG|NOMEGA}}.
{{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
{{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE.
{{TAGBL|NELM}} = number of iterations wanted
! eigenvalues only 
{{TAGBL|ALGO}} = EVGWR
! or eigenvalues and orbitals
{{TAGBL|ALGO}} = QPGWR


Here the same caveats can be expected as for {{TAG|ALGO}}=''EVGW'' or ''QPGW''.  
The default values are usually reasonable choices provided the tag {{TAG|MAXMEM}} is set correctly and we strongly recommend to set {{TAG|MAXMEM}} instead of {{TAG|NTAUPAR}}.  
</span>
{{NB|important|As of version 6.2, {{TAGBL|MAXMEM}} is estimated automatically (if not set) from the "MemAvailable" entry of the Linux kernel in "/proc/meminfo".}}


<span id="scGWR">
The required storage for a low-scaling RPA or GW calculation depends mostly on {{TAG|NTAUPAR}}, the number of MPI groups that share same imaginary time points. A rough estimate for the required bytes is given by
== Self-consistent GW caluclations ==
The space-time implementation allows for true self-consistent GW calculations. That is, the solution of the Dyson equation for the Green's function can be obtained with a modest computational effort. The main procedure of a self-consistent GW calculation consists of four main steps
* Obtain Green's function <math> G^{(j)}_{\bf G G'}({\bf q},i\omega_n) = \left[ \delta_{\bf GG'} - \Sigma^{(j-1)}_{\bf GG'}({\bf q},i\omega_n)G^{(j-1)}_{\bf G G'}({\bf q},i\omega_n)\right]^{-1}G^{(j-1)}_{\bf G G'}({\bf q},i\omega_n)</math>
* Compute irreducible polarizability <math>  \chi^{(j)}({\bf r},{\bf r}',i\tau_m) = -G^{(j)}({\bf r},{\bf r}',i\tau_m)G^{(j)}({\bf r}',{\bf r},-i\tau_m) </math>
* Determine screened potential <math>W^{(j)}_{{\bf G}{\bf G}'}({\bf q},\omega_m)=\left[\delta_{{\bf G}{\bf G}'}-\chi^{(j)}_{{\bf G}{\bf G}'}({\bf q},\omega_m)V_{{\bf G}{\bf G}'}({\bf q})\right]^{-1}V_{{\bf G}{\bf G}'}({\bf q})</math>
* Calculate GW self-energy  <math>  \Sigma^{(j)}({\bf r},{\bf r}',i\tau_m) =- G^{(j)}({\bf r},{\bf r}',i\tau_m)W^{(j)}({\bf r}',{\bf r},-i\tau_m) </math>


This procedure can be selected with the following {{TAG|INCAR}} settings
(NGX*NGY*NGZ)*(NGX_S*NGY_S*NGZ_S) / ( NCPU  / {{TAGBL|NTAUPAR}} ) * 16


  {{TAGBL|System}} = Si
where "NCPU" is the number of MPI ranks used for the job,"NGX,NGY,NGZ" denotes the number of FFT grid points for the exact exchange and "NGX_S,NGY_S,NGZ_S" the number of FFT grid points for the supercell. Note, both grids are written to the {{TAGBL|OUTCAR}} file after the lines
{{TAGBL|ISMEAR}} = 0 ; {{TAGBL|SIGMA}} = 0.05
{{TAGBL|LOPTICS}} = .TRUE. ; {{TAGBL|LPEAD}} = .TRUE. 
{{TAGBL|NELM}} = number of iterations wanted
! eigenvalues only 
{{TAGBL|ALGO}} = GW0R


The number of self-consistency steps can be set with the {{TAG|NELM}} tag and usually 4 cycles suffice to reach self-consistency.
FFT grid for exact exchange (Hartree Fock)
After each cycle VASP writes the following lines into the {{TAG|OUTCAR}} file
FFT grid for supercell:


band No. KS-energies  sigma(KS)    QP-e(linear)    Z        QP-e(zeros)    Z        occupation    Imag(E_QP)    QP_DIFF
The smaller {{TAG|NTAUPAR}} is set, the less memory per node the job requires to finish successfully.  
      1    -11.6864    -15.1191    -14.6792      0.8719    -14.0937      0.9442      2.0000      -2.7943      0.6135
      2      9.8536      9.1173      9.2220      0.8577      9.2221      0.8575      2.0000      0.0347      0.1150
      3      9.8536      9.1175      9.2222      0.8577      9.2223      0.8575      2.0000      0.0347      0.1150
      4      9.8536      9.1182      9.2228      0.8577      9.2229      0.8575      2.0000      0.0347      0.1150
      5      15.4136      15.8263      15.7690      0.8613      15.7689      0.8610      0.0000      0.0045      0.1762
      6      15.4136      15.8265      15.7692      0.8613      15.7691      0.8610      0.0000      0.0045      0.1762
      7      15.4136      15.8267      15.7694      0.8613      15.7694      0.8610      0.0000      0.0045      0.1762


The first column is the band index, the second column the quasi-particle energies from the previous cycle (DFT energies in the first cycle), the third column the diagonal matrix elements of the self-energy and  
The approximate memory requirement is calculated in advance and printed to screen and {{TAGBL|OUTCAR}} as follows:


</span>
min. memory requirement per mpi rank 1234 MB, per node 9872 MB


== Related Tags and Sections ==
= Related tags and articles =


* {{TAG|ALGO}} for response functions and ''GW'' calculations
* {{TAG|ALGO}} for response functions and ''GW'' calculations
* {{TAG|LMAXFOCKAE}}
* {{TAG|LOPTICS}}, derivative of wavefunction w.r.t. <math>k</math>
* {{TAG|LPEAD}}, derivative of wavefunction with finite differences
* {{TAG|LMAXFOCKAE}} overlap densities and multipoles
* {{TAG|MAXMEM}}, memory available to one mpi rank on each node
* {{TAG|NOMEGA}}, {{TAG|NOMEGAR}} number of frequency points
* {{TAG|NOMEGA}}, {{TAG|NOMEGAR}} number of frequency points
* {{TAG|LSPECTRAL}}: use the spectral method for the polarizability
* {{TAG|LSPECTRAL}}, use the spectral method for the polarizability
* {{TAG|LSPECTRALGW}}: use the spectral method for the self-energy
* {{TAG|LSPECTRALGW}}, use the spectral method for the self-energy
* {{TAG|OMEGAMAX}}, {{TAG|OMEGATL}} and {{TAG|CSHIFT}}
* {{TAG|OMEGAMAX}}, {{TAG|OMEGATL}} and {{TAG|CSHIFT}}
* {{TAG|ENCUTGW}}: energy cutoff for response function
* {{TAG|ENCUTGW}}, energy cutoff for response function
* {{TAG|ENCUTGWSOFT}}: soft cutoff for Coulomb kernel
* {{TAG|ENCUTGWSOFT}}, soft cutoff for Coulomb kernel
* {{TAG|ODDONLYGW}} and {{TAG|EVENONLYGW}}: reducing the ''k''-grid for the response functions
* {{TAG|ODDONLYGW}} and {{TAG|EVENONLYGW}}, reducing the ''k''-grid for the response functions
* {{TAG|LSELFENERGY}}: the frequency dependent self energy
* {{TAG|LSELFENERGY}}, the frequency dependent self energy
* {{TAG|LWAVE}}: selfconsistent ''GW''
* {{TAG|LWAVE}}, self-consistent ''GW''
* {{TAG|NOMEGAPAR}}: frequency grid parallelization
* {{TAG|NOMEGAPAR}}, frequency grid parallelization
* {{TAG|NTAUPAR}}: time grid parallelization
* {{TAG|NTAUPAR}}, time grid parallelization
* {{TAG|NATURALO}}, natural orbitals
* {{TAG|LALL_IN_ONE}}, all-in-one ''GW'' mode
* {{TAG|IALL_IN_ONE}}, all-in-one ''GW'' mode
* {{TAG|NBANDSEXACT}}, number of KS bands in all-in-one mode
* {{TAG|NBANDS_WAVE}}, number of bands written to {{FILE|WAVECAR}} in all-in-one mode
* {{TAG|LSINGLES}}, singles contribution to correlation energy
{{sc|GW calculations|Examples|Examples that use this tag}}


{{sc|GW calculations|Examples|Examples that use this tag}}
= References =


== References ==
<references>
<ref name="Hedin">[https://journals.aps.org/pr/abstract/10.1103/PhysRev.139.A796 L. Hedin, Phys. Rev. 139, A796 (1965)]</ref>
<ref name="shishkin-PRB74">[http://link.aps.org/doi/10.1103/PhysRevB.74.035101 M. Shishkin and G. Kresse, Phys. Rev. B 74, 035101 (2006).]</ref>
<ref name="shishkin-PRB75">[http://link.aps.org/doi/10.1103/PhysRevB.75.235102 M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 (2007).]</ref>
<ref name="shishkin:prl:07">[http://link.aps.org/doi/10.1103/PhysRevLett.99.246403 M. Shishkin, M. Marsman, and G. Kresse, Phys. Rev. Lett. 99, 246403 (2007).]</ref>
<ref name="fuchs-PRB76">[http://link.aps.org/doi/10.1103/PhysRevB.76.115109 F. Fuchs, J. Furthmüller, F. Bechstedt, M. Shishkin, and G. Kresse, Phys. Rev. B 76, 115109 (2007).]</ref>
<ref name="liu">[http://journals.aps.org/prb/abstract/10.1103/PhysRevB.94.165109 P. Liu, M. Kaltak, J. Klimes and G. Kresse, Phys. Rev. B 94, 165109 (2016).]</ref>
<ref name="rojas">[https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.74.1827 H. N. Rojas, R. W. Godby, R. J. Needs, Phys. Rev. Lett. 74, 1827 (1995)]</ref>
<ref name="HybertsenLouie">[https://journals.aps.org/prb/abstract/10.1103/PhysRevB.34.5390 M. S. Hybertsen, S. G. Louie Phys. Ref. B 34, 5390 (1986)]</ref>
<ref name="bruneval:prb:06">[http://link.aps.org/doi/PhysRevB.74.45102 F. Bruneval, N. Vast, and L. Reining, Phys. Rev. B 74, 45102 (2006).]</ref>
</references>
----
[[The_VASP_Manual|Contents]]


[[Category:INCAR]][[Category:GW]][[Category:Howto]]
[[Category:Many-body perturbation theory]][[Category:GW]][[Category:Howto]][[Category:Low-scaling GW and RPA]]

Latest revision as of 14:21, 18 October 2024

The GW approximation is an approximation to the self-energy. GW calculations are available as of VASP.5.X. For details on the implementation and use of the GW routines, we recommend the papers by Shishkin et al. [1][2][3] and Fuchs et al.[4]

Single step procedure: GW in one go

As of VASP.6.3 all GW approximations can be run in one single run by selecting the corresponding ALGO tag and omitting NBANDS), for instance like so

System = SiC
ALGO = EVGW0, QPGW0, EVGW, QPGW, GW0R or GWR # use an algorithgm described below
NELMGW = 1,2,.. # number of self-consistency cycles
ISMEAR = 0 ; SIGMA = 0.05  ! small sigma is required to avoid partial occupancies
LOPTICS = .TRUE.  # for insulators, omit for metals

Note, NBANDS must not be present in the INCAR to select this procedure.

Important: In older versions a two step procedure is required, where the first step is always a DFT calculation and the second step the actual GW calculation.

The two-step procedure is described below.

Caveats

The single-step GW procedure performs a DFT step internally with an exact diagonalization of the Kohn-Sham Hamiltonian using the maximum available NBANDS supported for the chosen ENCUT value. Consequently, a large number of unoccupied bands is initialized with random plane-wave coefficients. In rare cases, this yields two linearly dependent column vectors in the Hamiltonian and results in LAPACK errors like "ZPOTRF fails". These errors can be prevented using the two-step GW procedure as described below. Furthermore, one can "ramp up" NBANDS to the maximum value by repeatedly restarting the DFT calculation from a pre-converged WAVECAR with fewer bands.

First step: DFT calculation

GW calculations always require a one-electron basis set. Usually this set is obtained from a standard DFT calculation and written into the WAVECAR file and can be calculated for instance the following INCAR file:

System = SiC
NBANDS = 512
ISMEAR = 0 ; SIGMA = 0.05  ! small sigma is required to avoid partial occupancies
LOPTICS = .TRUE.

Note, that a significant number of empty bands is required for GW calculations, so that it might be better to perform the calculations in two steps: first a standard ground-state calculation with few unoccupied orbitals only,

System = SiC ground-state occupied orbitals
ISMEAR = 0 ; SIGMA = 0.05  ! small sigma is required to avoid partial occupancies
EDIFF = 1E-8               ! required tight tolerance for ground-state orbitals

and, second, a calculation of a large number of unoccupied orbitals

System  = SiC unoccupied orbitals
ALGO = Exact               ! use exact diagonalization of the Hamiltonian
NELM = 1                   ! since we are already converged stop after one step
NBANDS = 512               ! maybe even larger                
ISMEAR = 0 ; SIGMA = 0.05  ! small sigma is required to avoid partial occupancies
LOPTICS = .TRUE.

Furthermore, note that the flag LOPTICS=.TRUE. is required to write the file WAVEDER, which contains the derivative of the orbitals with respect to k. This derivative is used to construct the head and wings of the dielectric matrix employing k·p perturbation theory and is important to accelerate k-point convergence for insulators and semiconductors.

Warning: For metals, in general, we recommend omitting the LOPTICS tag and removing the WAVEDER file from the directory.

Optional: Use Hybrid functionals

Optionally, one can start a GW calculation from a hybrid functional, such as HSE. For hybrid functionals, the two step procedure will accordingly involve the following INCAR files. In the first step, converged HSE03 orbitals are determined (see here for a selection of available hybrid functionals):

System  = SiC ground-state occupied orbitals
ISMEAR = 0 ; SIGMA = 0.05
ALGO = Damped ; TIME = 0.5  ! or ALGO = Conjugate
LHFCALC = .TRUE. ; AEXX = 0.25 ; HFSCREEN = 0.3 
EDIFF = 1E-6      ! required tight tolerance for ground-state orbitals

Secondly, determine the HSE03 orbitals for unoccupied states:

System  = SiC unoccupied orbitals
NBANDS = 512      ! maybe even larger
ALGO = Exact
NELM = 1          ! since we are already converged stop after one step
ISMEAR = 0 ; SIGMA = 0.05
LHFCALC = .TRUE. ; AEXX = 0.25 ; HFSCREEN = 0.3 
LOPTICS = .TRUE. # for insulators

Second step: GW calculation

The actual GW calculation is done in a second step. Here different GW flavors are possible and are selected with the ALGO tag.

Note that as of VASP.6 the GW ALGO tags have been renamed, see here for VASP.5.X tags.

Single shot quasiparticle energies: G0W0

This is the simplest GW calculation and computationally the most efficient one. A single-shot calculation is often referred to as G0W0 and calculates the quasiparticle energies from a single GW iteration by neglecting all off-diagonal matrix elements of the self-energy and employing a Taylor expansion of the self-energy around the DFT energies . The corresponding equation becomes

with the renormalization factor

In VASP, G0W0 calculations are selected using an INCAR file such as

System  = SiC
NBANDS = 512
ISMEAR = 0 ; SIGMA = 0.05
NELMGW = 1 ! use NELM for VASP.6.2 and older 
ALGO = EVGW0 ! use "GW0" for VASP.5.X
NOMEGA = 50
Mind: Convergence with respect to the number of empty bands NBANDS and with respect to the number of frequencies NOMEGA must be checked carefully.

To avoid complicated inter-nested tests, we recommend calculating all orbitals that the plane-wave basis set allows to calculate (except for simple tests). For further reading, please consult the section on ENCUTGW.

After a successful G0W0 run, VASP will write the quasiparticle energies into the OUTCAR file for a set of NBANDSGW bands for every k-point in the Brillouin zone. Look for lines similar to

QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
for sc-GW calculations column KS-energies equals QP-energies in previous step 
and V_xc(KS)=  KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1  + <V_x>^1)

k-point   1 :       0.0000    0.0000    0.0000
 band No.  KS-energies  QP-energies   sigma(KS)   V_xc(KS)     V^pw_x(r,r')   Z            occupation Imag(sigma)

     1      -7.1627      -8.3040     -14.5626     -12.7276     -21.6682       0.6219       2.0000       1.2037
     2      -2.0901      -3.4347     -15.7660     -14.2799     -21.7439       0.9048       2.0000       0.6914
     3      -2.0901      -3.4347     -15.7660     -14.2799     -21.7439       0.9048       2.0000       0.6914
     4      -2.0901      -3.4347     -15.7660     -14.2799     -21.7439       0.9048       2.0000       0.6914
     5       0.4603      -0.4663     -13.7603     -12.5200     -18.1532       0.7471       2.0000       0.2167
     6       0.4603      -0.4663     -13.7603     -12.5200     -18.1532       0.7471       2.0000       0.2167

The first column is the band index and the third column denotes the quasiparticle energies . Column two, four, five and seven refer to the DFT energies , diagonal matrix elements of the self-energy , the exchange-correlation potential and the renormalization factor defined above, respectively.

Partially self-consistent calculations: EVGW0

In most cases, the best results (i.e., closest to experiment) are obtained by iterating only via the spectral representation

but keeping and the orbitals fixed to the initial DFT level. This method goes back to Hybertsen and Louie [5] and can be achieved in two ways.

If the spectral method is not selected (LSPECTRAL=.FALSE., requiring much more compute time), the quasiparticle (QP) shifts are iterated automatically four times, and one finds four sets of QP shifts in the OUTCAR file. The first one corresponds to the G0W0 case. The INCAR file is simply:

System = SiC
NBANDS = 512
ISMEAR = 0 ; SIGMA = 0.05
ALGO = EVGW0 ! use "GW0" in VASP.5.X
LSPECTRAL =.FALSE.
Tip: In self-consistent GW calculations, convergence with the number of updated bands NBANDSGW must be checked carefully.

For technical reasons, it is not possible to iterate in this manner if LSPECTRAL=.TRUE. is set in the INCAR file (this is the default). In this case, an iteration number must be supplied in the INCAR file using the NELMGW tag. Usually, three to four iterations are sufficient to obtain accurate QP shifts.

System = SiC
NBANDS = 512
ISMEAR = 0 ; SIGMA = 0.05
ALGO = EVGW0 ! use "GW0" in VASP.5.X
NELMGW = 4 ! use NELM in VASP.6.2 and older

The results are found again in the OUTCAR file

QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 4

k-point   1 :       0.0000    0.0000    0.0000
 band No. old QP-enery  QP-energies   sigma(KS)   T+V_ion+V_H  V^pw_x(r,r')   Z            occupation Imag(sigma)

     1      -8.6924      -8.7107     -14.2871       5.5647     -21.6681       0.6076       2.0000       1.1648
     2      -3.4692      -3.4806     -15.6742      12.1894     -21.7437       0.7304       2.0000       0.6351
     3      -3.4692      -3.4806     -15.6742      12.1894     -21.7437       0.7304       2.0000       0.6351
     4      -3.4692      -3.4806     -15.6742      12.1894     -21.7437       0.7304       2.0000       0.6351
     5      -0.6957      -0.7006     -13.6827      12.9802     -18.1531       0.7264       2.0000       0.2769
     6      -0.6957      -0.7006     -13.6827      12.9802     -18.1531       0.7264       2.0000       0.2769
     7      -0.6957      -0.7006     -13.6827      12.9802     -18.1531       0.7264       2.0000       0.2769

In contrast to single shot GW calculations, the second column represent now the QP-energies from the previous iteration.

Partially self-consistent quasiparticle calculations: QPGW0

If non diagonal components of the self-energy (in the orbital basis) should be included use ALGO=QPGW0. The following setting can be used:

System = SiC
NBANDS = 512
ISMEAR = 0 ; SIGMA = 0.05
ALGO = QPGW0 ! or "scGW0" for VASP.5.2.11 and older 
LOPTICS = .TRUE. ; LPEAD = .TRUE. ! ommit this lines for metals
NELMGW = 4 ! use NELM for VASP.6.2 and older

In this case, the orbitals are updated as well by constructing a hermitian (energy independent) approximation to the self-energy [3]. The "static" COHSEX approximation can be selected by setting NOMEGA = 1 [6]. To improve convergence to the ground-state, the charge density (and the charge density only) is mixed using a Kerker type mixing in VASP.5.3.2 and more recent versions (see IMIX). The mixing parameters AMIX, BMIX, AMIX_MAG, BMIX_MAG, AMIN can be adjusted, if convergence problems are encountered.

We strongly urge the user to monitor convergence by inspecting the lines

charge density residual

in the OUTCAR files.

Alternatively, the mixing may be switched off by setting IMIX=0 and controlling the step width for the orbitals using the parameter TIME (which defaults to 0.4). This selects a fairly sophisticated damped MD algorithm that is also used for DFT methods when ALGO=Damped. This method is generally more reliable for metals and materials with strong charge sloshing.

After every iteration, VASP writes the following lines into the OUTCAR file

QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
    GWSYM:  cpu time   15.8978: real time   15.9528

k-point   1 :       0.0000    0.0000    0.0000
 band No. DFT-energies  QP-energies  QP-e(diag)   sigma(DFT)    Z            occupation

     1      -7.1626      -8.4217      -8.3038      -8.9978       0.6219       2.0000
     2      -2.0899      -3.4394      -3.4347      -3.5765       0.9047       2.0000
     3      -2.0899      -3.4394      -3.4347      -3.5765       0.9047       2.0000
     4      -2.0899      -3.4394      -3.4347      -3.5765       0.9047       2.0000
     5       0.4604      -0.4787      -0.4663      -0.7800       0.7471       2.0000
     6       0.4604      -0.4787      -0.4663      -0.7800       0.7471       2.0000
     7       0.4604      -0.4787      -0.4663      -0.7800       0.7471       2.0000
     8       5.1013       4.1883       4.2149       3.9518       0.7711       2.0000

For the first iteration, here, the fourth column should be identical to the third column of the G0W0 results discussed above. The third column reports the quasiparticle energies obtained from including the off-diagonal matrix elements in the eigenvalue equation.

Caveats

The QPGW0 (or scGW0 in VASP.5.2.11 and older) must be used with great caution, particularly in combination with symmetry. Symmetry is handled in a rather sophisticated manner. Specifically, only the minimal number of required combinations of q and k points is considered. In this case, symmetry must be applied to restore the full star of q. This is done by determining degenerate eigenvalue/eigenvector pairs and restoring their symmetry according to their irreducible representation. Although the procedure is generally relatively reliable, it fails to work properly if the degenerate states do not possess eigenvalues that are sufficiently close due to insufficient convergence in the preceding DFT calculations. That is because states are treated as degenerate if, and only if, their eigenenergies are within 0.01 eV.

For large supercells with low symmetry, we strongly recommend switching off symmetry.

Self-consistent EVGW and QPGW calculations

Self-consistent QPGW calculations are only supported in a QP picture. As for QPGW0, it is possible to update the eigenvalues only (ALGO=EVGW or GW for VASP.5.X), or the eigenvalues and one-electron orbitals (ALGO=QPGW or scGW in VASP.5.2.11 and older). In all cases, a QP picture is maintained, i.e., satellite peaks (shake ups and shake downs) can not be accounted for in the self-consistency cycle. Self-consistent QPGW calculations can be performed by simply repeatedly calling VASP using:

System = SiC
NBANDS = 512
ISMEAR = 0 ; SIGMA = 0.05
ALGO = EVGW    ! "GW" in VASP.5.X, eigenvalues only  or alternatively
ALGO = QPGW    ! "scGW" in VASP.5.2.11 and older, eigenvalues and one electron orbitals

For QPGW0 or QPGW, nondiagonal terms in the Hamiltonian are accounted for, e.g. the linearized QP equation is diagonalized, and the one-electron orbitals are updated [3]. Alternatively (and preferably), the user can specify an electronic iteration counter using NELMGW (NELM in VASP.6.2 and older):

System = SiC
NBANDS = 512
ISMEAR = 0 ; SIGMA = 0.05
NELMGW = 3   ! use NELM in VASP.6.2 and older
ALGO = EVGW  ! "GW" in VASP.5.X 
# or  
ALGO = QPGW  ! "scGW" in VASP.5.2.11 and older

In this case, the one-electron energies (=QP energies) are updated 3 times (starting from the DFT eigenvalues) in both G and W. For ALGO=QPGW (or ALGO=scGW in VASP.5.2.11 and older), the one electron energies and one electron orbitals are updated 3 times [3]. As for ALGO = QPGW0 (or scGW0 in vasp.5.2.11 and older), the "static" COHSEX approximation can be selected by setting NOMEGA=1 [6].

To improve convergence to the ground-state, the charge density is mixed using a Kerker type mixing starting with VASP.5.3.2 (see IMIX). The mixing parameters AMIX, BMIX, AMIX_MAG, BMIX_MAG, AMIN can be adjusted, if convergence problems are encountered. Alternatively, the mixing may be switched off by setting IMIX=0 and controlling the step width for the orbitals using the parameter TIME (which defaults to 0.4). This selects a fairly sophisticated damped MD algorithm that is also used for DFT methods when ALGO=Damped. This method is generally more reliable for metals and materials with strong charge sloshing.

Additional information about this method is found here.

Caveats

Fully self-consistent QPGW calculations with an update of the orbitals in and [3] require significant care and are prone to diverge (QPGW0 calculations are usually less critical). As discussed, above, one can select this mode using:

System = SiC
NBANDS = 512
ISMEAR = 0 ; SIGMA = 0.05
ALGO = QPGW  ! or "scGW" in VASP.5.2.11 and older, eigenvalues and one-electron orbitals
NELMGW = number of steps ! use NELM for VASP.6.2 and older

However, one caveat applies to this case: when the orbitals are updated, the derivatives of the orbitals with respect to (stored in the WAVEDER file) will become incompatible with the orbitals. This can cause severe problems and convergence to the incorrect solution.

Warning: For metals, in general, we recommend omitting the LOPTICS tag and removing the WAVEDER file from the directory.

For insulators, VASP (version 5.3.2 or higher) can update the WAVEDER file in each electronic iteration if the finite difference method is used to calculate the first derivative of the orbitals with respect to :

System = SiC
NBANDS = 512
ISMEAR = 0 ; SIGMA = 0.05
ALGO = QPGW ! "scGW" in VASP.5.2.11 and older, eigenvalues and one-electron orbitals
NELMGW  = 10 ! use NELM in VASP.6.2 and older
LOPTICS = .TRUE. ; LPEAD = .TRUE.

The combination LOPTICS=.TRUE.; LPEAD=.TRUE. is required since is not available for GW like methods. LPEAD=.TRUE. circumvents this problem by calculating the derivatives of the orbitals using numerical differentiation on the finite k-point grid (this option is presently limited to insulators).

Vertex corrections are presently not documented. This is a feature still under construction, and we recommend collaborating with the Vienna group if you desperately need that feature.

Low scaling GW algorithms

The GW implementations in VASP described in the papers of Shishkin et al.[1] [2] avoids storage of the Green's function as well as Fourier transformations between time and frequency domain entirely. That is, all calculations are performed solely on the real frequency axis using Kramers-Kronig transformations for convolutions in the equation of and in reciprocal space.

As of VASP.6 a new cubic scaling GW algorithm [7] (called space-time implementation in the following) can be selected. This approach follows the idea of Rojas et al. [8] and performs the GW self-consistency cycle on imaginary time and imaginary frequency axes .

Tip: Using the low-scaling GW algorithm also calculates the total energy in the Random Phase approximation (RPA), which is described in a separate article.

Low scaling, single shot GW calculations: G0W0R

The low-scaling analogue of G0W0 is selected with ALGO=G0W0R. In contrast to the single-shot GW calculations on the real-axes, here the self-energy is determined on the imaginary frequency axis. To this end, the overall scaling is reduced by one order of magnitude and is cubic with respect to the system size, because a small value for NOMEGA can be used (usually <20).

This algorithm evaluates:

  • Single-shot GW quasiparticle energies (from an analytical continuation of the self-energy to the real axis)[7]
  • Natural orbitals from the first order change of the density matrix (i.e. ), see the NATURALO tag for more information [9].
Mind: This selection ignores NELMGW.

Following INCAR file selects the low-scaling GW algorithm:

System = SiC
ISMEAR = 0 ; SIGMA = 0.05
LOPTICS = .TRUE.  
ALGO = G0W0R
NOMEGA = 12 ! small number of frequencies necessary

Search the OUTCAR file for the following lines

  QP shifts evaluated in KS or natural orbital/ Bruckner basis
  k-point   1 :       0.0000    0.0000    0.0000
  band No.  KS-energies   sigma(KS)    QP-e(linear)    Z         QP-e(zeros)     Z        occupation    Imag(E_QP)    QP_DIFF TAG
 
       1      -7.1627      -8.6732      -8.2451       0.7166      -8.2346       0.7026       2.0000      -1.3101       0.0000   2
       2      -2.0901      -3.4155      -3.0350       0.7129      -3.0272       0.7011       2.0000      -0.5582      -0.0000   2
       3      -2.0901      -3.4155      -3.0350       0.7129      -3.0272       0.7011       2.0000      -0.5582       0.0000   2
       4      -2.0901      -3.4155      -3.0350       0.7129      -3.0272       0.7011       2.0000      -0.5582      -0.0000   2
       5       0.4603      -0.8219      -0.4904       0.7414      -0.4814       0.7273       2.0000      -0.1902       0.0000   2
       6       0.4603      -0.8219      -0.4904       0.7414      -0.4814       0.7273       2.0000      -0.1902      -0.0000   2

Here column four is obtained by a linearization of the self-energy around the Kohn-Sham energies (second column) and can be compared to the third column of single-shot GW calculations on the real axis. Column six represents another set of QP-energies that is obtained from the roots of the following equation

These roots represent the poles of the Green's function in the spectral representation.

Output description

The meaning of each column is explained briefly in the following.

  • band No. the band index of KS orbital at given k-point
  • KS-energies eigenenergies corresponding to band index
  • sigma(KS) diagonal matrix elements of self-energy evaluated at KS energies
  • QP-e(linear) quasiparticle energies obtained from linearizing frequency dependence of diagonal self-energy around KS energies
  • Z renormalization factor obtained from five-point stencil for derivative of self-energy w.r.t. frequency
  • QP-e(zeros) quasiparticle energies obtained from full frequency dependence of self-energy, i.e. real part of complex pole of Green's function
  • Z renormalization factor obtained from central difference for derivative of self-energy w.r.t. frequency
  • occupation occupation number for band at given k-point
  • Imag(E_QP) imaginary part of complex pole , i.e. measure for inverse lifetime of quasi-particle
  • QP_DIFF difference of QP energies (of linearized self-energy) obtained from Eq. 77 of Liu et. al.[7] and M. Grumets thesis[10].

Optional: RPA Forces

Optionally, RPA forces can be calculated by adding following line to the INCAR:

 LRPAFORCE = .TRUE. 

After the QP-energies, VASP performs a linear-response calculation that is required for the RPA forces.[11] Following data block in the OUTCAR file can be found after a successful run:

POSITION                                       TOTAL RPA FORCE (eV/Angst)
-----------------------------------------------------------------------------------
     0.17542     -0.22348      0.17542        -0.292069      7.581315     -0.292069
     1.12850      1.31044      1.12850         0.304683     -7.605527      0.304683
-----------------------------------------------------------------------------------
   total drift:                                0.012614     -0.024212      0.012614

SUGGESTED UPDATED POSCAR (direct coordinates)  step
-----------------------------------------------------------------------------------
 -0.00958461  -0.00958461   0.13485779         0.04179056   0.04179056   0.00283088
  0.25787833   0.25787833   0.22191754        -0.04337198  -0.04337198   0.00431513
Warning: Currently RPA forces for metallic systems are not supported.

Low scaling, partially self-consistent GW calculations: EVGW0R

Mind: available as of vasp 6.4.0

The low-scaling analogue of EVGW0 is selected with ALGO=EVGW0R. Following INCAR file selects this algorithm:

System = SiC
ISMEAR = 0 ; SIGMA = 0.05
LOPTICS = .TRUE.  
ALGO = EVGW0R 
NELMGW = 4  ! number of iterations in G
NOMEGA = 12 ! small number of frequencies necessary

After each iteration, a similar block of data as for ALGO=G0W0R calculations is written to OUTCAR showing the NBANDSGW updated quasi-particle energies (poles) of the Green's function.

Partially self-consistent GW calculations: GW0R

The space-time implementation allows for true self-consistent GW calculations. That is, the solution of the Dyson equation for the Green's function can be obtained with a modest computational effort. The main procedure of a self-consistent GW calculation consists of four main steps

  • Obtain Green's function
  • Compute irreducible polarizability
  • Determine screened potential
  • Calculate GW self-energy

This procedure can be selected with the following INCAR settings

System = SiC
ISMEAR = 0 ; SIGMA = 0.05
LOPTICS = .TRUE. ; LPEAD = .TRUE.  
NELMGW = number of iterations wanted ! NELM in 6.2 and older
ALGO = GW0R ! ALGO = scGW0R has the same effect here, that is self-consistency in G, no update in W

The number of self-consistency steps can be set with the NELMGW tag.

Due to efficiency, VASP performs each step in the Hartree-Fock basis. This is the reason why there are two sets of QP-energies found after the first iteration (one for the QP-energies in the KS-basis and one for the QP energies in the HF basis) After the second iteration, only the QP energies obtained in the HF basis are printed, and a similar output as follows is found in the OUTCAR file

QP shifts evaluated in HF basis
k-point   1 :       0.0000    0.0000    0.0000
band No.  KS-energies   sigma(KS)    QP-e(linear)    Z         QP-e(zeros)     Z        occupation    Imag(E_QP)    QP_DIFF TAG

     1      -7.1626      -8.6510      -8.2275       0.7154      -8.2173       0.7017       2.0000      -1.3177       0.0000   2
     2      -2.0899      -3.4157      -3.0348       0.7127      -3.0269       0.7008       2.0000      -0.5614       0.0000   2
     3      -2.0899      -3.4157      -3.0348       0.7127      -3.0269       0.7008       2.0000      -0.5614      -0.0000   2
     4      -2.0899      -3.4157      -3.0348       0.7127      -3.0269       0.7008       2.0000      -0.5614       0.0000   2
     5       0.4604      -0.8170      -0.4857       0.7407      -0.4768       0.7266       2.0000      -0.1945       0.0000   2
     6       0.4604      -0.8170      -0.4857       0.7407      -0.4768       0.7266       2.0000      -0.1945      -0.0000   2
     7       0.4604      -0.8170      -0.4857       0.7407      -0.4768       0.7266       2.0000      -0.1945       0.0000   2
     8       5.1013       4.0069       4.2594       0.7693       4.2645       0.7598       2.0000      -0.0602       0.0000   2

Here the meaning of each column is the same as for the other low-scaling GW algorithms.

Fully self-consistent GW caluclations: GWR

If the screened potential should be updated during the self-consistency circle [12] the following INCAR file can be used

System = SiC
ISMEAR = 0 ; SIGMA = 0.05
LOPTICS = .TRUE. ; LPEAD = .TRUE.  
NELMGW = number of iterations wanted ! use NELM in VASP.6.2 and older
ALGO = GWR ! ALGO = scGWR has the same effect here, that is self-consistency in G and W

The output is similar to partially self-consistent GW calculations, with the difference that KS-energies are replaced by the QP energies from previous iteration.

Caveats

Using this option, similar caveats can be expected as for ALGO=EVGW and QPGW calculations and we recommend to leave out the LOPTICS and LPEAD line for metals.

The cubic scaling space-time GW algorithm requires considerably more memory than the corresponding quartic-scaling implementations, two Green's functions have to be stored in real-space. To reduce the memory overhead, VASP exploits Fast Fourier Transformations (FFT) to avoid storage of the matrices on the (larger) real space grid, on the one hand. The precision of the FFT can be selected with PRECFOCK, where usually the values Fast sufficient.

On the other hand, the code avoids storage of redundant information, i.e., both the Green's function and polarizability matrices are distributed as well as the individual imaginary grid points. The distribution of the imaginary grid points can be set by hand with the NTAUPAR and NOMEGAPAR tags, which splits the imaginary grid points NOMEGA into NTAUPAR time and NOMEGAPAR groups. For this purpose both tags have to be divisors of NOMEGA.

The default values are usually reasonable choices provided the tag MAXMEM is set correctly and we strongly recommend to set MAXMEM instead of NTAUPAR.

Important: As of version 6.2, MAXMEM is estimated automatically (if not set) from the "MemAvailable" entry of the Linux kernel in "/proc/meminfo".

The required storage for a low-scaling RPA or GW calculation depends mostly on NTAUPAR, the number of MPI groups that share same imaginary time points. A rough estimate for the required bytes is given by

(NGX*NGY*NGZ)*(NGX_S*NGY_S*NGZ_S) / ( NCPU  / NTAUPAR ) * 16

where "NCPU" is the number of MPI ranks used for the job,"NGX,NGY,NGZ" denotes the number of FFT grid points for the exact exchange and "NGX_S,NGY_S,NGZ_S" the number of FFT grid points for the supercell. Note, both grids are written to the OUTCAR file after the lines

FFT grid for exact exchange (Hartree Fock)
FFT grid for supercell:

The smaller NTAUPAR is set, the less memory per node the job requires to finish successfully.

The approximate memory requirement is calculated in advance and printed to screen and OUTCAR as follows:

min. memory requirement per mpi rank 1234 MB, per node 9872 MB

Related tags and articles

Examples that use this tag

References