Practical guide to GW calculations

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. and Fuchs et al. The GW approach gained popularity in the community after Hybertsen and Louie showed that the GW band gaps and spetral properties yield very good agreement with measurements for semiconductors and insulators. This method can be understood in terms of the following eigenvalue equation

$(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}})$

Here $T$ is the kinetic energy, $V_{ext}$ the external potential of the nuclei, $V_{h}$ the Hartree potential and $E_{n{\bf {k}}}$ the quasiparticle energies with orbitals $\phi _{n{\bf {k}}}$. In contrast to DFT, the exchange-correlation potential is replaced by the many-body self-energy $\Sigma$ and should be obtained together with the Green's function $G$, the irreducible polarizability $\chi$, the screened Coulomb interaction $W$ and the irreducible vertex function $\Gamma$ in a self-consistent procedure. For completeness, these equations are

$G(1,2)=G_{0}(1,2)+\int d(3,4)G_{0}(1,3)\Sigma (3,4)G(4,2)$

$\chi (1,2)=\int d(3,4)G(1,3)G(4,1)\Gamma (3,4;2)$

$W(1,2)=V(1,2)+\int d(3,4)V(1,3)\chi (3,4)W(4,2)$

$\Sigma (1,2)=\int d(3,4)G(1,3)\Gamma (3,2;4)W(4,1)$

$\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)$

and provide an alternative to the Schrödinger equation for the many-body problem. Here the common notation $1=({\bf {r}}_{1},t_{1})$ was adopted and $V$ denotes the bare Coulomb interaction.

The GW approximation is obtained by neglecting the equation for the vertex function and using the bare vertex $\Gamma (1,2;3)=\delta (1,2)\delta (1,3)$ instead reducing the equations for the polarizability and self-energy to

$\chi (1,2)=G(1,2)G(2,1)$

$\Sigma (1,2)=G(1,2)W(2,1)$

The equations for the Green's function and the screened potential, however, remain the same, but usually are solved in reciprocal space

