[WIP] State refactor

[Technici4n]

No description provided.

[Technici4n]

For exact exchange, would the psi go into the densities as well?

[antoine-levitt]

We discussed this with @mfherbst. Ideally the codensities would go there, but they're probably too expensive to store. Maybe we could do a kind of lazy structure...

[antoine-levitt]

Actually we probably should think about this a bit more because it violates the assumption that the "direct arrow" is quadratic

[Technici4n]

I was thinking we might keep that quadratic assumption for direct arrows and instead go for an indirect arrow:

orbitals -> densities.orbitals -> E_x_exact

(with the first map being the identity)

It should allow for initial guessing, mixing, preconditioning, and saving (expensive but necessary for exact exchange?) of the orbitals. But I don't know if it's the best fit. For example for response you'd get both an orbital change hiding as a density change and a "normal" orbital change from the response solver.

[antoine-levitt]

Actually I'm not even sure about codensities I think I was wrong there. Having orbitals in densities kind of defeats the whole purpose of the thing. We should maybe just drop the quadratic assumption, or have three kinds of terms, quadraticterm, densitiesterm and "misc" terms.

[antoine-levitt]

I think hybrids are sufficiently different (and rate limiting, since they're the most expensive part by far) that they have to be special-cased anyway

[Technici4n]

Alright, pushed a big update! It's still quite messy and a lot of functionality is likely broken, but at least a basic SCF works for silicon, for Hubbard and for meta-GGAs 😉

Note 1: Having a Hubbard potential is quite unusual.
Note 2: I haven't touched other solvers such as potential mixing yet, but note that mixing the potentials will be a lot cleaner, as there won't be a need for hacks such as total_local_potential.

[antoine-levitt]

Note 1: Having a Hubbard potential is quite unusual.

Really? It should be something like a Fock matrix.

Note 2: I haven't touched other solvers such as potential mixing yet, but note that mixing the potentials will be a lot cleaner, as there won't be a need for hacks such as total_local_potential.

That's a nice side effect!

[Technici4n]

Really? It should be something like a Fock matrix.

Here it's just 1/T(2) * term.U * (I - 2n). A bit unusual maybe 😄

But it's probably the right expression if we want to mix the hubbard_n (an interesting direction!)

[antoine-levitt]

Yeah I want to look at hubbard mixing more closely, I don't know what people do currently.

[mfherbst]

Yeah I want to look at hubbard mixing more closely, I don't know what people do currently.

I agree, we briefly discussed this here, too. At least QuantumEspresso seems to do nothing (i.e. just rebuild every iteration with the current state)

update: see comment by @jaemolihm below !

[jaemolihm]

Hi, Quantum ESPRESSO does mix the Hubbard occupations. If I understand correctly, it essentially concatenates the charge density at each real-space grid and the Hubbard occupations (in a TYPE mix_type), and applies Broyden mixing.

You can take a look at https://gitlab.com/QEF/q-e/-/blob/master/PW/src/mix_rho.f90 and search for lda_plus_u.

But I am not aware of any rigorous analysis of the convergence or stability.

[antoine-levitt]

Thanks! Do you know if they precondition it in any way?

[jaemolihm]

It seems preconditioning is applied only to the charge density, but not to the Hubbard occupation.

[Technici4n]

A related question is mixing and preconditioning for meta-GGAs. (May or may not be more relevant than DFT+U.) I think that QE might be missing meta-GGA DFPT, for example. (DFPT, preconditioning, convergence analysis... all are related :))

[Technici4n]

Some new thoughts: the state refactor is purely a 2P (two-point) concept, to help us keep track of various 2P quantities. Contrast that to hybrids or other orbital-dependent functionals which need the full 4P picture in the SCF, and where the 2P picture is an internal detail.

[antoine-levitt]

Some new thoughts: the state refactor is purely a 2P (two-point) concept, to help us keep track of various 2P quantities. Contrast that to hybrids or other orbital-dependent functionals which need the full 4P picture in the SCF, and where the 2P picture is an internal detail.

By that you mean 1P and 2P respectively? (2P and 4P correspond to operators between 1P and 2P quantities).

As done currently, I agree with what you said and in this formulation I agree it's orthogonal to hybrids.

My question is if we can generalize the notion of state to include "stuff computable from the density matrix and storable in memory that is needed to build the hamiltonian". Eg for +U the occupation matrices fit this description. For hybrids, the set of all codensities would be nice to have there, but it's (probably) too big to store. I have to read the pcdiis and ace papers in some more detail to see if the projected DM in the set of current orbitals or something like that can go into the notion of state.

[mfherbst]

I have to read the pcdiis and ace papers in some more detail to see if the projected DM in the set of current orbitals or something like that can go into the notion of state.

I think what happens here is that we have the concept of an exact representation of the states (which is e.g. all orbitals or all codensities for hybrid) or an quasi-exact representation of states (which ACE-like compressions or the PCdiis ideas can yield).

I say quasi-exact here because these compressed storage methods rely on projections into the exact "occupied" orbitals to be lossless. So if we say the contract is that the state object captures the state to an accuracy similar to the accuracy of the psis you have computed, than I think this is feasible. If we want more than than (i.e. true lossless state from which a calculation can be continued later without harm) that may be trickier.

[Technici4n]

By that you mean 1P and 2P respectively? (2P and 4P correspond to operators between 1P and 2P quantities).

Indeed yes.

For hybrids, the set of all codensities would be nice to have there, but it's (probably) too big to store. I have to read the pcdiis and ace papers in some more detail to see if the projected DM in the set of current orbitals or something like that can go into the notion of state.

Does the orbital representation of the density matrix (i.e. psi and occupations) not suffice? If you are willing to perform the SCF over that I don't see a conceptual reason to store anything else, hence the thought of having a dedicated 2P SCF.

(warning: rambling)

1P SCF

The density-to-potential map (or a generalization thereof for mgga, hubbard, ...) is entirely 2P, whereas the potential-to-density map is internally always a 4P operator.

2P SCF

Now we embrace a 4P "density-to-potential" (in reality: density-matrix-to-hamiltonian) map, and we use the 4P potential-to-density map directly. In this picture the 1P terms become an implementation detail of the 4P density-to-potential map. So for example you would lift the 2P density-to-XC-potential map to a 4P equivalent.

I think it makes sense to have the two as separate algorithms. After all the meat of the computation (the 4P potential-to-density map) is shared between them. Very interestingly, it should be possible to take the 2P SCF and make a 2P DFPT out of it which could work for hybrids.

an quasi-exact representation of states (which ACE-like compressions

ACE should be an implementation detail of the density-matrix-to-hamiltonian map, no? It sounds like it should not be part of the state.

[mfherbst]

Does the orbital representation of the density matrix (i.e. psi and occupations) not suffice?

and the related

ACE should be an implementation detail of the density-matrix-to-hamiltonian map, no? It sounds like it should not be part of the state.

My thinking is along our discussion to split up this map into two steps density matrix -> state -> hamiltonian.

In that setting for both your remarks it depends how you think of the state, i.e. what you want to be able to do with it. In our current philosophy (where all is one step) you are right, but since our whole goal is to insert these "state" intermediates it is not so clear any more at which of the two mappings compressions like the ACE should be applied. E.g. if the state is the object we base Anderson on we potentially have quite a few of them floating around in memory (to make a linear combination and from that construct the Hamiltonian), so storage may become a concern and compressing in the first map is thus a necessity.