Atomic Forces from Electron Density

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.

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.

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 Hellmann-Feynman contribution to the forces, which do I need?

potential_grad - The one-electron gradient integral associated with the electron-nuclear potential
overlap_grad - The one-electron gradient integral associated with how the perturbation changes the overlap between the LCAO-MOs, because the perturbation will change the AOs. Even Hartree-Fock theory needs this term.
kinetic_grad - The one-electron 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 Hellmann-Feynman derivative of the Hamiltonian does produce overlap and two-electron 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 two-electron terms are necessary because even though there are no nuclei involved in the electron-repulsion interaction, you have changed the two-electron operator itself by changing the AOs. The Wikipedia example of using Hellmann-Feynmann 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.

  1. The Hellmann-Feynman contribution. This is what is spelled out in the Wikipedia article you linked.

  2. 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.

  3. The non-variational 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 Hellmann-Feymann forces as defined in the linked Wikipedia article. I think this just looks like using computing the the OEI electron-nuclear gradient integrals, which psi4 already has machinery to to.

I understand that to the extent that my training data is non-variational and my basis set is non-complete, 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 non-variational, 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 nuclear-ao gradient for each ao)?