[xtp] QPsolver for CDA and grid inefficiency

[baumeier]

TLDR
We should be more clever with the QP solver, including:

• check for stability of solution (grid/fixedpoint)
• combine that to more flexibility in the grid: (i) extending it, if no stable point is found; (ii) stopping the grid evaluation if a stable point is found
• be a bit more clever with the fixed point methods for the CDA, and only revert to an (improved) grid search if it fails. Are there algorithms for multiple fixed points?

The default method for solving the QP equations is a grid search, followed by bisection. As default, the grid is constructed with 1000 points covering the interval of [-0.5 : +0.5] Hartree around the input QP energy.

General issue:
The range covered by the interval is too small and the reported QP energy is not a stable fixed-point of Sigma_c(w) +(Sigma_x + w_0 -V_xc) = w.
Example: H2O core level

... ...  Roots found for qplevel:0 (qpenergy:qpweight)
                -19.105:-0.0012064 -19.103:0.0063074 Root chosen -19.103

The QP weight is defined as Z=(1-dSigma_c(w)/dw)^-1 and for "the" solution should be on the order of 1.
This can also be considered from the criterion of stability for a fixed point, |dSigma_c(w)/dw| < 1, which would translate to 1/2 < Z < infty. However, there is also evidence that dSigma_c(w)/dw < 0, so 1/2 < Z < 1.
In the example above, Z~10^-3 for both found solutions, so we can conclude that none of them is what we are looking for, which includes the reported solution.

Possible solution:
Add a check for stability (e.g., Z>0.3), and if none of the found roots is stable, extend the grid search.

Specific issue for the CDA
The grid search becomes extremely costly for the CDA. As a rough estimate, the dielectric matrix has to be computed with the RPA, and then inverted, at each grid point for as many poles as are included in the contour. Even for the HOMO, this means way more than one pole, when the grid includes frequencies up to -0.5 Hartree. This escalates quickly when the molecules get bigger, and the spectrum denser.
Example: HOMO value for MADN (59 atoms)

HOMO & LUMO from CDA G0W0 on grid
... ...   HOMO  =  116 DFT = -0.2039 VXC = -0.3982 S-X = -0.4312 S-C = -0.0076 GWA = -0.2445
... ...   LUMO  =  117 DFT = -0.0644 VXC = -0.3892 S-X = -0.2597 S-C = -0.0868 GWA = -0.0217
Time for preparation of 100 dielectric matrices for numerical integration (56 threads): ~30m
Time for QPeq solution on grid for *2* levels (56 threads): *~29hours*

From the timing, it seems that around 6000 matrix calculations and inversions (dimension 3555) need to be calculated. This is only for two QP states.

Fixed point methods
For comparison, I played around with different accelerated fixed point methods (Aitken's Delta Squared Process, and Anderson mixing).
Example: HOMO and LUMO for MADN (with exact method)

.. ... 2020-9-6 21:41:8 Solving QP equations  Anderson -0.228594
 Anderson -0.245539
 Anderson -0.245445
 Anderson -0.245446
Found stable QP at -0.245445 Hrtr with weight 0.836778
 Anderson -0.0392835
 Anderson -0.0216834
 Anderson -0.0217725
 Anderson -0.021772
Found stable QP at -0.0217725 Hrtr with weight 0.84345

Basically four steps are enough here. Note that this one is already checking the QP weight, if this is a stable solution.

The problem with these methods is of course that

• convergence can be difficult even with accelerators
• they can find unstable fixed points

I played around with reinitialized the fixed point procedure in case an unstable FP is found. This is at this point a bit arbitrary, it worked for the H2O case mentioned above, but isn't guaranteed to work always.

Thoughts? @gtirimbo @JensWehner @felipeZ

[JensWehner]

So, yes the rootfinding is not super. I like your approaches. I do not understand what the accelerators do, we also implemented the Newton-Method for rootfinding, which converges quickly,but only locally. Basically you have to find all the roots of a very unpleasant function and a local scheme will not help you there. Additionally the fact that our eta is not zero, but has to be a finite number is also makes the evaluation of the quasiparticle weight more difficult.

When I implemented the current method I was not happy with the runtime but it was the most robust at the moment, also because for ppm and exact the individual evaluation is pretty cheap and only the setup is expensive.

Specific question:

What should you do if no stable solution with Z>0.3 is found?
What do you do if Z is inf, how do you protect against a nan there? (Because of eta in ppm and exact, Z can take basically all possible values)
This is used very deep inside the scGW, how do you properly blackbox it?
If you use a lower bound for Z, which one? Typcially low Z indicate that physically the QP approximation failed.

What exactly do you mean by stable fixed points?

[baumeier]

Newton Method and CDA is again a bit meh, since we need to take the derivative numerically. Probably still better than going the full grid but this is why I was looking at some of these derivative-free methods.

It is true that for ppm and exact the current method is not that much of an issue (besides the fixed range maybe). The CDA suffers massively. (no free lunches)

What should you do if no stable solution with Z>0.3 is found?

Depends what method we are talking about. Let's assume we run a grid solver and find nothing that fits in the prescribed range. Range can be extended but I would say only up to a certain physically sensible (yet to be defined) amount. If any fixedpoint method does not yield anything after N attempts, let it fall back to grid as is done now.

What do you do if Z is inf, how do you protect against a nan there? (Because of eta in ppm and exact, Z can take basically all possible values)

Do you mean Z infinite or dSigma/dw? I don't think the former will happen, the latter could be large but then Z will be zero.

This is used very deep inside the scGW, how do you properly blackbox it?

If it "works" for G0W0, which is a special case of evGW, isn't that blackboxing enough?

If you use a lower bound for Z, which one? Typcially low Z indicate that physically the QP approximation failed.

Not sure at which point a low Z means QP failure. As I said the choices above are very adhoc.

What exactly do you mean by stable fixed points?

That's maybe a bit handwavy and inadequately borrowed from proper mathematics.

[JensWehner]

Do you mean Z infinite or dSigma/dw? I don't think the former will happen, the latter could be large but then Z will be zero.

Z which is zero if dSigma/dw is close to 1. I think dSigma/dw is a better criterion.

If it "works" for G0W0, which is a special case of evGW, isn't that blackboxing enough?

true

Not sure at which point a low Z means QP failure. As I said the choices above are very adhoc.

Yeah that is the part with blackboxing, where I do not know what to do. E.g. in scGW iteration you might not care too much about it.