HF and GW bandgap different from absorption onset

[a.ugolotti]

Dear Yamboers,

I am studying a 2D insulator, the bandgap of which I calculate from HF or GW doesn't agree with the edge of optical adsorption I calculate both for IP or LF approximations. I am not talking about the comparison with experimental results (to which the optical simulation are in good agreement) but somehow consistency beteween to approaches of the same code.
To give an example, for a not-so-badly-converged system, the calculated gap lies in the range 8-10 eV, where the GW result points towards 8 eV while the HF one is closer to 10 eV. Altough this may be expected, the optical spectrum has the onset between 2 and 3 eV, peaking around 3-4 eV, where the IP result is closer to the lower ranges and the LF result is closer to the higher values.

Additionally, I tried to get the HF band gap of the same system using a different software, namely CRYSTAL, which gives a result ~ 3 eV.

I wonder what I am missing here: I mean, other than the discrepancy with the expected result, I am puzzled by the difference between the two calculated values of the bandgap. As for the former, I have tried to converge all the paramers:

- EXX/VXC cutoffs
- with/without rim+cutoff
- HF band ranges
- scf/nscf K point mesh
- scf (Espresso) pseudo cutoffs

I am attaching an archive with an example of the results for the convergency tests, along with the output of a HF, GW and abs calculations.
I would be very grateful if you could give me any hint on this problem.

Thanks in advance,

[Daniele Varsano]

Dear Aldo,

Altough this may be expected, the optical spectrum has the onset between 2 and 3 eV, peaking around 3-4 eV, where the IP result is closer to the lower ranges and the LF result is closer to the higher values.

please note you are comparing HF and GW electronic quasiparticle gap with optical gap calculated in RPA(IP) calculated on top of KS gap, this is not consistent. In order to be consistent you should either solve the BSE equation on top of GW, or consider the GW energies in the RPA calculations.

About the discrepancies among the two codes it is something you need to investigate and it is probably related with the different used basis set (PW vs gaussians).
I have the impression that your cell is rather small in the z direction for a 2D system, moreover the cut you used is low compared to the cell size in the z direction, it should be just slightly smaller e.g. 28 bohr

Here some suggestions:

1) You performed a spin-polarised calculation and it seems it is not strictly needed, you can do a spin non-polarised calculation saving half of the computational resources.

2) Check the E_xc contribution reported in the yambo report

Code: Select all

E_xc :  -246.4162 [Ry]

this should be equal to the one reported by QE, otherwise it means you are not converged with the number of plane wave, but it should be ok as you considered all the G vectors needed to represent the density.

3) There is no need to check convergence wrt bands in HF calculations as it is obtained summing up the contribution to all occupied states.

4) 10eV vs 3eV it is an enormous difference, you can try to perform an HF calculation using QE (note you will not obtain the same results of Yambo as the QE is a self-consistent HF while Yambo calculated the HF level on top of PBE in a perturbative way) anyway you can check if you find something similar to Crystal (self-consistent), otherwise it means there is some problem in your unit cell setup. I suggest you also check the literature if similar calculations were performed on your system.

Best,
Daniele

[a.ugolotti]

Dear Daniele,

thanks for your reply. Here are some questions/comments on it:

please note you are comparing HF and GW electronic quasiparticle gap with optical gap calculated in RPA(IP) calculated on top of KS gap, this is not consistent.

I understand your concern in the case of the HF result, however I was expecting to obtain, for the case of the spectrum in the RPA-LFE approximation, a bandgap qualitatively similar to that calculated from GW quasiparticle energies as both take into account some degree of correlation (and the perturbation is along the plane of the layer, not perpendicular).

In order to be consistent you should either solve the BSE equation on top of GW, or consider the GW energies in the RPA calculations.

The BSE spectrum is being calculated, however could you please clarify the last part of the sentence? If the RPA calculations you are referring to are those for calculating the spectrum, how can I get the GW energies?

I have the impression that your cell is rather small in the z direction for a 2D system, moreover the cut you used is low compared to the cell size in the z direction, it should be just slightly smaller e.g. 28 bohr

The whole cell is ~28 bohr and the potential was quite flat already at 22 bohrs, hence I assumed the wavefunction was already converged with such cell, while I didn't worry about spurious interactions as I included the cutoff in yambo. However that is something I can carefully test, thanks.

2) Check the E_xc contribution reported in the yambo report

