All requests for technical support from the VASP group must be addressed to: vasp.materialphysik@univie.ac.at

# GW calculations

## Contents

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]}

# The GW approximation of Hedin's equations

The GW method can be understood in terms of the following eigenvalue equation^{[5]}

Here is the kinetic energy,
the external potential of the nuclei,
the Hartree potential and
the quasiparticle energies with orbitals
. In contrast to DFT, the exchange-correlation potential is replaced by the many-body self-energy
and should be obtained together with the Green's function
, the irreducible polarizability
, the screened Coulomb interaction
and the irreducible vertex function
in a self-consistent procedure. For completeness, these equations are
^{[6]}

Here the common notation was adopted and 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 instead.

This means that the equations for the polarizability and self-energy reduce to

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

Here, in the first iteration is the non-interacting Green's function

with being a set of one-electron orbitals and 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.

In the following we discuss those supported by VASP and how the calculations are performed in practice.

# 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 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 = Si NBANDS = 96 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 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: GW step

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 quasi-particle energies: G_{0}W_{0}

This is the simplest *GW* calculation and computationally the most efficient one.
A single-shot calculation is often referred to as G_{0}W_{0} 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 . The corresponding equation becomes

with the renormalization factor

In VASP, G_{0}W_{0} calculations are selected using an INCAR file such as

