Query Regarding K-Point Dependence of GW0 Band Gap

[alok_shukla1]

Dear Experts,

I am calculating the band gap of a 6-atom heterostructure using the GW0 method in VASP. I am currently checking the convergence of the calculated band gap with respect to the k-point sampling.

However, I find that the GW0 band gap changes significantly when I change the k-point mesh. The band gaps I obtained are:

k-point mesh GW_0 band gap (eV)
6*6*1 2.2066
7*7*1 2.4651
9*9*1 2.1481
11*11*1 2.0854
15*15*1 2.0610
17*17*1 1.79
21*21*1 2.0006

As can be seen, the band gap does not show a clear convergence trend with increasing k-point density. The system contains only 6 atoms, but the GW0 calculations are computationally expensive. With my available computational resources, a 21*21*1 mesh is approximately the maximum k-point density that I can practically use.

I would therefore like to ask whether there may be something in my calculation setup that I should check or modify.

I have attached my POSCAR and INCAR files. I need suggestions to obtain a reliable GW0 band gap.
https://drive.google.com/drive/folders/ ... drive_link

[christopher_sheldon1]

Dear alok_shukla1,

Thank you for your question. I cannot access your Google Drive. Could you attach the files here so that I can read them? Or permit me to access the Google Drive.

At a glance, it looks like the 17x17x1 step is out of step with the series, so I'd check to see if anything has gone wrong there.

Best wishes,

Chris

[christopher_sheldon1]

Could you also upload the OUTCAR files? This will help to diagnose the problem.

Best,

Chris

[alok_shukla1]

Dear Chris,
Yes, files are uploaded. Please check.

[christopher_sheldon1]

Dear Alok,

Thank you. I try to repeat the calculations and see if I can find out what has been going on.

Best wishes,

Chris

[alok_shukla1]

Dear Chris,
Did you find anything? It would be helpful if you could please share.

[christopher_sheldon1]

Dear Alok,

17x17 being out of place is due to it not including the high-symmetry K point in you k-point mesh. This is the point where the valence band minimum is (VBM). Without K, you find the gap at a nearby point, and increase it, hence 7x7 being so large. These are the Gamma, M, K, and K' high-symmetry points in the following k-meshes:

Mesh Γ M K, K′
3×3 yes no yes
6×6 yes yes yes
7×7 yes no no
8×8 yes yes no
9×9 yes no yes
10×10 yes yes no
11×11 yes no no
12×12 yes yes yes
14×14 yes yes no
15×15 yes no yes
17×17 yes no no
18×18 yes yes yes
21×21 yes no yes

Since the VBM is at K, being a hexagonal cell, you need to converge this with a k-mesh that is a multiple of 3, i.e., 3x3, 6x6, 9x9. This gives me the following trend:

gap_vs_kmesh_3n.png

So it converges more smoothly. Mine are G0W0 tests but the same is true for GW0, as you can see in your results (green line).

I can't find any specific issue with the 17x17, but it does not include the K point, so does not look at the VBM directly. I'd need the files of the other calculations to try to understand better. Let me know if this answers your question. Sorry about the delay, the calulations took a while.

Best wishes,

Chris

[alok_shukla1]

Dear Chris,

Thank you for your response and for sharing your results. They are very helpful.

Following your suggestion, I performed single-shot G0W0 calculations using 18*18*1 and 21*21*1 k-point meshes, both of which are multiples of 3. I obtained band gaps of 1.776 eV and 1.751 eV, respectively. However, for the 18*18*1 mesh, your reported value is approximately 1.81 eV.

As per the literature survey, your result appears to be somewhat closer to the reported fundamental gap than mine.

For my calculations, I used optB86b-vdW for structural relaxation. Since optB86b-vdW is not directly supported in the standard G0W0 workflow. Therefore, after taking the CONTCAR from the relaxation step, I removed optB86b-vdW from the SCF and from the rest of the steps.

Also, should I need to use GW-specific PAW datasets (the '_GW' POTCAR files) for the G0W0 calculations? Could the choice of PAW potentials or any other computational settings be responsible for the difference in the calculated band gaps?

If possible, could you please share your input files, particularly the INCAR settings and POTCAR dataset information, so that I can compare the computational setups and better understand the discrepancy?

Thank you again for your time and help.

Best regards

[christopher_sheldon1]

Dear Alok,

I've attached my input files and OUTCARs. This was using standard POTCARs but for GW calculations, you should use GW POTCARs; sorry, I missed this earlier. You will need to repeat these convergence tests with GW POTCARs. A few things about your settings that should be updated:

• use GW POTCARs - the unoccupied states are required for GW calculations - they describe high-energy scattering states better
• if you want to use GW0, rather than G0W0, then we recommend using ALGO = EVGW0 with NELMGW = 1, rather than ALGO = GW0, which is deprecated since 6.6 (it will still work mind)
• NOMEGA = 192 is large, 50 is more normal and will reduce the calculation by around two-thirds - this should be tested for convergence - check 50 vs 100 once
• you use a rather small number of NBANDS = 96 - the default is to use all unoccupied bands - you can use less but would need to do convergence tests (this can be done for a sparse k-mesh) with respect to the band gap. I tend to put the following in the final script:

  Code: Select all

  grep "maximum number of plane-waves:"  OUTCAR.1 >nbands
  NBANDS=`awk '{ printf "%i \n", $5 }' < nbands`  ;
  echo "NBANDS = $NBANDS">>INCAR

• unnecessary tags - LVTOT, LVHAR, LORBIT, EDIFFG and the Wannier90 - these are not required for the GW bandgap, and can only add to the calculation time, so I'd remove them unless you specifically want them for other properties
• spin-orbit coupling - this may be important for your system - you should check it (for LSORBIT, you need to use vasp_ncl, rather than vasp_std. It will be more expensive, so bear that in mind)

  Code: Select all

  LNONCOLLINEAR = .TRUE.
  LSORBIT       = .TRUE.
  MAGMOM        = 6*0      ! 3 components per atom, all zero for a non-magnetic system
  

Then a few points about your structure:

• Vacuum size - the band gap with respect to vacuum should be checked - using the smallest vacuum possible while still achieving convergence will make this calculation more affordable
• symmetry - a Janus MoSeTe layer should have threefold symmetry - but VASP only finds a mirror plane (The static configuration has the point symmetry C_1h. in the OUTCAR). If you can symmetrise the system, then this reduces the number of k-points required and will make the calculation cheaper.

I hope this helps.

Best wishes,

Chris

[alok_shukla1]

Dear Chris,
Thank you for explaining things in detail. Indeed, this explanation is really helpful. I have one more query:
I am using optB86b-vdW in my calculations, and G0W0 calculations are not compatible with optB86-vdW. I used this vdW in the relaxation step along with '_GW' POTCAR files, but I removed it during the SCF run. Now, my SCF run is not converging if I am using '_GW' POTCAR files. SCF run is converging well with normal POTCAR files (I removed optB86-vdW here as well). I need some suggestions on how to perform the calculations in this case.

[christopher_sheldon1]

Dear Alok,

I would use PBE orbitals for GW. So perform an SCF with PBE and then use those for the GW values. This is the standard approach for RPA, and post-DFT methods. You can use other functionals to generate the orbitals but you may encounter difficulties. An alternative is to use a nearly converged density (CHGCAR) and start from those; that should help converge the orbitals for optB86b-vdW.

Best wishes,

Chris

[alok_shukla1]

Dear Chris,
Thank you for your suggestions.

Thanks.