Yes, I did and from the output of QE I get : -269.58085122 Ry. The E_xc values reported by Yambo is connected to the EXXcutoff I chose? Would that be a useful was to check the convergence, other than the value of the bandgap?

4) 10eV vs 3eV it is an enormous difference, you can try to perform an HF calculation using QE (note you will not obtain the same results of Yambo as the QE is a self-consistent HF while Yambo calculated the HF level on top of PBE in a perturbative way) anyway you can check if you find something similar to Crystal (self-consistent), otherwise it means there is some problem in your unit cell setup.

I will try to set up such a calculation, I have not experience in HF-only systems with Espresso; I have already calculated the results with B3LYP functional though and they give a gap of 3-3.6 eV, which agrees to the results available in the literature (which are most commonly PBE or hybrid).

Bests,
Aldo

[Daniele Varsano]

Dear Aldo,

I was expecting to obtain, for the case of the spectrum in the RPA-LFE approximation, a bandgap qualitatively similar to that calculated from GW quasiparticle energies as both take into account some degree of correlation

Electronic, or transport gap (GW) and optical gap (RPA) are two different concepts. If you want to calculate optical gap in Many-Body-Perturbation Theory you need to solve the BSE. Next the BSE and RPA will be probably different as the first take into account e-h interaction while RPA only local field effect. Already in silicon the two spectra are totally different and if they are similar it is for an error cancellation reason. See. e.g. Onida, Rening, Rubio (REVIEWS OF MODERN PHYSICS, VOLUME 74, APRIL 2002).

The BSE spectrum is being calculated, however could you please clarify the last part of the sentence? If the RPA calculations you are referring to are those for calculating the spectrum, how can I get the GW energies?

You need to calculate GW energies then you can plug them in a BSE calculation by using the KfnQPdb variable or as a scissor using KfnQP_E.
If you want to do a GW-RPA calculation you need to plug them by using XfnQPdb or XfnQP_E.

The whole cell is ~28 bohr and the potential was quite flat already at 22 bohrs, hence I assumed the wavefunction was already converged with such cell, while I didn't worry about spurious interactions as I included the cutoff in yambo.

Convergence wrt vacuum in GW is much harder than in DFT where the density only need to be converged. The coulomb cutoff in order to work need to fulfil some geometrical requirements (in particular enough vacuum space needs to be present) , please see Rozzi, Varsano et al PHYSICAL REVIEW B 73, 205119 􏰪2006􏰀 for both points.

Yes, I did and from the output of QE I get : -269.58085122 Ry. The E_xc values reported by Yambo is connected to the EXXcutoff I chose? Would that be a useful was to check the convergence, other than the value of the bandgap?

It is rather different and this is a bad sign. The E_xc value is connected to VXCRLvcs, and this has to be set to the maximum number go G vector (density cutoff). Check carefully this value as the two needs to agree with more accuracy the one you found.

I will try to set up such a calculation, I have not experience in HF-only systems with Espresso;

you need to set input_dft='hf' and pay attention to the divergence treatment (as you did for B3LYP).

I have already calculated the results with B3LYP functional though and they give a gap of 3-3.6 eV, which agrees to the results available in the literature (which are most commonly PBE or hybrid).

Note your PBE gap is smaller than that (1.84 eV).

Best,
Daniele

[a.ugolotti]

Dear Daniele,

It is rather different and this is a bad sign. The E_xc value is connected to VXCRLvcs, and this has to be set to the maximum number go G vector (density cutoff). Check carefully this value as the two needs to agree with more accuracy the one you found.

even keeping the EXX anc VXC cutoffs at maximum, namely the same number of Gvecs I have in the database input, the E_xc contribution is significantly different from QE output :

QE : xc contribution = -269.58092146 Ry
Yambo : E_xc : -246.4161 [Ry]

it might be related to the RIM+cut? That's the only operation I am doing on my system.

Bests,

Aldo

[Daniele Varsano]

Dear Aldo,

RIM+cut do not enter in the evaluation of the xc energy.
Are you using a pseudopotential with non-linear core corrections?
You can check it in the pseudo upf format searching for:

Code: Select all

core_correction="F"

If this is set to T, you need to uncomment/add the following flag in the yambo input:

Code: Select all

UseNLCC 

If this is not the case, we will need to reproduce your error.

Best,
Daniele

[a.ugolotti]

Dear Daniele,