$W_{{\bf {G}}{\bf {G}}'}({\bf {q}},\omega )=(\delta _{{\bf {G}}{\bf {G}}'}-\chi _{{\bf {G}}{\bf {G}}'}({\bf {q}},\omega )V_{{\bf {G}}{\bf {G}}'}({\bf {q}}))^{-1}V_{{\bf {G}}{\bf {G}}'}({\bf {q}})$

$G_{{\bf {G}}{\bf {G}}'}({\bf {q}},\omega )=(\delta _{{\bf {G}}{\bf {G}}'}-\Sigma _{{\bf {G}}{\bf {G}}'}({\bf {q}},\omega )G_{{\bf {G}}{\bf {G}}'}^{0}({\bf {q}}))^{-1}G_{{\bf {G}}{\bf {G}}'}^{0}({\bf {q}})$

Due to the complexity further approximations are usually performed giving rise to different GW flavors.

General outline of a GW calculation

As of VASP.6 all GW approximations can be selected directly via the ALGO tag (omitting 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 WAVECAR file in an initial step, using for instance the following INCAR file:

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

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,

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

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

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

Furthermore note that the flag LOPTICS=.TRUE. is required in order 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. For metals, in general, we recommend to omit the LOPTICS tag and remove 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 (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):

System  = Si groundstate 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 groundstate orbitals

Secondly, determine the HSE03 orbitals for unoccupied states:

System  = Si unoccupied orbitals
NBANDS = 96
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.

Step 2: quasi-particle GW

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

G₀W₀: single shot quasi-particle energies

This is the simplest GW calculation, but computationally the most efficient one. A single-shot G₀W₀ calculation calculates the quasi-particle energies from a single GW iteration by neglecting all off-diagonal matrix elements of the self-energy and linearizing the eigenvalue equation around the DFT eigenvalues $E_{n{\bf {k}}}^{(0)}$. The corresponding equation reads

$E_{n{\bf {q}}}=E_{n{\bf {q}}}^{(0)}+Z_{n{\bf {q}}}\langle \phi _{n{\bf {q}}}^{(0)}|\Sigma (E_{n{\bf {q}}}^{(0)})-V_{xc}|\phi _{n{\bf {q}}}^{(0)}\rangle ,\quad Z_{n{\bf {q}}}={\frac {1}{1-{\rm {Re}}\langle \phi _{n{\bf {q}}}^{(0)}|\Sigma '(E_{n{\bf {q}}}^{(0)})|\phi _{n{\bf {q}}}^{(0)}\rangle }}$

In VASP, G₀W₀ calculations can be selected using an INCAR file such as

System  = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
NELM = 1 ; ALGO = EVGW0 ! use "GW0" for VASP.5.X
NOMEGA = 50

Convergence with respect to the number of empty bands NBANDS and with respect to the number of frequencies NOMEGA must be checked carefully.

As of VASP.6 the corresponding space-time algorithm can be selected with the following setup

System  = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
NELM = 1 ; ALGO = EVGW0R ! space-time GW algorithm, available as of VASP.6
NOMEGA = 16 ! less frequency points necessary on imaginary axes

and is recommended for larger unit cells.

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 ENCUTGW.

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

 band No.  KS-energies  QP-energies   sigma(KS)   V_xc(KS)     V^pw_x(r,r')   Z         occupation Imag(sigma) 
     
     1      -9.8139     -10.2490     -11.6986     -11.0666     -17.7440       0.6884       2.0000       1.1996
     2       2.9032       2.3351     -12.1377     -11.3852     -13.3214       0.7550       2.0000       0.0987
     3       2.9032       2.3351     -12.1377     -11.3852     -13.3214       0.7550       2.0000       0.0987
     4       2.9032       2.3351     -12.1377     -11.3852     -13.3214       0.7550       2.0000       0.0987
     5       3.4884       4.3792      -9.8672     -11.0531      -5.9700       0.7511       0.0000      -0.1035
     6       6.5738       6.9295      -8.9534      -9.4305      -4.7050       0.7455       0.0000      -0.1463
     7       6.5738       6.9295      -8.9534      -9.4305      -4.7050       0.7455       0.0000      -0.1463
     8       6.5738       6.9295      -8.9534      -9.4305      -4.7050       0.7455       0.0000      -0.1463

The first column is the band index, the second column are the quasi-particle energies $E_{n{\bf {q}}}$.

Partially self-consistent EVGW₀ calculations

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

$G({\bf {r}},{\bf {r}}',\omega )=\sum _{n{\bf {k}}}{\frac {\phi _{n{\bf {k}}}^{*(0)}({\bf {r}})\phi _{n{\bf {k}}}^{(0)}({\bf {r}}')}{\omega -\epsilon _{n{\bf {k}}}}}$

but keeping $W$ fixed to the initial DFT level W₀.

This can be achieved in two ways. If the spectral method is not selected (LSPECTRAL=.FALSE., requiring much more compute time), the quasi-particle (QP) shifts are iterated automatically four times, and you will find four sets of QP shifts in the OUTCAR file. The first one corresponds to the G₀W₀ case, the final one to the GW₀ results. The INCAR file is simply:

System = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = EVGW0 ! use "GW0" in VASP.5.X
LSPECTRAL =.FALSE.

For technical reasons, it is not possible to iterate G 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 NELM-tag. Usually three to four iterations are sufficient to obtain accurate QP shifts.

System = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = EVGW0 ! use "GW0" in VASP.5.X
NELM = 4

VASP.6 supports the space-time implementation of this procedure and can be selected via

System = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = EVGW0R ! space-time algorithm, available as of VASP.6
NOMEGA = 16   ! less points necessary on imaginary axes
NELM = 4 

Again, we recommend to employ the the 'R' option for larger unit cells.

Self-consistent QPGW₀ calculations

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

System = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = QPGW0 ! or "scGW0" for VASP.5.2.11 and older 
NELM = 4

In this case, the orbitals are updated as well by constructing a hermitian (energy independent) approximation to the self-energy. The "static" COHSEX approximation can be selected by setting NOMEGA = 1. 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 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. In general, this method is more reliable for metals and materials with strong charge sloshing.

As of VASP.6 the corresponding space-time implementation can be selected with an INCAR file such as

System = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = QPGW0R ! space-time algorithm, available as of VASP.6 
NOMEGA = 16   ! less points necessary on imaginary axes 
NELM = 4

After every iteration VASP writes following lines into the OUTCAR

 band No. DFT-energies  QP-energies  QP-e(diag)   sigma(DFT)    Z            occupation
 
     1      -9.8139     -10.2500     -10.2490     -10.4460       0.6884       2.0000
     2       2.9032       2.3350       2.3352       2.1508       0.7550       2.0000
     3       2.9032       2.3350       2.3350       2.1506       0.7550       2.0000
     4       2.9032       2.3352       2.3350       2.1506       0.7550       2.0000
     5       3.4884       4.3792       4.3792       4.6743       0.7511       0.0000
     6       6.5738       6.9297       6.9295       7.0510       0.7455       0.0000
     7       6.5738       6.9297       6.9295       7.0509       0.7455       0.0000
     8       6.5738       6.9297       6.9295       7.0509       0.7455       0.0000

Here the forth column should be identical to the third column of the EVGW0 results discussed above, whereas 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, 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.

Self-consistent QPGW calculations

Self-consistent QPGW calculations are only supported in a QP picture. As for QPGW₀, 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 = Si
NBANDS = 96
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

Analogously to the other GW algorithms, the low scaling GW implementation can be selected by using EVGWR or QPGWR respectively.

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. Alternatively (and preferably), the user can specify an electronic iteration counter using NELM:

System = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
NELM = 3
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. As for ALGO = QPGW0 (or scGW0 in vasp.5.2.11 and older), the "static" COHSEX approximation can be selected by setting NOMEGA=1.

To improve convergence to the groundstate, 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. In general, this method is more reliable for metals and materials with strong charge sloshing.

Caveats

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

System = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = QPGW  ! or "scGW" in VASP.5.2.11 and older, eigenvalues and one-electron orbitals
NELM = number of steps

However, one caveat applies to this case: when the orbitals are updated, the derivatives of the orbitals with respect to k (stored in the 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 WAVEDER file alltogether (LOPTICS=.TRUE. should not be used in the preparatory DFT runs).

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 $k$:

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

The combination LOPTICS=.TRUE.; LPEAD=.TRUE. is required since δ(H-ε_nkS)/δk_i is not available for GW like methods. 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).

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.

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

Self-consistent GW

The GW implementations in VASP described in the papers of Shishkin et al. avoids storage of the Green's function $G$ 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 $\chi$ and $\Sigma$ in reciprocal space.

As of VASP.6 a new cubic scaling GW algorithm (called space-time implementation in the following) can be selected. This approach follows the idea of Rojas et al. and allows to solve the Dyson equation for the Green's function self-consistently providing access for true self-consistent GW calculations. The main procedure for a self-consistent GW calculation consists of four steps

• Determine the Green's function $G={\frac {G_{0}}{1-\Sigma G_{0}}}$
• Evaluate the polarizability $\chi =GG$
• Calculate the screened potential $W={\frac {V}{1-\chi V}}$
• Compute the self-energy $\Sigma =GW$
• Repeat until self-consistency reached
• Solve the eigenvalue problem

This algorithm can be selected via

ALGO = scGWR ! self-consistent GW 
NELM = 4     ! number of self-consistentcy steps 

and updates the Green's function as well as the screened potential in every self-consistency step. Replacing scGWR by scGW0R updates the Green's function only, similar to the QPGW and EVGW options. After NELM iterations the self-energy is analytically continued on the real frequency axis using a continued fraction and the same linearized eigenvalue equation of the QPGW routines is solved. The resulting quasi-particle energies can be found in the OUTCAR file as discussed above.

Furthermore, the space-time implementation can be used to perform QPGW and EVGW calcualtions as well by adding a R to the ALGO GW options described above.

Large systems

Note, that the cubic scaling GW algorithm works on the imaginary time and frequency axis and makes extensive use of Fast Fourier Transformations (FFT). The frequent usage of FFTs gives a relatively large overall scaling pre-factor, so that for small unit cells the QPGW and EVGW ALGO options are often more efficient than the QPGWR and EVWGR options. Only for larger unit cells (more than ~30 atoms) we recommend to use the space-time implementation.

Here similar settings for PRECFOCK, NOMEGAPAR, NTAUPAR and NOMEGA should be chosen as for ACFDTR/RPAR calculations, with one difference. OMEGAMIN(OMEGAMAX) should be smaller(larger) than the HOMO(LUMO) w.r.t. the Fermi energy.

The reason for this is the additional Fourier transform (FT) has to capture all possible one-particle energies in the Green's function

$G(i\tau )\propto e^{-(\epsilon _{n{\bf {k}}}-\mu )\tau }(\Theta (\tau )(1-n_{n{\bf {k}}})-\Theta (-\tau )n_{n{\bf {k}}})$

Only then, the FT $G(i\tau )W(i\tau )=\Sigma (i\tau )\to \Sigma (i\omega )$ is carried out with small error.

Recipes

• G0W0 calculations: single-shot EVGW0 (G₀W₀) calculations
• EVGW0, QPGW0 calculations: partially self-consistent, with respect to G (ALGO=EVGW0 or ALGO=QPGW0)
• EVGW and QPGW calculations: self-consistent QPGW (ALGO=EVGW or ALGO=QPGW)
• DM calculations with GW routines: Determination of the frequency dependent dielectric matrix (DM) using the GW routines

Related Tags and Sections

• ALGO for response functions and GW calculations
• LMAXFOCKAE
• NOMEGA, NOMEGAR number of frequency points
• LSPECTRAL: use the spectral method for the polarizability
• LSPECTRALGW: use the spectral method for the self-energy
• OMEGAMAX, OMEGATL and CSHIFT
• ENCUTGW: energy cutoff for response function
• ENCUTGWSOFT: soft cutoff for Coulomb kernel
• ODDONLYGW and EVENONLYGW: reducing the k-grid for the response functions
• LSELFENERGY: the frequency dependent self energy
• LWAVE: selfconsistent GW
• NOMEGAPAR: frequency grid parallelization
• NTAUPAR: time grid parallelization

Examples that use this tag

References

---

Contents