First draft of Pcdiis acceleration for SCF iterations including exact exchange

[j-mangold]

Summary

Adds a rudimentary implementation of the Pcdiis accelerator for SCF iterations including exact exchange.

Pcdiis stands for Projected-Commutator-DIIS and is an attempt to make CDIIS feasible for large basis sets.
Instead of computing e = [H,ρ] in the full basis, this is only done in a subspace: ē = <ψ_ref|e|ψ_ref>
Instead of mixing full density or Fock matrices, the occupied states are mixed after gauge fixing (gf): ψ_gf = |ψ_occ><ψ_occ|ψ_ref_occ>

Solving for the coefficients is done just like for the original DIIS:

┌ ┐ ┌ ┐ ┌ ┐
│ 0 1 1 … 1 │ │ λ │ │ 1 │
│ 1 B₁₁ B₁₂ … B₁ₙ│ │c₁ │ │ 0 │
│ 1 B₂₁ B₂₂ … B₂ₙ│ · │c₂ │ = │ 0 │
│ ⋮ ⋮ ⋮ ⋱ ⋮ │ │⋮│ │ ⋮│
│ 1 Bₙ₁ Bₙ₂ … Bₙₙ│ │cₙ│ │ 0 │
└ ┘ └ ┘ └ ┘

where Bᵢⱼ = ⟨eᵢ|eⱼ⟩

I copied the layout of the anderson.jl file for the implementation.

For now, there is no history management implemented except for a maximal depth.

##Changes

• Added src/scf/pcdiis.jl containing the Pcdiis implementation
• Added examples/pcdiis_acceleration.jl for demonstration
• Modified src/scf/scf_solvers.jl to add corresponding solver function
• Modified src/DFTK.jl to include src/scf/pcdiis.jl

[antoine-levitt]

Thanks, looking nice! I'd like to try and merge this discussion with that of #1142, so we don't actually special case exact exchange, meta-GGA, +U, etc (which goes against the DFTK spirit of being relatively agnostic to what precise terms are in the system and have a clear terms API) but rather have DFTK ask the terms what they think is the minimal state they need to keep around to build a hamiltonian and do mixing, and then operate generically on this "state". I have to read up more on PCDIIS and ACE to see what the terms boundary should be, and if this is even possible. Maybe let's set up a zoom call with all interested parties (@Technici4n , @mfherbst @toschaefer ... ?)

[j-mangold]

Thanks! I'd be happy to discuss with you! Do you need anything from me to merge this with #1142, or should we talk before that anyway?
Best regards,
Jona

[antoine-levitt]

The state stuff is more of a long term project. Since it doesn't actually change any "core" DFTK file and it's of interest to you I'm fine with merging as is (maybe just move them to a separate exchange folder?) Or maybe it's simpler for the code to just live outside of DFTK (you can make a separate package, or just a repo with the files, it's quite easy to do) until we figure out the "right" way to do it.

[j-mangold]

Then it's maybe really the easiest to leave it outside for now. I'll set up a repo.

[toschaefer]

Thanks for you feedback, Antoine! I would actually prefer to take your offer and merge this (into a separate exchange folder or so...). We really want to make a performant EXX available to everyone quickly, without complicated external packages. Now that EXX is available in DFTK, it would be a shame if a good mixer had to wait outside the code.

Regarding the refactoring later to fit the long-term API project in #1142, this will be just peanuts with AI coding agents and we promise to do it as soon as such an API is ready.

[Technici4n]

For now I would just keep the state refactor #1142 in mind as a nice conceptual tool and not focus on getting it implemented. Hybrids are a big question mark with the state refactor, and I think we should continue to merge more hybrids developments. That will provide us with enough information to successfully perform the state refactor in the future.

[Technici4n]

I forgot if this executed or not in the tests like this.

[mfherbst]

I think yes.

In any case I would suggest to merge this into the HartreeFock examples and tests we already have.

[j-mangold]

I replaced this example with a cleaner one. I think there was a problem with an unnecessary package (NPZ) being used. I am sorry!

[Technici4n]

Doing this and inserting it in a density-based SCF loop is IMO quite a hack. I think it would be preferable to introduce a new function that serves as an alternative self_consistent_field, for example orbital_mixing_scf (see also the potential_mixing_scf) that fully embraces $\psi$ as the main variable and does not attempt to perform any mixing on $\rho$, but rather recomputes it from the new $\psi$. In other words I am suggesting a proper 4P SCF (4 point i.e. density-matrix-to-density-matrix, of course using an orbital representation).

More work, but also more forward-thinking I think. This would then be a nice playground for various types of orbital-dependent XC functionals, and a nice alternative to hubbard_n / $\tau$ mixing for +U and meta-GGAs.

[mfherbst]

While I think the current implementation of @j-mangold is indeed good for a first exploration, I agree with @Technici4n that the cleaner solution would be to have an scf_orbital_mixing (use this order of the name segments to follow our current naming convention) where the occupations and the orbitals are the prime state variables. After the state refactor we may merge all these various SCFs, but for now I think indeed the best solution to have two.

So my concrete suggestion is to not have the scf_pcdiis_solver function, but introduce a scf_orbital_mixing where one directly uses the PcdiisAcceleration (a bit like it's done with the scf_potential_mixing. I would also not be surprised if this gives some opportunities to share more code with AndersonAcceleration or even generalise the current AndersonAcceleration implementation to cover the PcdiisAcceleration case as well.

But one step at a time. If @j-mangold and @toschaefer you want to work towards getting this into DFTK asap, I'd suggest the first step is to draft an orbital-based SCF. Let me know if I can be of help for that (e.g. for a pair coding session or similar).

[j-mangold]

I do not have a strong opinion on the further approach here. I suggest to talk about it in our meeting next week and then decide how to continue.

[mfherbst]

Thanks @j-mangold. Good to see you have managed to get a first working version so quickly. Now it's all about massaging it and bringing it into a form it can live as part of the main code.

[mfherbst]

I think yes.

In any case I would suggest to merge this into the HartreeFock examples and tests we already have.

[mfherbst]

While I think the current implementation of @j-mangold is indeed good for a first exploration, I agree with @Technici4n that the cleaner solution would be to have an scf_orbital_mixing (use this order of the name segments to follow our current naming convention) where the occupations and the orbitals are the prime state variables. After the state refactor we may merge all these various SCFs, but for now I think indeed the best solution to have two.

So my concrete suggestion is to not have the scf_pcdiis_solver function, but introduce a scf_orbital_mixing where one directly uses the PcdiisAcceleration (a bit like it's done with the scf_potential_mixing. I would also not be surprised if this gives some opportunities to share more code with AndersonAcceleration or even generalise the current AndersonAcceleration implementation to cover the PcdiisAcceleration case as well.

But one step at a time. If @j-mangold and @toschaefer you want to work towards getting this into DFTK asap, I'd suggest the first step is to draft an orbital-based SCF. Let me know if I can be of help for that (e.g. for a pair coding session or similar).