Implement support for auxiliary density DFT

[susilehtola]

The input data for LDAs and GGAs can be evaluated faster with the help of Coulomb fitting.

First, instead of evaluating the density in the orbital basis

$$ n({\bf r}) = \sum_{\mu \nu} P_{\mu \nu} \chi_\mu({\bf r}) \chi_{\nu}({\bf r}) $$

one can use the auxiliary basis expansion of the density

$$ n({\bf r}) \approx \sum_\alpha c_\alpha \chi_\alpha({\bf r}) $$

for an update that is linearly scaling in the number of auxiliary functions, instead of quadratic scaling in the orbital functions.

Also the Fock matrix contributions can be evaluated with the help of the auxiliary basis.

An implementation within GauXC would be great, since only a linear amount of data needs to be passed in and out: the density expansion coefficients $c_\alpha$ and the resulting Fock matrix expansion coefficients $f_\alpha$, for each spin.

The algorithm is discussed in Chem. Phys. Lett. 281, 151 (1997) and J. Chem. Phys. 121, 3417 (2004), and is extensively used in e.g. the deMon2k program. It has also been recently used for protein-ligand scoring.

[wavefunction91]

Yea, this is how the CRYSTAL folks do their calcs as well (at least from what I can gather from their paper from the 90s). Easy enough to add, but I'd need a reference calculation to implement. This might generally require access to some of the "guts" of the implementation to ensure the right "knobs" get turned to yield a meaningful comparison.

[susilehtola]

Which paper is this?

[evaleev]

Should be easy to use the existing code to produce reference answer by using DFBS AO basis + a unit (zero exponent) AO as the OBS AO basis and making the AO density matrix only include off-diagonal blocks corresponding to the products of DFBS and unit AO.

[wavefunction91]

@susilehtola Comp. Phys. Comm. 98, 181 (1996) see Eq (10) and surrounding text. They expand the potential in the aux basis, but the principle is the same.

@evaleev That seems reasonable. I'm more curious how much this will save in practice. I think there's definitely an argument on the GPU side as it will decrease memory movement by quite a bit on the density eval side, and if we went ahead and implemented a postiori screening on the potential / Z eval, we could reduce the backend BLAS eval as welll.

[evaleev]

@evaleev That seems reasonable. I'm more curious how much this will save in practice. I think there's definitely an argument on the GPU side as it will decrease memory movement by quite a bit on the density eval side, and if we went ahead and implemented a postiori screening on the potential / Z eval, we could reduce the backend BLAS eval as welll.

the quality of fitted density with regular density fitting basis sets is pretty poor (since we are fitting operator, not density), so I would expect a substantial accuracy hit ... but for movie-making / ML "training" might be OK

[susilehtola]

Thanks; I have had a vague recollection that Laikov was not the first to propose this, yet he fails to cite the preceding work. Towler et al indeed suggest a LDA level scheme which is identical to that of Laikov.

The preceding work is by Zheng and Almlöf, who suggested in Chem. Phys. Lett. 214, 397 (1993) to employ matrix algebra to evaluate density functionals, bypassing numerical integration altogether.

[susilehtola]

the quality of fitted density with regular density fitting basis sets is pretty poor (since we are fitting operator, not density), so I would expect a substantial accuracy hit ... but for movie-making / ML "training" might be OK

The accuracy should be fine since the auxiliary basis set already can describe the Coulomb potential. But yes, I think deMon2k may use larger auxiliary basis sets than the usual ones in quantum chemistry, since they use autogenerated even-tempered sets.

[wavefunction91]

@evaleev @susilehtola I'm fine adding it either way, it's relatively simple to set up as it really only touches one part of the code. However, the more I look at things like the NEO PR, I think the current infra for basis specification is a bit limiting - I think a refactor is likely in order (i.e. specify an OBD,ABS,PBS, etc). I think I have a few ideas.

[susilehtola]

@wavefunction91 while you're at it, you should refactor the code to support arbitrary particles: a formalism analogous to that in NEO can be used to handle positrons, positive and negative muons, etc...

[wavefunction91]

@susilehtola without concrete use cases, I'm probably not going to worry about that at this time. The problem is not "specify an interface for every particle imaginable", it's in the specification of the system of interactions. Right now, you could run any purely fermionic (actually, I suppose that doesn't matter either) system of a single particle type given the proper definition of the XC functional. The challenging bit for NEO is that you also have to evaluate the EPC coupling, which ideally can be implemented by reusing intermediates. Charge/mass don't enter anywhere in GauXC, so as long as we support collections of particles and keep track of their mutual and cross-class coupling interactions, it should just be fine - completely separate and apart from spec of ABS/OBS though, make another issue with concrete use cases for a more general API if it's something you want to see.