rigw_contour vs rigw vs GW

[martijn]

Hi,

I have run G0W0 calculations on top a B3LYP/def2-TZVPP reference for a (MgO)4 cluster using gw, rigw (ri-ac-gw) and rigw_contour (ri-cd-gw), as implemented in Turbomole 7.4.1 and observe something I was not expecting:

gw
  10b1      -5.992    -7.807   -22.642   -24.913     2.272   -20.386     0.805     -0.243   
  10b3      -5.992    -7.807   -22.642   -24.913     2.272   -20.386     0.805     -0.243   
  10b2      -5.992    -7.807   -22.642   -24.913     2.272   -20.386     0.805     -0.243   
 ----------------------------------------------------------------------------------------
  11a       -2.060    -0.518    -3.280    -2.385    -0.894    -4.930     0.934     -0.702E-01

ri-ac-gw
  10b1      -5.992    -7.823   -22.642   -24.913     2.272   -20.386     0.812     -0.231   
  10b3      -5.992    -7.823   -22.642   -24.913     2.272   -20.386     0.812     -0.231   
  10b2      -5.992    -7.823   -22.642   -24.913     2.272   -20.386     0.812     -0.231   
 ----------------------------------------------------------------------------------------
  11a       -2.060    -0.513    -3.280    -2.385    -0.894    -4.930     0.938     -0.666E-01

ri-cd-gw
  10b1      -5.992    -8.247   -22.642   -24.913     2.272   -20.386     1.000       0.00   
  10b3      -5.992    -8.247   -22.642   -24.913     2.272   -20.386     1.000       0.00   
  10b2      -5.992    -8.247   -22.642   -24.913     2.272   -20.386     1.000       0.00   
 ----------------------------------------------------------------------------------------
  11a       -2.060    -0.410    -3.280    -2.385    -0.894    -4.930     1.000       0.00

I would have expected ri-cd-gw to lie as close or closer than ri-ac-gw to gw but in contrast ri-cd-gw results lie further away than ri-ac-gw, which in it self lies quite close to the gw results. I am also surprised that for ri-cd-gw the z column is exactly 1 and the dS/de column is exactly 0 for all quasiparticle levels. Was I wrong in my expectations and is this the expected behavior and/or am I missing something?

Thanks,

Martijn

[chris.hol]

Dear Martijn,

the difference between ri-cd-gw and gw/ri-ac-gw is mainly that ri-cd-gw does not use linearization in G0W0. Therefore, all z's are 1. You calculations only differ by this linearization as far as I can see.
To use ri-cd-gw for G0W0 it is recommended to iteratively solve the quasiparticle equation (by adding for example qpeiter 10 to the $rigw keyword group) - note that this, unlike for ri-ac-gw or linearization in general, will also work properly for non-valence orbitals (though the results may be different to linearized ones, especially if linearization is not well-suited). Further ri-cd-gw is the recommended method for the iterative evGW approach when the molecule is too large for full GW ($gw). But for IP/EA G0W0 (so HOMO/LUMO) I just recommend to use ri-ac-gw, as it is the fastest by far for large molecules - also quite some effort was invested into yielding good linearized values in this case (we employ a unique double-grid double-Pade expansion for this case, which is simply theoretically not compatible with ri-cd-gw).

Concerning the accuracy of ri-ac-gw and ri-cd-gw in general:

• For HOMO and LUMO in linearized G0W0: ri-ac-gw will in nearly any case yield quasiparticle energies VERY close to those of full G0W0. Therefore it is brilliant if one only needs to know the IP/EA, but not much else. This is often the case though.
• Beyond HOMO and LUMO the quality of ri-ac-gw will degrade significantly, as does the quality of linearized results. 
• ri-cd-gw will yield correlation self-energies (Sigma_C) which are always VERY close to the of full GW. This holds for valence and core orbitals, making it the best suited for itertive approaches as evGW or G0W0 where quasiparticle energies are solved iteratively. Basically if one treats all orbitals within ri-cd-gw and an iterative treatment the results will be nearly identical to those of the corresponsing full $gw approach. However, the trick is of course to limit the orbital range for making calculations feasible.

