G0W0 in buk water

[ioachim_dusa]

Hello everyone,

I am having problems trying to reproduce the literature bulk water band gap at G0W0 level starting from PBE. The input structure consists of 32 water molecules and is thermalised at ambient conditions. The computations are all done at gamma point and I am using the O_GW and H_GW pseudo-potentials. On this structure I am able to reproduce the band gap at PBE level of theory (about 4.2 eV), however no matter what I try, the G0W0 calculation produces band gaps of about 13 eV, which are way off compared to the literature value of 8.0 eV.

My whole setup is as suggested in VASP tutorials.

1. Start a PBE calculation with a the exact algorithm with the INCAR file
ENCUT = 800
ISMEAR = 0; SIGMA = 0.01; NBANDS = 2048
ALGO = Exact
GGA = PE
KPAR = 1
NPAR = 128
LOPTICS = .TRUE.
LPLANE=.FALSE.
EDIFF = 1E-8

2. Perform G0W0
NOMEGA = 200; LSPECTRAL=.TRUE.
ENCUT = 800
ENCUTGW = 530
ISMEAR = 0; SIGMA = 0.01; NBANDS = 2048
ALGO = EVGW0
KPAR = 1
NPAR = 128
LPLANE=.FALSE.
EDIFF = 1E-8

I have varied the number of bands up to 4000 (the system has 130 occupied bands). Moreover, I have varied the number of omega points up to 400 and the ENCUT value up to 800. I played with the ENCUTGW as well, however it seems like I am not approaching any convergence to the desired 8.0 eV value.

Any help would be appreciated! Thank you!

[michael_wolloch]

Dear Ioachim Dusa,

Please attach zipped folder containing your input files (INCAR, KPOINTS, POSCAR, POTCAR) and relevant output files (at least OUTCAR) as per the forum guidelines.
This makes analyzing your issue and giving you advice much easier.

Thanks a lot,
Michael Wolloch

[ioachim_dusa]

Dear Michael Wolloch,

thank you for your response. Here are the relevant input files. Any advice is highly appreciated!

Best wishes,
Ioachim

[alexey.tal]

Dear Ioachim,

I have some experience with calculating GW band gap in water, so let's see if I can help you with your question.

Let me point out first that there has been a number of papers that reported band gap of water calculated via G0W0 and these numbers have to be compared with a lot of caution as the results strongly depend on the chosen approximation and details of the calculation.

• The first thing I noticed is the low density in your system. Your water is around 0.86 g/cm3, which is too low compared to the other papers, which used a value much closer to experiment ~ 1 g/cm3. I think that your value is found through a calculation without vdW corrections. See what density was used in SI to these papers: Chen et al., Phys. Rev. Lett., 186401 (2016) and Gaiduk et al., Nat. Commun., 247 (2018).
• Furthermore, according to Kharche et al., Phys. Rev. Lett., 176802 (2014) (see Erratum), hard potentials (O_h_GW and H_h_GW) open the gap by 0.3 eV compared to the standard GW potentials and thus should be used in G0W0 calculations of water.
• According to Ziaei et al., J. Chem. Phys., 064508 (2016) the use of plasmon-pole approximation (PPA) results in reduction of the gap by 0.5 eV. In VASP we do not use PPA which should be considered when comparing the band gaps from literature.
• In principle, all plane waves should be included in the response function calculation. In your OUTCAR the number of plane waves is maximum number of plane-waves: 28682, hence NBANDS should be set to this values (see wiki) Of course, this makes the calculation very demanding. So one can instead extrapolated the gap based on calculations with reduced NBANDS. This should be considered when results are compared to other papers.
• VASP provides a low-scaling GW algorithm that is much faster for large systems, such as the one you are studying, so you can try to use it.
• The dielectric constant in your OURCAR is around 4.4671 which is much higher than the value found via RPA at PBE or experiment, i.e., ~ 1.8 (Garbuio et al., Phys. Rev. Lett., 137402 (2006)).
• You should double check if such a large number of frequencies is required. From my experience with NOMEGA=60-100 the gap should be well converged within 50 meV.
• I would point out that the PBE gap should be around 4.35, so your value appears to be too small. See Chen et al., Phys. Rev. Lett., 186401 (2016) and Gaiduk et al., Nat. Commun., 247 (2018).

I hope this helps.

Best wishes,
Alexey

[ioachim_dusa]

Dear Alexey,

