Is there any straightforward way to compute atomic forces from a density? Specifically, I am looking for a way to take a density represented in an auxiliary, density fitting basis set and compute Hellmann-Feynman forces on each atom.

# Atomic Forces from Electron Density

It’s easy to get the derivative integrals in any basis, which can then be dotted with the density to get the gradients. The problem is figuring out how to handle the overlap and two electron terms. In SCF the two electrons don’t explicitly see each other so the two particle density matrix is easily expressed as products of the one particle densities; these can then be contracted with the two electron integral derivatives. These TEI derivatives are also available through MintsHelper. The TLDR answer is that the machinery all exists in terms of integral helpers, it’s just a case of thinking about what the derivative expression should look like, which depends on the nature of the density you’re handling (frozen finite width multipoles, I’m guessing?).

Maybe I’m thinking about this too simplistically, but I don’t think I need any overlap terms. I’m representing the density in a single auxiliary basis set, so I think all I need is Fi = integral(Zi * density / r^3). Am I missing something fundamental here?

Sounds very reasonable and we have the integrals that are needed. I wasn’t sure if the question was about just a standalone HIPPO-like entity or one in the presence of a QM molecule. The latter situation is doable, but harder. Just handling a standalone entity should be very easy. This Psi4Numpy tutorial gives a complete working example of how it’s done for SCF, and I suspect your problem will turn out to be simpler because of the lack of reorthonormalization gradients (the overlap terms) and the lack of two electron terms, if I’m understanding the problem correctly.