need help with HSE and GW band gap calculations

[panda]

Dear Admin,

I am working on band gap calculations for one of my research projects. I want to compare the PBE exchange-correlation functional with the HSE hybrid functional and the GWo method. To ensure I am getting the correct results, I decided to perform these calculations on the Si system first, as the band gap is well known. From a recent paper by Kresse the PBE value was reported as 0.65 eV and the HSE and GWo values were around 1.1 eV. However, when I perform the calculations, I get 0.65 eV for PBE, HSE, and GWo. I don't know what I am doing wrong, because the HSE exchange portion is printed in my output from the HSE calculation, and I see the dielectric matrix in the GWo output. For you to help me with this problem, I am going to paste in my INCAR file for each of the calculations, as well as the details of the calculations. In advance I would like to say that I really appreciate your help in this matter, as I have been trying to resolve it myself for some time with no avail.

INCAR file for PBE Si run:

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
Automatic mesh
0 ! number of k-points = 0 ->automatic generation scheme
G !
16 16 16 ! subdivisions N_1, N_2 and N_3 along recipr.
0. 0. 0. ! optional shift of the mesh (s_1, s_2, s_3)

After this step completes, I run a PBE band structure calculation:

ISTART = 1
ISYM = 0
ICHARG = 11
ISMEAR = 0
SIGMA = 0.05
LORBIT = 11
NSW = 0
IBRION = -1
LORBIT=11
NEDOS=2000
EMIN=-30.0
EMAX=30.0

with KPOINTS file:
k-points along high symmetry lines
20 ! 20 intersections
Line-mode
rec
0 0 0 ! gamma
0.5 0.0 0 ! X

0.5 0.0 0 ! X
0.5 0.5 0.0 ! M

0.5 0.5 0.0 ! M
0.0 0.0 0.0 ! gamma

0.0 0.0 0.0 ! gamma
0.5 0.5 0.0 ! M

0.5 0.5 0.0 ! M
0.5 0.5 0.5 ! R

0.5 0.5 0.5 ! R
0.0 0.0 0.0 ! gamma

I also perform a HSE calculation from the original PBE calculation (not continuing from the PBE band structure, but copying the original PBE directory after the first calculation finishes and starting from there).

INCAR for HSE calculation:

ISTART = 1
LHFCALC = .TRUE. ; HFSCREEN = 0.2
NBANDS = 51
ALGO = All ; TIME = 0.4


KPOINTS
Automatic mesh
0 ! number of k-points = 0 ->automatic generation scheme
G !
16 16 16 ! subdivisions N_1, N_2 and N_3 along recipr.
0. 0. 0. ! optional shift of the mesh (s_1, s_2, s_3)

After this calculation completes, I make a new KPOINTS file. I take the kpoints along ditrections of high symmetry from vasprun.xml in the PBE band gap calculation (120 kpoints for the lines I chose) and make a file called Si.kpts. Then I move that file to the HSE directory where the first HSE calculation completed and do
cat IBZKPT Si.kpts > KPOINTS

then I go in and make the weights of all the 120 kpoints from Si.kpts be 0 and change the number at the top of the new KPOINTS file to be the sum of the kpoints from IBZKPT + Si.kpts. With the grid I chose this comes to be 265 KPOINTS.

INCAR is:
ALGO = N
NELMIN = 5
IMIX = 1
ISTART = 1
LORBIT = 11

when I analyze the results, Si is s2p2 so the HOMO is at ' 4 ' and the LUMO is at ' 5 '. I can do grep ' 4 ' EIGENVAL or grep ' 4 ' PROCAR and get all the energies. If I take the difference between the highest HOMO and lowest LUMO that gives me the band gap, which is 0.65 for PBE and HSE.

I went back and tried a different INCAR for the first step of the HSE calculation to see if that would make a difference:
# OUTPUT AND START
NWRITE=2
ISTART=1

# ELECTRONIC CONVERGENCE
NELM =400
ISMEAR=0
SIGMA=0.2
EDIFF=1E-6
PREC=Medium
ALGO=Damped
TIME=0.4

1#HSE
LHFCALC=.TRUE.
HFSCREEN=0.2
ENCUTFOCK=0

# TYPE OF RUN
IBRION=-1
NSW=0
GGA=PE

with