thank you very much for your answer!
I have followed some of your points and yet I still do not obtain the right band gap value.
- I have picked a structure to have a density closer to experiment (In particular I have chosen a structure created with the PBE-rVV10 density functional from the data set [https://archive.materialscloud.org/records/vactg-v2r97], used in this GW study [Chen et al., Phys. Rev. Lett. 117, 186401 (2016)]). Choosing one of these structures gives me a higher PBE gap than before.
- I have used the O_h_GW and H_h_GW hard potentials.
- I have increase the number of NBANDS (as far as I could) and I get a change of less than 0.1eV when increasing up to 10,240 bands (see attached plot).

However, I am obtaining now an even higher band gap of about 14eV, which seems to be approaching convergence w.r.t. ENCUT. (see attached plot).
I am attaching also the input and output files for the computation at the highest ECUT and NBANDS (1400eV and 10,240 bands).

Thank you again for your help! Any further advice is appreciated.
Best,
Ioachim

[alexey.tal]

- I have picked a structure to have a density closer to experiment (In particular I have chosen a structure created with the PBE-rVV10 density functional from the data set [o], used in this GW study [Chen et al., Phys. Rev. Lett. 117, 186401 (2016)]). Choosing one of these structures gives me a higher PBE gap than before.

That is a good idea. Did you take a classical trajectory or with NQE? Keep in mind that the NQEs gaps should be smaller.

- I have used the O_h_GW and H_h_GW hard potentials.

Very good.

- I have increase the number of NBANDS (as far as I could) and I get a change of less than 0.1eV when increasing up to 10,240 bands (see attached plot).

The QP should be fairly well converged with that many bands, so clearly your issue is not related to that.

However, what I find it strange is that you use parallelization over plane waves, i.e., NPAR<number of cores, and the calculation does not stop immediately with an error.
The problem is that our GW code doesn't support parallelization over plane waves, so I think that it might be a bug in the code.

Could you try to run this calculation without setting LPLANE and NPAR?

[ioachim_dusa]

Dear Alexey,

Thank you for your help so far. I have rerun the computations without LPLANE and NPAR and I am now getting a band gap of around 8.42 eV, which seems more reasonable. Therefore, it seems like the computation involving the NPAR flag combined with GW might contain a bug.
I have not been able to perform these computations with the same parameters as before (since I get memory issues with GW). However, I think I have been able to approach convergence at around ENCUT = 1100 eV (however, with a smaller ENGUTGW = 2.5 / 9 * ENCUT = 305 eV compared to the default of 1/2 * ENCUT, otherwise I again get memory issues) and a number of bands of 4096. See attachment for the convergence plots and the input and output files with the parameters from above.

Best,
Ioachim

[alexey.tal]

It is unfortunate that we didn't catch this bug sooner. Apologies for the inconvenience. We have merged the bug fix for this and it should be release with the upcoming version.

I think your band gap results look reasonably converged. You could also use the extrapolation to the infinite basis set described in Grüneis et al., Phys. Rev. Lett., 096401 (2014) and Grumet et al., Phys. Rev. B, 155143 (2018).

Another important point here is that you are calculating your G0W0 band gap for a single snapshot, so it might be hard to compare to the literature. For example, in SM here: Gaiduk et al., Nat. Commun., 247 (2018), a gap of 8.2 eV is reported for classical trajectories. However, in this reference the position of VBM was determined via linear extrapolation.

In here Ziaei et al., Phys. Rev. B, 245109 (2017) a gap of 8.5 eV was determined based on 5 snapshots.

[ioachim_dusa]

Thanks for your help so far!

I have another question regarding the pre computation of the wave function with the EXACT algorithm. As you have seen, I have computed the G0W0 quasi particle energies starting from the PBE level of theory. For the 32 water system these computations already take a substantial amount of time (in the worst case 5 hours, having 4096 bands and a cutoff of 1100eV). However, I would also like to compute quasi particle energies starting from the hybrid level like PBE0. Right now, I do not seem to be able to bring the computation to converge in a reasonable amount of time, since the GGA was already heavy enough. My setup on the cluster so far was to have 4 nodes while per node I'd have 128 cpus and 1 cpu per task. Once again the INCAR file reads

ENCUT = 1100
ISMEAR = 0; SIGMA = 0.05; NBANDS = 4096
ALGO = Exact
GGA = PE
LHFCALC = T
KPAR = 1
NPAR = 32
LOPTICS = .TRUE.
LPLANE=.FALSE.
EDIFF = 1E-8

Do you have any tips how I could speed up the computation, maybe I am missing something once again. Thank you!

Best wishes,
Ioachim

[alexey.tal]

Indeed, for large systems the exact diagonalization (ALGO=Exact) can be too computationally demanding. In such cases, one can switch to ALGO=Normal or ALGO=Damped (for hybrids) but use a more stringent convergence criterion and/or more iterations. The problem here is that in iterative diagonalization the convergence energy is evaluated based on the occupied bands, but the unoccupied states might still be not properly converged. By using a more stringent convergence criterion we essentially increase the number of iterations and improve the empty states.