Point-charge-interaction integrals for each point

[PaulWAyers]

For the electrostatic potential and several other things we are doing, it would be helpful to be able to compute the 1-electron "electron-nuclear" attraction integrals for each point separately. That is, for points $\mathbf{r}_k$ and atomic orbitals (AO basis functions) $\phi_a(\mathbf{r})$, we'd like to be able to evaluate the quantity,

$$ V_{abk} = \int \frac{\phi_a(\mathbf{r}) \phi_b(\mathbf{r})}{|\mathbf{r} - \mathbf{r}_k|} d \mathbf{r} $$

It seems that the "raw" integral we need is GB2AttractionIntegral

horton/horton/gbasis/ints.cpp

Line 198 in e8f9124

GB2AttractionIntegral::GB2AttractionIntegral(long max_shell_type, double* charges,

and that overlay (with sum over "atomic" centers?) is compute_nuclear_attraction.

horton/horton/gbasis/gbasis.cpp

Line 232 in e8f9124

void GOBasis::compute_nuclear_attraction(double* charges, double* centers, long ncharge,

Unfortunately, this part of the code is not (very well) documented. @tovrstra do you have any advice on how to do/expose this? We could call the high-level function one time for every grid point, but that seems wasteful since (I think) a lot of the recursive work gets redone over and over.... @leila-pujal volunteered to delve into whatever implementation we need to add.

gbasis does this too, but we were trying to make it so we can do this with both "old" and "new" HORTON.