System = Si NBANDS = 64 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.

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_{0}W_{0} 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

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 -6.4074 -6.9316 -11.1153 -10.4469 -17.4566 0.7843 2.0000 1.3376 2 5.5804 5.1918 -11.9505 -11.4385 -13.4545 0.7589 2.0000 0.0124 3 5.5804 5.1918 -11.9505 -11.4385 -13.4545 0.7589 2.0000 0.0124 4 5.5804 5.1918 -11.9505 -11.4385 -13.4545 0.7589 2.0000 0.0124 5 8.1035 8.3477 -9.7546 -10.0823 -5.3852 0.7449 0.0000 -0.0459 6 8.1035 8.3477 -9.7546 -10.0823 -5.3852 0.7449 0.0000 -0.0459 7 8.1035 8.3477 -9.7546 -10.0823 -5.3852 0.7449 0.0000 -0.0459 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 . 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: EVGW_{0}

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 quasi-particle (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 G_{0}W_{0} case. The INCAR file is simply:

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

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 NELM-tag. Usually three to four iterations are sufficient to obtain accurate QP shifts.

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

### Partially self-consistent quasi-particle calculations: QPGW_{0}

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 = Si NBANDS = 64 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 NELM = 4

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.^{[7]} 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.

charge density residualin 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.

After every iteration VASP writes following lines into the OUTCAR file

k-point 1 : 0.0000 0.0000 0.0000 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 2 5.5804 5.1885 5.1918 5.0683 0.7589 2.0000 3 5.5804 5.1885 5.1918 5.0683 0.7589 2.0000 4 5.5804 5.1885 5.1918 5.0683 0.7589 2.0000 5 8.1035 8.3022 8.3477 8.4313 0.7449 0.0000 6 8.1035 8.3022 8.3477 8.4313 0.7449 0.0000 7 8.1035 8.3022 8.3477 8.4313 0.7449 0.0000 8 8.9494 9.3052 9.3129 9.4313 0.7543 0.0000

For the first iteration, here the forth column should be identical to the third column of the G_{0}W_{0} 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, 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 EVGW and QPGW calculations

Self-consistent *QPGW* calculations are only supported in a QP picture. As for *QPGW*_{0}, 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 = 64 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 non diagonal 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 NELM:

System = Si NBANDS = 64 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.^{[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.^{[7]}

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 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 = Si NBANDS = 64 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 (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 :

System = Si NBANDS = 64 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 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.

## 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^{[8]} (called space-time implementation in the following) can be selected. This approach follows the idea of Rojas *et al.*^{[9]} and performs the GW self-consistency cycle on imaginary time and imaginary frequency axes
.

Following the 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

Here is the step function and
the occupation number of the state
. Because the Green's function is non-oscillatory on the imaginary time axis it can be represented on a coarse grid
, where the number of time points can be selected in VASP via the
NOMEGA tag. Usually 12 to 16 points are sufficient for insulators and small band gap systems, while metals require finer girds with roughly 24 points.^{[10]}

Subsequently, the irreducible polarizability is calculated from a contraction of two imaginary time Green's functions

Afterwards, the same compressed Fourier transformation as for the low scaling ACFDT/RPA algorithm is employed to obtain the irreducible polarizability in reciprocal space on the imaginary frequency axis .
^{[10]}
The next step is the computation of the screened potential

followed by the inverse Fourier transform and the calculation of the self-energy

### Low scaling, self-consistent quasi-particle calculations: QPGW0R

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.^{[8]} The previously discussed self-consistent quasi-particle GW approximations (QPGW[0]) can be performed within this formalism by adding an additional *R* to the ALGO option, *e.g.* single-shot GW calculations or partially self-consistent quasi-particle energies can be selected via the following lines

System = Si ISMEAR = 0 ; SIGMA = 0.05 LOPTICS = .TRUE. ; LPEAD = .TRUE. NELM = number of iterations wanted ! set 1 for G0W0 ALGO = QPGW0R NOMEGA = 12 ! small number of frequencies necessary

Note, that a preceding DFT calculation is not necessary with this setting, but required if one wants to start from a hybrid functional calculation. Due to the analytical continuation of the self-energy small deviations in the quasi-particle energies and orbitals can be expected. however, deviations of the QP energies around the Fermi energy are usually smaller than 0.05 eV.^{[8]}

After a successful run, VASP writes the following lines into the OUTCAR file

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 1 -6.4074 -6.8617 -6.7178 0.6832 -6.7168 0.6789 2.0000 -1.1944 0.0058 2 5.5804 5.0902 5.2064 0.7629 5.2073 0.7591 2.0000 0.0001 0.0018 3 5.5804 5.0894 5.2058 0.7629 5.2067 0.7591 2.0000 0.0001 0.0018 4 5.5804 5.0893 5.2057 0.7629 5.2066 0.7591 2.0000 0.0001 0.0018 5 8.1035 8.4039 8.3336 0.7660 8.3332 0.7632 0.0000 -0.0190 0.0031 6 8.1035 8.4026 8.3326 0.7660 8.3322 0.7632 0.0000 -0.0192 0.0034 7 8.1035 8.4023 8.3324 0.7660 8.3320 0.7632 0.0000 -0.0192 0.0035 8 8.9494 9.4257 9.3101 0.7574 9.3089 0.7524 0.0000 -0.0410 0.0145

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 forth column in ALGO=*QPGW0* calculations, while column six is obtained from the roots of following equation

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

### Low scaling, fully self-consistent quasi-particle calculations: QPGWR

Similarly, low scaling self-consistent quasi-particle calculations can be chosen where the screened potential is updated as well

System = Si ISMEAR = 0 ; SIGMA = 0.05 LOPTICS = .TRUE. ; LPEAD = .TRUE. NELM = number of iterations wanted ALGO = QPGWR ! update W during self-consistency cycle as well NOMEGA = 12 ! small number of frequencies necessary

Here the same caveats can be expected as for ALGO=*EVGW* and *QPGW*. As for ALGO=*EVGW* and *QPGW*, we suggest to leave out the LOPTICS and LPEAD line for metals.

### Partially self-consistent GW caluclations: scGW0R

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 = Si ISMEAR = 0 ; SIGMA = 0.05 LOPTICS = .TRUE. ; LPEAD = .TRUE. NELM = number of iterations wanted ALGO = scGW0R ! full self-consistent GW calculations, no update in W

The number of self-consistency steps can be set with the NELM tag and usually 4 cycles suffice to reach self-consistency. After each cycle VASP writes the following lines in to the OUTCAR file

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 1 -6.4074 -6.9582 -6.7744 0.6664 -6.7723 0.6584 2.0000 -0.7917 0.5691 2 5.5804 5.0717 5.1897 0.7680 5.1910 0.7629 2.0000 -0.0026 0.0062 3 5.5804 5.0710 5.1891 0.7680 5.1905 0.7629 2.0000 -0.0026 0.0062 4 5.5804 5.0709 5.1891 0.7680 5.1904 0.7629 2.0000 -0.0026 0.0062 5 8.1035 8.3674 8.3042 0.7605 8.3038 0.7572 0.0000 0.0036 0.0129 6 8.1035 8.3663 8.3033 0.7605 8.3029 0.7572 0.0000 0.0036 0.0129 7 8.1035 8.3661 8.3032 0.7605 8.3028 0.7572 0.0000 0.0036 0.0129 8 8.9494 9.3762 9.2698 0.7508 9.2689 0.7465 0.0000 -0.0237 0.0076

Here the meaning of each column is the same as for ALGO=*QPGW0R*.

### Fully self-consistent GW caluclations: scGWR

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

System = Si ISMEAR = 0 ; SIGMA = 0.05 LOPTICS = .TRUE. ; LPEAD = .TRUE. NELM = number of iterations wanted ALGO = scGWR ! full self-consistent GW calculations, no update in W

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.

## Related Tags and Sections

- ALGO for response functions and
*GW*calculations - LOPTICS, derivative of wavefunction w.r.t.
- LPEAD, derivative of wavefunction with finite differences
- LMAXFOCKAE overlap densities and multipoles
- MAXMEM, memory available to one core
- 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

## References

- ↑
^{a}^{b}M. Shishkin and G. Kresse, Phys. Rev. B 74, 035101 (2006). - ↑
^{a}^{b}M. Shishkin and G. Kresse, Phys. Rev. B 75, 235102 (2007). - ↑
^{a}^{b}^{c}^{d}^{e}M. Shishkin, M. Marsman, and G. Kresse, Phys. Rev. Lett. 99, 246403 (2007). - ↑ F. Fuchs, J. Furthmüller, F. Bechstedt, M. Shishkin, and G. Kresse, Phys. Rev. B 76, 115109 (2007).
- ↑
^{a}^{b}M. S. Hybertsen, S. G. Louie Phys. Ref. B 34, 5390 (1986) - ↑ L. Hedin, Phys. Rev. 139, A796 (1965)
- ↑
^{a}^{b}F. Bruneval, N. Vast, and L. Reining, Phys. Rev. B 74, 45102 (2006). - ↑
^{a}^{b}^{c}P. Liu, M. Kaltak, J. Klimes and G. Kresse, Phys. Rev. B 94, 165109 (2016). - ↑ H. N. Rojas, R. W. Godby, R. J. Needs, Phys. Rev. Lett. 74, 1827 (1995)
- ↑
^{a}^{b}M. Kaltak, J. Klimeš, and G. Kresse, Journal of Chemical Theory and Computation 10, 2498-2507 (2014).