A critical bug in the QP equation (g0w0_eigvals) #36
Description
Hello @kovalp @mbarbry
My recent benchmark on QP spectra for several medium-size molecules shows a big difference between the value of Sig_c calculated in Pyscf vs. molgw. As you can find in the figure below, everything in mf level is identical, both HF eigenvalues and Fock exchange energies. But Sig_c significantly differs especially for states above fermi([47]), leading to different QP spectra. This is clearly a problem in pyscf must be resolved.
You can find input/output files of the above plots for anthracene here. Several points that I should add:
1- both molgw and pyscf use the same basis. however, molgw benefits from an auxiliary one which shouldn't be the origin on such a big discrepancy.
2- pyscf calculations are done in a reasonable range: nocc=15 and nvrt=15
3- all keys which are recently added to W_c in gw_iter class are False. i.e. initguess, limited_bnd, duplicate,...
4- I doublechecked the matrix element of W_c and confirm: get_snmw2sf=get_snmw2sf_iter. Also gw_corr_res_iter()= gw_corr_res().
5- pyscf has not achieved convergency even after 1000 iterations due to these wrongly calculated energies. Even the first iteration /pyscf/pyscf.out clearly shows our code is sick:
Iteration #1 Relative Error: 0.1519962
Spin1: [-14.54872 -14.4106 -13.81212 -13.56293 -13.22433 -13.22686 -12.33767 -12.18217 -11.84294 -11.82644 -10.44084 -10.25039 -9.14067 -8.50937 -7.24384 0.05593 1.69312 2.45949 -5.95737 3.89415 16.55898 6.98528 -2.52392 17.60815 4.87657 8.02117 6.16321 9.66957 -15.02137 38.41561]
5- Earlier I had benchmarked EA and IE of 42 small molecules and never face with such a huge difference. It seems this bug won't show up effectively in the small calculations. You can find a recent benchmark for CH3 using gw_iter in the small_mol folder and figure below:
I spent some time on this issue and finally found the issue.
In both GW and GW_iter, iterating the QP-equation def g0w0_eigvals(self) sums the correction obtained from integral (sn2i) and residue (sn2r) parts plus mf eigenenergies (get_h0_vh_x_expval) as follows:
for s,(evhf,n2i,n2r,nn) in enumerate(zip(self.get_h0_vh_x_expval,sn2i,sn2r,self.nn)):
sn2eval_gw.append(evhf[nn]+n2i+n2r)
Surprisingly, get_h0_vh_x_expval for big systems makes a huge error, while it makes an error in the range of 1e-04 (Ha) for small systems. The following figure shows the benchmark of H_0 components for Anthraene w.r.t MOLGW reference.
Briefly (get_h0_vh_x_expval != mo_energy[0]). Moreover, one can see that the calculated Hartree for high energies (above state [80]) is different wrt MOLGW while Fock is fine.
In pyscf_master, mo_energy,mo_coeff=eig(get_hcore(mol), get_ovlp(mol)) based on the standard Roothaan approximation, i.e. HC = SCE, not like ours which is too simplistic. Moreover, in get_h0_vh_x_expval the non-diagonal elements are ignored and just diagonal elements are taken into account by einsum. This is why it works for small systems with low overlaps, but for medium-sized, it doesn't.
Therefore, just by substituting get_h0_vh_x_expval by mo_energy[0] or ksn2e[0] in above-mentioned code, everything will be satisfying w.r.t. MOLGW's QP and convergence has been reached at 6th iteration. as you can find in figure below.
I will commit these changes in addition few typos if you both agree!