关注：

1） GW0、GW、G0W0方法的区别及各自物理含义

2）各方法计算态密度和能带结构的步骤与异同

手册解读：

一、ALGO for response functions and GW calculations

http://cms.mpi.univie.ac.at/vasp/vasp/ALGO_response_functions_GW_calculations.html#7667

ALGO= CHI | GW0 | GW | scGW |  scGW0 | QPGW | QPGW0

(1) Default: none.

(2) ALGO=CHI 

calculates the frequency dependent response functions in the independent particle approximation and in the  RPA (or DFT). The computational effort is fairly small for bulk systems, but possibly huge for large supercells  with significant vacuum.

Usually one wants to determine the response function at the point only (optically allowed transitions with momentum transfer ), and this can be achieved by setting NKRED to the number of k-points in all three direction.

  For a k-point grid, NKRED=4 will restrict the calculation of the frequency dependent response function to  the important transitions with momentum transfer  as measured by optical experiments.

For a grid, one needs to set  NKREDX=4, NKREDY=6 and NKREDZ=8 (compare Sec. 6.73.9).

Furthermore, we note that ENCUTGW can be usually set to fairly small values (see Sec. 6.73.7).

Often one third or one quarter  of ENCUT (or less) will suffice.

Note that the  independent particle response function is independent of ENCUTGW, and the RPA converges reasonably fast with ENCUTGW. Reducing ENCUTGW, reduces the storage  requirements and compute time significantly compared to the defaults.

(3) For ALGO=GW and ALGO=GW0

 the orbitals of the previous groundstate calculations are maintained, and single shot GW calculations are performed.

If NELM is set as well, several iterations are performed,  and the eigenvalues are updated in the calculation of G (ALGO=GW0)  or W and G (ALGO=GW).

A full update of the orbitals can be performed by specifying  ALGO=SCGW and ALGO=SCGW0 (QPGW  and QPGW0  are synonymous to these setting, and available in VASP.5.2.13).

In the former case, the orbitals and eigenvalues are updated in G and W, whereas in the latter case the orbitals and eigenvalues are only updated in G.  