One final remark: You molecule features degenerate orbitals. This is the single case where more integration points (npoints 256 or 512) really helps in improving the results of ri-ac-gw and ri-cd-gw. This behaviour originates from the sharp step in the Green's function at the poles, were the description actually suffers from numerical noise when they become degenerate.
Otherwise the default of 128 integration points will yield deviations between gw & ri-ac-gw for HOMO and LUMO that are usually well below 0.01 eV. Though the deviation of 0.015 eV between gw & ri-ac-gw in your calculation is probably also not critical.

All the best,
Christof

[martijn]

Hi Christof,

Thank you very much! I think I have figured out the origin of the difference between ri-cd-gw and GW/ri-ac-gw, I hadn't realized I had to set qpeiter (and at least define in 7.4.1 seems to not set it automatically when switching on contour). The results above were for the first iteration, after 10 iteration the ri-cd-gw results lie very close to their GW/ri-ac-gw counterparts:

  10b1      -5.992    -7.756   -22.150   -24.913     2.763   -20.386     1.000       0.00   
  10b3      -5.992    -7.756   -22.150   -24.913     2.763   -20.386     1.000       0.00   
  10b2      -5.992    -7.756   -22.150   -24.913     2.763   -20.386     1.000       0.00   
 ----------------------------------------------------------------------------------------

As an aside how do you calculate unoccupied QP energies with ri-cd-gw, when not using gap/ips+? While gw/ri-ac-gw seem to by default give the lowest 5 GW unoccupied QP energies, ri-cd-gw, assuming I haven't done something (else) wrong, only calculates occupied QP-energies with GW. It prints values for all possible unoccupied orbital energies but with values equaling the corresponding KS energies.

Best,

Martijn

[chris.hol]

Hi Martijn,

for ri-cd-gw an enhanced selection may be used, where the orbital range may be defined (note that this does NOT work for ri-ac-gw, as we strictly recommend the latter only for HOMO and LUMO):

$rigw
 npoints 256
 contour start=10 end=11 irrep=1 ! a
 contour start=10 end=11 irrep=2 ! b1
 contour start=10 end=11 irrep=3 ! b2
 contour start=10 end=11 irrep=4 ! b3

This should get you the highest occupied and lowest unoccupied states for each irrep in you calculation. If you just want HOMO/LUMO you may also specify:

$rigw
 npoints 256
 contour start=11 end=11 irrep=1 ! a
 contour start=10 end=10 irrep=2 ! b1
 contour start=10 end=10 irrep=3 ! b2
 contour start=10 end=10 irrep=4 ! b3

In an open shell calculation further nspin=1 or nspin=2 can be added to each line. The number refers to the number reported in the GW output (and is consistent with that of "eiger") Note that only 1 specification line per irrep and spin is currently possible. Multiple definitions of the same IRREP will lead to only the last one being applied.

All the best,
Christof

P.S.: We tried to align the ips/gap keywords in TM 7.5 between the different methods (and improve the logics behind incorporating gaps). However, especially when the molecule belongs to a point group with degenerate states, a manual choice is sometimes inevitable...

[martijn]

Thanks Christof! That was very useful. Much appreciated.

Maybe it would be good if the need for qpeiter in the case of ri-cd-gw calculations is also mentioned in section 14.1.3 of the manual and/or if define would set say 'qpeiter 10" automatically when the contour option is turned on (this might be the case in 7.5 already, haven't had a chance to try it yet).

Btw, and forgive more for my ignorance, is it correct that the effect of using qpeiter is small/negligible for the other GW implementations?

Thanks again,

Martijn

[chris.hol]

Dear Martijn,

setting qpeiter to some default value in define when ri-cd-gw is definitely a good idea, I will do so in the next release (it is unfortunately too late for TM 7.5).

For HOMO/LUMO the qpeiter option is usually very close to the linearized results and therefore for other GW implementations the effect is often really small or negligible. This may no longer be the case for core/non-valence states, were linearization may lead to quite off results sometimes - in these case qpeiter will be closer to the "true" quasiparticle energy.

All the best,
Christof