I am using pseudopotentials including non-linear core corrections. Indeed, adding the flag you specified the xc energy is correctly evaluated.
The HF bandgap is however quite unaffected: considering the size of the vacuum region, I can see that 15 A is a little small, but already at 20 A the HF bandgap is practically the same as that at 30 A: ~ 9.37 eV, with the full set of G vecs for EXX and VXC parameters, which is still much grater the the expected size.

Do you have any additional suggestion?

[Daniele Varsano]

Dear Aldo,
the use of the NLCC is discouraged as it is included in the local part of the potential but then it is neglected in the non-local one (e.g. HF), this is the same argument why it is discouraged in QE when dealing with hybrid functionals.

Anyway, when the flag is considered, the value of the local part is calculated correctly and should have been changed wrt the case it was not included (you should check the entity) and hence also the value of the HF and GW gap. BTW, looking again at your old report files you have a GW gap of about 6.2 eV and not 8 eV (see kpoint n. 2).

As suggested before, before passing to GW I would go for an HF calculation with QE and try to understand the discrepancy you found wrt a localized basis set code (if any).

Best,
Daniele

[a.ugolotti]

Dear Daniele,

thank you for all your help.

Let me just try to wrap the topic up, to be sure I got the correct understandig:

1) technical side: non-local core corrections are discouraged, therefore other than handling them correclty by setting up the UseNLCC flag in input, it would be reccomended to re-generate the pseudo without them a check again the results (HF/optical spectra and so-on)

2) theoretical side: IP and RPA-LF spectra are calculated on top of KS energies. Local fields only account for a redistribution of intensities. To re-calculate the energies, the BSE hamiltonian is required, which allows to include the GW energies.

Still, I have two open questions :

1) can the GW energies be given in input to the RPA-LFE case ?
2) is it possible to use the BSE hamiltonian, therefore getting the absorption in eh space, but including only RPA-LFE interactions, as a hartree kernel? Would that correspond to calculating the RPA-LFE spectrum?

Bests,
Aldo

[Daniele Varsano]

Dear Aldo,

1) technical side: non-local core corrections are discouraged, therefore other than handling them correclty by setting up the UseNLCC flag in input, it would be reccomended to re-generate the pseudo without them a check again the results (HF/optical spectra and so-on)

Yes, by using UseNLCC you are consistent with the QE calculations, i.e. you are subtracting correctly the <Vxc> term, but NLCC are anyway not taken into account in the Fock term (as in the QE hybrids calculations), so in order to be consistent it is better to avoid.I think that you will not have major changes, but it is something to verify. If you do not want to generate a pseudopotential the sg15 set should not have NLCC.

2) theoretical side: IP and RPA-LF spectra are calculated on top of KS energies. Local fields only account for a redistribution of intensities. To re-calculate the energies, the BSE hamiltonian is required, which allows to include the GW energies.

Local field affects take into account the disomogeneity of the system, so you also have shift in absorption peak depending indeed by the disomogeneity.
BSE: according to Hedin equation you need to build the excitonic matrix considering QP energies e.g. GW energies, then in the kernel you have the electron-hole attraction which renormalize your sepctrum.

1) can the GW energies be given in input to the RPA-LFE case ?

Yes, by using the following keywords:

Code: Select all

XfnQPdb= "E < ./SAVE/ndb.QP" 

where you include the path of the QP database previously calculated in GW.
or

Code: Select all

% XfnQP_E
scissor value | 1.000000 | 1.000000 |        # [EXTQP Xd] E parameters  (c/v) eV|adim|adim
%

if you want to take into account GW shift by a scissor operator and eventually some stretching of valence and conduction bands.

2) is it possible to use the BSE hamiltonian, therefore getting the absorption in eh space, but including only RPA-LFE interactions, as a hartree kernel? Would that correspond to calculating the RPA-LFE spectrum?

Yes, yon need to build the BSE by incuding Hartree term only, this is done by setting:

Code: Select all

BSENGexx= Ng

where Ng are the same number of block you included in the RPA-LF (NGsBlkXd) and

Code: Select all

BSENGBlk= 0 eV    

here you negelct the e-h attraction.

The two calculation will be equivalent if you consider the same number of bands in BSEBands and BndsRnXd.
Actually in order to be exactly the same you should incude couplng term in the BSE or set the RPA Green function as resonant (GrFnTpXd= "r").
Of course, in order to be the same you need to consider in both either KS or QP energies. QP energies in the BSE are included similarly as wrote before using:
KfnQPdb or KfnQP_E

Best,
Daniele