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 HellmannFeynman 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 HIPPOlike 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.
I’m just getting back to this. That tutorial is very helpful, thanks!
I see three functions, potential_grad
, overlap_grad
, and kinetic_grad
, that look like they’re doing what I want. I realize this is probably pretty elementary, but can you spell out what each of these is doing? If I want just the HellmannFeynman contribution to the forces, which do I need?
potential_grad
 The oneelectron gradient integral associated with the electronnuclear potential
overlap_grad
 The oneelectron gradient integral associated with how the perturbation changes the overlap between the LCAOMOs, because the perturbation will change the AOs. Even HartreeFock theory needs this term.
kinetic_grad
 The oneelectron gradient integral associated with the electronic kinetic energy
I can’t say “which of these you need” because I still don’t quite understand your big picture problem. (I have less experience with things like this than @andysim, so bear with me.) For computations done in a finite basis set, the HellmannFeynman derivative of the Hamiltonian does produce overlap and twoelectron terms. As you move the atoms, you change the basis your AOs are in, so you’ve changed the operators in the Hamiltonian themselves. The overlap terms are then necessary to make sure the MOs you’ve defined your Hamiltonian in are still orthonormal, and the twoelectron terms are necessary because even though there are no nuclei involved in the electronrepulsion interaction, you have changed the twoelectron operator itself by changing the AOs. The Wikipedia example of using HellmannFeynmann for molecular forces does not include these terms arising from a finite basis set.
Thanks for the help, @jmisiewicz!
The way I think about it, there are 3 contributions to the force.

The HellmannFeynman contribution. This is what is spelled out in the Wikipedia article you linked.

The finite basis set contribution. This would go away if you were using sufficiently complete basis set, but in reality you need to account for the fact that the basis functions move slightly when you move the nuclei.

The nonvariational contribution. This come from the fact that HF theorem assumes that your wavefunction is a minimum with respect to your Hamiltonian. If it is not, then you have expensive chain rule terms to compute.
What I’m trying to do is figure out how big each of these contributions are. The big picture here is this: I have a molecular electron density represented as a linear combination of auxiliary basis functions. I want to calculate contributions 1 and 2.
Let me add some more detail on my use case since it looks like that might be helpful.
I have a machine learning algorithm that I am training to reproduce very accurate electron densities. The output of this ML model is a molecular density represented in an auxiliary basis set: density = sum(ci * Xi), where Xi are the basis functions of the auxiliary basis set and ci are the learned coefficients. I want to be able to read in this density to psi4 and compute forces from it.
Specifically, I want to start with the simple HellmannFeymann forces as defined in the linked Wikipedia article. I think this just looks like using computing the the OEI electronnuclear gradient integrals, which psi4 already has machinery to to.
I understand that to the extent that my training data is nonvariational and my basis set is noncomplete, there will be an error. I would like to understand how big those errors are by computing each of those terms independently.
A secondary, but related goal is independent of the ML task. If I do a CCSD(T) calculation with a very large basis set, I would like to know how big the contributions to the force are from the 3 classes I listed above. My hunch is that even though CCSD(T) is nonvariational, it’s density is so close to the FCI density, that the contributions from 2 and 3 will be small. I think it would be great to test this hypothesis.
It looks like for the auxiliary basis set density problem, what I actually need is ao_oei_deriv1
. The documentation looks like this:
ao_oei_deriv1
( self: psi4.core.MintsHelper , arg0: str , arg1: int ) → List[psi4.core.Matrix]
What are arg0 and arg1 here? I thought arg1 is the atom index, but putting in an integer greater than the number of atoms still gives a result.
It also appears that this method gives a result per atom that is 3 matrices, presumably x, y and z. But each of these matrices in num_basis_funcXnum_basis_func, which I did not expect. I was expecting it to return just a single num_basis_func vector (ie. the nuclearao gradient for each ao)?