Practical guide to GW calculations
- 1 First step: DFT calculation
- 2 Step 2: GW step
- 3 Low scaling GW algorithms
A brief summary of the theoretical background can be found here.
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. Furthermore, cubic scaling algorithms are available.
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 grounstate calculation with few unoccupied orbitals only,
System = SiC 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 = 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 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 (see here for a selection of available hybrid functionals):
System = SiC 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 = 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.
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.
Single shot quasi-particle 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 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, G0W0 calculations are selected using an INCAR file such as
System = SiC NBANDS = 512 ISMEAR = 0 ; SIGMA = 0.05 NELM = 1 ; ALGO = EVGW0 ! use "GW0" for VASP.5.X NOMEGA = 50
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.
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 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: 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 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 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.
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.
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 form previous iteration.
Partially self-consistent quasi-particle 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 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.
After every iteration VASP writes 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 forth column should be identical to the third column of the G0W0 results discussed above, while the third column reports the quasi-particle energies obtained from including the off-diagonal matrix elements in the eigenvalue equation.
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 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 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 = SiC NBANDS = 512 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.
Additional information about this method is found here.
Fully self-consistent QPGW calculations with an update of the orbitals in and  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 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 = SiC NBANDS = 512 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. 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 (called space-time implementation in the following) can be selected. This approach follows the idea of Rojas et al. and performs the GW self-consistency cycle on imaginary time and imaginary frequency axes .
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 quasi-particle energies (from an analytical continuation of the self-energy to the real axis)
- Natural orbitals from the first order change of the density matrix (i.e. ), see the NATURALO tag for more information.
N.B. This selection ignores NELM.
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 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 an other set of QP-energies that is obtained from the roots of following equation
These roots represent the poles of the Green's function in the spectral representation.
Partially self-consistent GW caluclations: GW0R or 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 = SiC ISMEAR = 0 ; SIGMA = 0.05 LOPTICS = .TRUE. ; LPEAD = .TRUE. NELM = number of iterations wanted 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 NELM 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 or scGWR
System = SiC ISMEAR = 0 ; SIGMA = 0.05 LOPTICS = .TRUE. ; LPEAD = .TRUE. NELM = number of iterations wanted ALGO = GWR ! ALGO = scGWR has the same effect here, that is self-consistency in G and W
The cubic scaling space-time RPA as well as GW algorithm require considerably more memory than the correspondong 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. The optimum value of MPI groups that share the same time points will then be set internally to an optimum value.
The required storage for an 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
NKPTS * (NGX*NGY*NGZ)^2 / ( NCPU / NTAUPAR ) * 16
where "NKPTS" is the number of irreducible q-points, "NCPU" the number of MPI ranks used for the job and "NGX,NGY,NGZ" the number of FFT grid points for the supercell, which is written in the OUTCAR file in the line
FFT grid for supercell: NGX = 32; NGY = 32; NGZ = 32
The smaller NTAUPAR is set, the less memory the job will require to finish successfully. Note that VASP finds the optimum value of NTAUPAR based on MAXMEM, the freely available memory per MPI rank on each node. Thus it is recommended not to set NTAUPAR in the INCAR, but to set MAXMEM instead and allow VASP to find the optimum NTAUPAR.
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 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 mpi rank on each node
- 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
- NATURALO, natural orbitals
- ↑ 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).
- ↑ M. S. Hybertsen, S. G. Louie Phys. Ref. B 34, 5390 (1986)
- ↑ a b F. Bruneval, N. Vast, and L. Reining, Phys. Rev. B 74, 45102 (2006).
- ↑ a b 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)
- ↑ B. Ramberger, Z. Surkuma, T. Schäfer, G. Kresse, J. Chem. Phys. 151, 214106 (2019).
- ↑ M. Grumet, P. Liu, M. Kaltak, J. Klimeš and Georg Kresse, Phys. Ref. B 98, 155143 (2018).