关注：

1） GW方法的物理含义；

2） VASP 计算过程及参数设置；

4） 手册是最好的老师，很多问题手册上都有答案：

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

1、To obtain band structure plots from the GW method, first you have to perform a DFT calculation. 

The INCAR file I used for DFT on Si looks like this:

ISTART = 0

ISYM = 0

ICHARG = 2

ISMEAR = 0

SIGMA = 0.1

LORBIT = 11

NSW = 200

IBRION = 2

ISIF = 3

EDIFF = 0.0001

EDIFFG = -0.01

KPOINTS file for the run is:

Monkhorst Pack

0

 Gamma

 4 4 4

 0 0 0  

 POSCAR file is:

  Si     5.38936000000000  

    0.5000000000000000 0.5000000000000000  0.0000000000000000

    0.0000000000000000 0.5000000000000000  0.5000000000000000

    0.5000000000000000 0.0000000000000000 0.5000000000000000

    2

    Direct  

    0.0000000000000000 0.0000000000000000 0.0000000000000000

     0.2500000000000000 0.2500000000000000 0.2500000000000000

POTCAR file is supplied with VASP 5.2, remember to use paw pseudopotentials so that in the band structure calculation the Wigner-Seitz radius for Si need not be directly specified in INCAR.

2.After this run completes you are ready to perform the HF-type (HSE06) calculation. Change

the INCAR file to:

ISTART=1

ICHARG=2

EDIFF=0.0001

EDIFFG=-0.01

ENCUT=400

ENAUG=800

LREAL=.FALSE.

LWAVE=.TURE.

NELM=200

NSE=0

IBRION=-1

LMAXMIX=4

ISMEAR=0

SIGMA=0.1

NSIM=4

IALGO=48

SIYM=0

LHFCALC=.TRUE.

HFSCREEN=0.2

ALGO=Damped

TIME=0.4

ENCUTFOCK=0NKRED=2

AEXX=0.25

keep the other files constant for this run (e.g. KPOINTS, POSCAR, POTCAR need not be changed)

3. once this run has completed, add to the bottom of the INCAR file from the HSE06 calculation:

LOPTICS=.TRUE.

This step calculats the dielectric matrix

4. Once this run has completed ,change

ALGO=Damped

to

ALGO=GWo

and

ICHARG=2

to ICHARG=11

and add:

NOMEGA=50

ENCUTGW=100

LORBIT=11

and change the KPOINTS file to:

k-points along high symmetry lines

20  !20 intersections

Line-mode

rec

0 0 0  !gamma

0.5 0.0 0.0 !X

.....

These values are the points of high symmetry along the lines of high symmetry within the first Brillouin Zone of the Si unit

cell structures in reciprocal space. And that's it !

  Once complete you can plot the PROCAR, EIGENVAL, and DOSCAR files from the output of this run as described in the previous posts.

  The band structure plot produce from the plotting utility is shown below.

附手册解读 ：Frequency dependent GW calculations

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

Frequency dependent GW calculations

Available as of VASP.5.X. For details on the implementation and use of the GW routines we recommend the following references: Ref. [111,112,113,114].

Subsections

• ALGO for response functions and GW calculations

• LMAXFOCKAE

• NOMEGA, NOMEGAR number of frequency points

• LSPECTRAL: use the spectral method

• OMEGAMIN, OMEGAMAX, OMEGATL and CSHIFT

• NBANDSGW Number of orbitals updated in GW

• ENCUTGW energy cutoff for response function

• ENCUTGWSOFT soft cutoff for Coulomb kernel

• ODDONLYGW and EVENONLYGW and NKRED: reducing the -grid for the response functions

• LSELFENERGY: the frequency dependent self energy

• LWAVE: selfconsistent GW

• Recipe for GW calculations

• Recipe for partially selfconsistent GW calculations

• Recipe for selfconsistent GW calculations

• Caveats for selfconsistent quasiparticle GW calculations

• Using the GW routines for the determination of frequency dependent dielectric matrix

• scGW0 caveats

112M. Shishkin and G. Kresse, ``Self-consistent calculations for semiconductors and insulators'', Phys. Rev. B 75, 235102 (2007).

113F. Fuchs, J. Furthmüller, F. Bechstedt, M. Shishkin, and G. Kresse, ``Quasiparticle band structure based on a generalized Kohn-Sham scheme'', Phys. Rev. B 76, 115109-1-8 (2007).

114M. Shishkin, M. Marsman, and G. Kresse, ``Accurate quasiparticle spectra from self-consistent GW with vertex corrections'', Phys. Rev. Lett. 99, 246403 (2007).

