Evaluation of density variables based on orbital coefficients

[susilehtola]

As far as I understand, GauXC currently operates by evaluating the electron density on the grid from the density matrix

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

This approach gets high FLOPs, since it can be formulated with efficient intermediates: first compute the matrix multiplication $p_\mu({\bf r}) = P_{\mu \nu} \chi_\nu({\bf r})$ and then get the electron density from $n({\bf r}) = \sum_\mu p_\mu ({\bf r}) \chi_\mu({\bf r})$ with $N_\text{AO}^2 N_\text{grid}$ operations.

However, when you have a large basis set and few occupied orbitals (extreme case is Perdew-Zunger self-interaction correction where you need to evaluate Fock matrices for single occupied orbitals), you can compute the density faster from

$$ n({\bf r}) = \sum_i f_i |\psi_i ({\bf r})|^2 = \sum_i f_i \left|\sum_\mu C_{\mu i} \chi_\mu({\bf r})\right|^2. $$

Evaluating the orbitals takes $N_\text{MO} N_\text{AO} N_\text{grid}$ effort, which is the rate determining step. We therefore save a factor of $N_\text{MO}/N_\text{AO}$ operations. One gets the speedup in a dense basis set, where sparsity is not significant.

[wavefunction91]

This, in fact, brings up two points that require careful consideration and a likely refactor of the input interface:

1. The stated change of density -> MO input. This will require some thought for two-component bases, but should be relatively straight forward
2. The inclusion of flexible occupation specification ($f_i$) - I can see a lot of potential applications here (finite temp, etc).

Practically, these changes are easy to make, but getting the interface right (w/o replicating a ton of code) will likely be a bit challenging. N.B. #76 also requires some thought about how we specify inputs, so it probably warrants a full refactor of the interface.