Add support for mGGA, SAP guess, orbital rotation gradient, DIIS mixer

[alexmaryewski]

This PR contains multiple features that had to be implemented simultaneously in order to make mGGA work.

These include:

• General support for mGGA functionals is added; 6 new functionals are implemented (TPSS, SCAN, r2SCAN, r4SCAN, TASK, TASK+CC).
• Superposition of atomic potentials (SAP) guess is implemented and is now the default.
• A single convergence criterion is replaced by two: SCF total energy change and occupied-virtual orbital rotation gradient norm.

All tests except B3LYP pass: in the latter, N-N and C-N comparison fails; this requires attention from a knowledgeable person.

Additional thanks go to @vanderhe for helping with stylistic issues.

[vanderhe]

Additional thoughts supplementing alexmaryewski#1:

• I would prefer not to change the convergence criterion for the existing functionals
• I would prefer not to expose the mixer choice to the user, if not absolutely necessary

To do:

• discuss the initial guess, as it is highly desirable to not having to dump the SAP data into the code base
• either correctly implement or block the ZORA code path for MGGAs
• add regression tests

[alexmaryewski]

Replying to @vanderhe points:

1. The mixed quantity was not changed in this PR; I now think the potential is the optimal choice for us.
2. @bhourahine, what do you think?
3. During my experiments with mGGA, I found it necessary to be able to switch to simple mixer in certain problematic cases, like hydrogen atom with numerically ill-behaved functionals.

Speaking of SAP: would you prefer it to be read dynamically from a data file?

Regression tests will be added as soon as I confirm that the two-centre part works as intended.

[vanderhe]

Replying to @vanderhe points:

1. The mixed quantity was not changed in this PR; I now think the potential is the optimal choice for us.

My bad, sorry.

[vanderhe]

Regarding the B3LYP testcase:
Unfortunately the HOAO sometimes indeed is sensitive w.r.t. the initial potential guess, which depending on the situation results in Hubbard $U$ deviations (~1e-06 a.u.) that exceed the testing tolerance. Checking the convergence of the eigenspectrum might help, but there is of course no guarantee that both starting points converge to the same eigenspectrum...

[alexmaryewski]

@vanderhe FYI

I added 446976a to see if putting a check on convergence of eigenvalues would solve this. Sadly, it seems like eigenvalues don't want to converge to better than 1e-08 even for the most stable functionals, e.g. PBE.

[vanderhe]

@vanderhe FYI

I added 446976a to see if putting a check on convergence of eigenvalues would solve this. Sadly, it seems like eigenvalues don't want to converge to better than 1e-08 even for the most stable functionals, e.g. PBE.

@alexmaryewski I would expect the unoccupied part of the eigenspectrum to be numerically unstable, hence suggest to restrict the comparison to the occupied states only.

[alexmaryewski]

@vanderhe, thanks for suggestion. I just updated the subroutine to check eigenvalue difference only for the occupied part of the spectrum. This does solve the problem of convergence, yet does not fix the B3LYP case.

[vanderhe]

@alexmaryewski Coming back to the convergence of the MGGAs: I think we should plot the Thomas-Fermi and SAP potential for one of the pathological cases (e.g. spin-pol. H). There must be some qualitative feature that the SCF iterations cannot restore (?), especially if we consider that SAP also performs well for the compressed runs.

[aradi]

Just a quick remark on this PR. As far as I can see, it contains a lot of different not necessarily strictly related changes. That is fine for a draft PR. But they should go into individual PRs at the very end before being merged into the main code.