Maximilian Scheurer

Just an additional comment to the feature request: To do this with TD-DFT in practice, we'd also need excited state gradients 😊

For g2g, maybe just [this](https://github.com/psi4/psi4/blob/d9093c75c71c2b33fbe86f32b25d138675ac22eb/psi4/src/psi4/libfock/points.cc#L767) line needs to be changed?

> For g2g, maybe just [this](https://github.com/psi4/psi4/blob/d9093c75c71c2b33fbe86f32b25d138675ac22eb/psi4/src/psi4/libfock/points.cc#L767) line needs to be changed? Confirmed to work locally. This fixes *a lot* of the failed tests above 🎉

> Those ints will be rate limiting, so it's important that they are threaded somehow. Yes -- and I assume they will be evaluated for *a lot* of points, which...

I think for your purpose, I'd use `ElectrostaticInt` from C++ (I'm using `ao_multipole_potential` because of tests from Python). `ElectrostaticInt` uses L2, whereas multipole potential ints will use M-D in the...

Need to put this on hold until https://github.com/evaleev/libint/issues/242 is clarified 😬

P.S.: Just as "proof" a screenshot of the WolframAlpha result: 

> Since you suspect an error in the Boys function, https://github.com/evaleev/libint/pull/243 may fix this. Only if `FmEval_Taylor` is used. I'm not sure but I think the default is `FmEval_Chebyshev7`...

I'm quite certain now that `erf_coulomb_gm_eval`/`erfc_coulomb_gm_eval` are correct, otherwise, the `erf_coulomb`/`erfc_coulomb` two-electron integrals would be wrong, too. Any hints on how to get to the bottom of this, @evaleev?

Maybe [this issue](https://github.com/maxscheurer/pycontact/issues/55) helps?