Convergence of the eingevalues with respect to ENCUTGW (see Sec. 6.73.7 and the number of unoccupied bands NBANDS is usually fairly slow and should be checked carefully. 【必须检验NBANDS参数】

We strongly recommend to read the following literature before performing calculations using VASP [111,112,113,114].

二、Recipe for GW calculations

http://cms.mpi.univie.ac.at/vasp/vasp/Recipe_GW_calculations.html

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.

【ISMEAR = 0,选0还是选-5？】

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:

（1） 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

（2）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 the k-points ; more precisely the matrix [compare (6.80)]

Calculation of this matrix requires the knowledge of the Hamiltonian, and therefore needs to be done in the preparatory DFT or hybrid functional run. 【所以需要之前的第一步DFT计算？】 

  （3） The actual GW calculations are performed in a second step using an INCAR file such as

(it is convenient to add  a single line):

System  = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
LOPTICS = .TRUE.
ALGO = GW0 ; NOMEGA = 50  

#与LSPECTRAL = .TRUE. 的区别 

The head and wings of the dielectric matrix are constructed using k.p perturbation theory (this requires the file WAVEDER).

In the present release the interaction between the core and the valence electrons is always treated on the Hartree Fock level [111].

For hybrid functionals, the three step procedure will accordingly involve the following INCAR files.

In the first two steps,

（1）converged HSE03 orbitals are determined (usually HSE03 calculations should be preceeded by standard DFT calculations,

we have not documented this step here, see Sec. 6.71.11):  

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

#参数释疑  HFSCREEN = 0.3 or 0.6?

HFSCREEN determines the range separation parameter in range separated hybrid functionals. In combination with PBE potentials, attributing a value to HFSCREEN will switch from the PBE0 functional (in case LHFCALC=.TRUE.) to the closely related HSE03 or HSE06 functional [93,94,95].

To conform with the HSE06 functional you need to select (HFSCREEN=0.2) [93,94,95].

参考：http://cms.mpi.univie.ac.at/vasp/vasp/HFSCREEN_LTHOMAS.html#7072

（2）Second determine the HSE03 orbitals for unoccupied states:

System  = Si unoccupied orbitals
NBANDS = 96
ALGO   = Exact    ! perform exact diagonalization
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.

 #还是需要设置LOPTICS参数？

As before, in the GW step, the head and the wings of the response matrix are  determined by reading the required data from the WAVEDER file.  

System  = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = GW0 ; NOMEGA = 50

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

【NBANDS和NOMEGA的收敛测试很重要】

三、Recipe for partially selfconsistent GW calculations 

http://cms.mpi.univie.ac.at/vasp/vasp/Recipe_partially_selfconsistent_GW_calculations.html

In most cases, the ``best'' results (i.e. closest to experiment)  are obtained by iterating only G, but keeping W fixed to the initial DFT W.

This can be achieved in two ways.

 (1) LSPECTRAL=.FALSE.

If the spectral method is not selected (LSPECTRAL=.FALSE. requiring much more compute time),  the 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  GW case, the final one to the GW results. The INCAR file is simply:

System  = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = GW0 ; NOMEGA = 50 ; LSPECTRAL=.FALSE.

（2）LSPECTRAL=.TRUE.

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 = GW0 ; NOMEGA = 50
NELM = 4

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

System  = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = scGW0 ; NOMEGA = 50  | or ALGO = QPGW0
NELM = 4

In this case, the orbitals are updated as well by constructing a Hermitian (energy independent) approximation to the self-energy [114].

The ``static'' COHSEX approximation can be selected by setting NOMEGA = 1 [115]. To improve convergence to the groundstate, the charge density (and the charge density only)  is mixed using a Kerker type mixing in VASP.5.2.13 and more recent versions (see Sec. 6.49). 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 is set the ``Damped''.  In general, this method is more reliable for metals and materials with  strong charge sloshing.

四、 Recipe for selfconsistent GW calculations

http://cms.mpi.univie.ac.at/vasp/vasp/Recipe_selfconsistent_GW_calculations.html

Selfconsistent GW calculations are only supported in a QP picture. 

As for GW0, it is possible to update the eigenvalues only  (ALGO=GW),  or the eigenvalues and one-electron orbitals (ALGO=scGW).

In all cases, a quasiparticle picture is maintained, i.e. satellite peaks  (shake ups and shake downs) can not be accounted for in the selfconsistency cycle. Selfconsistent GW calculations can be performed by simply repeatedly calling VASP using:  

System  = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = GW      # eigenvalues only  or alternatively
ALGO = scGW    # eigenvalues and one electron orbitals

For scGW0 or scGW non diagonal terms in the Hamiltonian are accounted for, e.g. the linearized QP equation is diagonalized, and the one electron orbitals are updated.[114]

Alternatively (and preferably),   the user can specify an electronic iteration counter using  NELM:  

System  = Si
NBANDS = 96
ISMEAR = 0 ; SIGMA = 0.05
ALGO = GW  # or  ALGO = scGW
NELM = 3

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 = scGW (or ALGO = QPGW in VASP.5.2.13),  the one electron energies and one electron orbitals are updated 3 times.

[114]

As for ALGO = scGW0, the ``static'' COHSEX approximation can be selected by setting  NOMEGA = 1 [115]. To improve convergence to the groundstate, the charge density is mixed using a Kerker type mixing starting with VASP.5.2.13 (see Sec. 6.49). 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 is set the ``Damped''.  In general, this method is more reliable for metals and materials with  strong charge sloshing.

五、参数释疑：NBANDSGW， ECUTGW

1. NBANDSGW Number of orbitals updated in GW

NBANDSGW  = [integer] twice the number of occupied states  

The flag NBANDSGW determines how many QP energies are calculated and updated in GW type calculations. This value usually needs to be increased somewhat for partially or fully selfconsistent calculations. Very accurate results are only obtained when NBANDSGW approaches NBANDS, although this  dramatically increases the computational requirements.

2. ENCUTGW= [real] (energy cutoff for response function)

Default: ENCUTGW=ENCUT

The parameter ENCUTGW controls the basis set for the response functions in exactly the same manner as ENCUT does for the orbitals. In the GW case,  updates of the response function dominate the computational work load:

(6.81)

The ENCUTGW  controls how many vectors are included in the  the response function .

Tests have shown that choosing ENCUTGW=ENCUT yields essentially exact results.

In principle, however, the response function contains contributions up to twice the plane wave cutoff (see Sec. 7.2). Since the diagonal of the dielectric matrix converges rapidly to one, such a large cutoff is never actually required (the present release has only been tested for  ENCUTGWENCUT, and might crash if  ENCUTGWENCUT).

Furthermore, in most cases, it is even possible to set ENCUTGW to a value between 150 to 200 eV, and even 100 eV gives usually QP shifts that are accurate to within a few hundreds of an eV (0.01-0.02 eV). This can help to speed up the calculations significantly and reduces the memory requirements substantially.

The flag PRECFOCK (Sec. 6.71.5), determines the FFT grid in all  GW (and Hartree-Fock) related routines. For small systems (which are often dominated by  FFT operations),  it can have a significant impact on the compute time for  QP calculations. For large systems, the FFT's usually do not dominating the computational work load and savings are expected to be small for PRECFOCK = fast .  QP shifts are usually not very sensitive to the setting of PRECFOCK (and it therefore does not harm to set  PRECFOCK = fast ), whereas for RPA calculations we recommend to set PRECFOCK= normal to avoid numerical errors.

 3.  LSPECTRAL: use the spectral method

LSPECTRAL= .FALSE. | .TRUE.

Default: LSPECTRAL=.TRUE. if NOMEGA2.

If LSPECTRAL=.TRUE. is set, the imaginary part of the independent particle  polarizability is calculated first, and afterwards the full independent particle polarizability  is determined using a Kramers-Kronig (or Hilbert) transform [111]. This reduces the computa

tional work load by almost a factor . The downside of the coin is that the response function must be kept in memory for all considered frequencies, which can cause excessive memory requirements.

VASP therefore distributes the dielectric functions among the available compute nodes.  

A similar trick is used when the QP-shifts are calculated. In general it is strongly recommended to set LSPECTRAL=.TRUE., except if memory requirements are too excessive

4. LOPTICS: frequency dependent dielectric matrix

LOPTICS= .TRUE. | .FALSE.

Default: LOPTICS=.FALSE.

If  LOPTICS=.TRUE., VASP calculates the frequency dependent dielectric matrix after the electronic ground state has been determined. The imaginary part is determined by a summation over empty states using the equation:

Note that local field effects, i.e. changes of the cell periodic part of the potential are neglected in this approximation.

These can be evaluated using either the implemented density functional perturbation theory (see Sec. 6.72.4) or the GW routines (see Sec. 6.73).

 Furthermore the method selected using  LOPTICS=.TRUE. requires an appreciable number of empty conduction band states. Reasonable results are usually only obtained, if the parameter NBANDS is roughly doubled or tripled in  the INCAR file with respect to the VASP default.

Furthermore it is emphasized that the routine works properly even for Hartree-Fock and screened exchange  type calculations and hybrid functionals.

In this case, finite differences are used to determine the derivatives of the Hamiltonian with respect to .

Note that the number of frequency grid points is determined by the parameter NEDOS (see Sec. 6.37). In many cases it is desirable to increase this parameter significantly from its default value. Values around 2000 are strongly recommended.

5.NELM, NELMIN and NELMDL-tag

NELM= [integer] NELMIN= [integer] NELMDL= [integer]  

NELM gives the maximum number of electronic SC (selfconsistency)  steps which may be performed. Normally, there is no need to change the default value: if the self-consistency loop does not converge within 40 steps, it will probably not converge at all.

In this case you should reconsider the tags IALGO, (ALGO), LDIAG , and the mixing-parameters.