KPOINTS
Automatic mesh
0 ! number of k-points = 0 ->automatic generation scheme
G !
16 16 16 ! subdivisions N_1, N_2 and N_3 along recipr.
0. 0. 0. ! optional shift of the mesh (s_1, s_2, s_3)

again, making a copy of the initial PBE run and performing this HSE calculation in that directory, followed by the same procedure to calculate the HSE band gap:

INCAR is:
ALGO = N
NELMIN = 5
IMIX = 1
ISTART = 1
LORBIT = 11


Again, after this calculation completes, I make a new KPOINTS file. I take the kpoints along ditrections of high symmetry from vasprun.xml in the PBE band gap calculation (120 kpoints for the lines I chose) and make a file called Si.kpts. Then I move that file to the HSE directory where the first HSE calculation completed and do
cat IBZKPT Si.kpts > KPOINTS

then I go in and make the weights of all the 120 kpoints from Si.kpts be 0 and change the number at the top of the new KPOINTS file to be the sum of the kpoints from IBZKPT + Si.kpts. With the grid I chose this comes to be 265 KPOINTS.

I analyze EIGENVAL and PROCAR and still get 0.65 for the band gap.

Same exact thing for GWo calculations.

I make a copy of the first PBE run directory and use this INCAR for thoptical properties calculation:
NELMIN = 5
IMIX = 1
ISTART = 1
LORBIT = 11
LOPTICS = .TRUE.
ISMEAR = 0
SIGMA = 0.05

with

KPOINTS
Automatic mesh
0 ! number of k-points = 0 ->automatic generation scheme
G !
16 16 16 ! subdivisions N_1, N_2 and N_3 along recipr.
0. 0. 0. ! optional shift of the mesh (s_1, s_2, s_3)

Once the optical properties run completes, I add:

ALGO = GWo
NOMEGA = 50
ENCUTGW = 100

to the INCAR file, so now INCAR is:

ALGO = N
NELMIN = 5
IMIX = 1
ISTART = 1
LORBIT = 11
ALGO = GWo
NOMEGA = 50
ENCUTGW = 100

After the GWo calculation completes, I perform the GWo band structure calculation, with the HSE KPOINTS file with 265 kpoints, high symmetry points weighted to 0, as described above. Again, I get the band gap to be 0.65 eV. I went back and tried again with optical properties and GWo calculation in the same calculation step using the 16x16x16 gamma grid, then recomputed the band structure with the 265 KPOINTS file with high symmetry points weighted to zero with no difference in the band gap energy result. There is a slight difference in the numerical values of the energies, but it is in the fourth decimal place. The HSE and GWo band gaps should be around 1.1 eV. I have tried running these tests on at least three different computing clusters to make sure it was not a compliation, linking, or math library error. I also tried running the calculations with a different software package, and I got the correct result when I used that software (0.65 eV with PBE, 1,1eV with HSE) but I WANT TO USE VASP TO GET THIS RESULT CORRECT!!! This is seriously driving me insane, please help!


<span class='smallblacktext'>[ Edited ]</span>

[aron]

Hi Panda.

Just starting with the PBE / HSE comparison. It appears that you may have switched off HF for the final band structure calculation? i.e. there is no LHFCALC in:

"INCAR is:
ALGO = N
NELMIN = 5
IMIX = 1
ISTART = 1
LORBIT = 11"

So you ended up doing a PBE calculation at the last step by mistake? In any case, you should check the eigenvalues for your initial HSE SCF calculation.

I hope that helps - it sounds like you have spent a long time on this.

[panda]

Wow Aron, that was the problem. Oops! I went back and ran with the same INCAR only added
LHFCALC = .TRUE.
HFSCREEN = 0.2
ENCUTFOCK = 0

and got 1.17 eV for the Si band gap. I don't know why I forgot to put that in the band structure calculation, I guess I assumed the wave function was generated from the HSE calculation in the previous step but apparently I was wrong, because like you said - you wind up doing a PBE calculation in the last step if you don't keep the HSE tags in the INCAR.

Now the only problem I nave is in larger systems - I tried the same thing and they never finish (15 atom unit cells). I suppose what I will have to do in this case is just pick specific KPOINTS of interest and calculate the band gap at only those points rather than trying to recompute the entire band structure.

Now also I will go back and test GW for Si starting with my HSE calculation. Thanks, I can't believe it was something this simple.