115F. Bruneval, N. Vast, and L. Reining, Phys. Rev. B 74, 45102 (2006).

（1）Recipe for selfconsistent GW calculations

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 orbitalsFor 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 = 3In 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.

附--网言摘录

Recipe for G0W0 calculations

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 orbitalsand 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 precissely the matrix [compare (6.76)] 

Calculation of this matrix requires the knowledge of the Hamiltonian, and therefore needs to be done in the preparatory DFT or hybrid functional run. 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 = 50The 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 [104].

For hybrid functionals, the three step procedure will accordingly involve the following INCAR files. In the first two steps, 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.69.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 orbitalsSecond 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.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 = 50Convergence with respect to the number of empty bands NBANDS and with respect to the number of frequencies NOMEGA must be checked carfully.

附. 参数说明

ENCUTFOCK: FFT grid in the HF related routines

ENCUTFOCK = real number (energy cutoff determining the FFT grids in the HF related routines)

default: none

The ENCUTFOCK parameter controls the FFT grid for the HF routines. The only sensible value for ENCUTFOCK is

ENCUTFOCK=0. This implies that the smallest possible FFT grid, which just encloses the cutoff sphere corresponding to the

plane wave cutoff, is used. This accelerates the calculations by roughly a factor two to three, but causes slight changes in the

total energies and a small noise in the calculated forces. The FFT grid used internally in the Hartree Fock routines is written

to the OUTCAR file. Simply search for lines starting with

FFT grid for exact exchange (Hartree Fock)

In many cases, a sensible approach is to determine the electronic and ionic groundstate using ENCUTFOCK = 0, and to

make one final total energy calculation without the flag ENCUTFOCK.

网络摘录：

http://emuch.net/html/201312/6800622.html

 NOMEGA, NOMEGAR number of frequency points

NOMEGA= [integer] (number of frequency points)

NOMEGAR= [integer] (number of frequency points along real axis)

Default:
NOMEGA	=50	for calculations

谢谢你的回帖。我依照你的建议再做一次Si的测试。
我在计算过程中注意了CPU的要求，NBAND和NOMEGA和CPU核心数目是一样的。Si的赝势，Si_GW，ENMAX只有245.345 eV，我使用的300 eV应该已经足够了，不过我会调大再看看结果。
现在看来，您提到的空带数，K网格，频率采样密度是最可能的问题原因。可惜，这三个参数都和计算量有直接的关系。现在的参数设置下，GW的计算已经非常慢了，实际考察的复杂结构所需的时间估计也会很长。我看了Abinit GW计算部分的介绍，使用Eective Energy Technique (EET) 可以避免使用上百条的空带也能得到很好的结果。如果VASP速度实在不可接受，最后也考虑试一下Abinit  

liqizuiyang(站内联系TA) 

  对于Si这样的元素，可以用Plasmon-Pole模型来算介电函数，这样就不用针对每个频率点都算一次介电函数了，可以节约不少计算量。
但是vasp只能算全频率GW，因此计算量很大。abinit支持4种PPA模型，对这四种模型的测试可以看下PRB 88, 125205(2013)这篇文章。

wekow(站内联系TA)

重新优化结构计算了Si的能隙，NBAND=240，ENCUT=300，KPOINTS 8*8*8，NOMEGA = 80, GW计算出来的能隙仍然是在1.8 eV。我搜索到有些Si GW计算的结构，即使使用高得多的收敛参数，仍然能得到接近实验值的结果。我的计算结构差距太大了，现在很怀疑是不是哪里做错了。
这是在GW计算上一步中计算介电矩阵的INCAR。计算介电矩阵这一步使用了静态计算的结果。
 PREC = Accurate
  ISTART =     1
  ICHARG =     2  
  ISPIN =      0  
  ENCUT = 400  
  ISYM = 2  
  GGA =  PS    
  NBANDS = 240
!GW0
  LOPTICS = .TRUE.
  NEDOS = 2000
  LPEAD = .TRUE.
!Electronic Relaxation
  ALGO = Exact
  EDIFF =  1E-05    stopping-criterion for ELM
!  LREAL = AUTO      real-space projection
!Ionic relaxation
  EDIFFG =  -1E-03  stopping-criterion for IOM
  ADDGRID= .TRUE.   reduce the noise in the forces
!DOS related values:
  ISMEAR =  0      gaussian smear 0 tetrahedron method -5
  SIGMA = 0.05

ALGO = Exact or  ALGO = Diag performs an exact diagonalization (IALGO = 90), and we recommend to use this if more than 30-50 % of the states are calculated (e.g. for or RPA calculations).