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.
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)?
We are quickly getting to the end of my ability to help, because I don’t know much about electron density based theories. Gradients with reduced density matrices, I can work with. It sounded like @andysim may be more knowledgeable here?
This is a case where it’s faster to look at the code yourself.
string: What type of derivatives do you want?
int: Which atom number? Because C++ is C++, array indices loop around. That’s why atom numbers grater than the number of atoms don’t (obviously) break anything.
As for return types, that’s because each oneelectron integral in the Hamiltonian is associated with two AOs, so of course the derivative of an integrals will be associated with two AOs.
And yes, the three matrices are indeed X, Y, and Z.
Thanks for the explanation. It turns out that this is not the function I want, so sorry for the runaround. All that I’m really looking for is the overlap integral derivative of each individual AO on a nuclear center. I know this is not directly an operation that shows up in mainstream QM methods, but I was hoping I might be able to highjack something.
If nothing else, thanks for helping me understand the code, @jmisiewicz. I will leave a concise statement of the problem below in case anyone happens to come along with some new insight into my problem.
Problem Statement:
I want to compute the following integral to compute the force on a nucleus, i, due to an electron density represented in an auxiliary atomic basis:
F_i = sum_j(integral(Z_i * c_j * X_j / R dv) )
Here, the sum is over the j basis functions, X, and coefficients, c, and the integral is over all space.
There are two problems to be solved here:

Read in the auxiliary basis and coefficients. This is referenced in this post.

Compute the integral above.
The equation looks like you want the derivative of a single GTO w.r.t to nuclear displacement?
Is this the equation you want to solve?
latex: F_i = \sum_j(\int (Z_i \cdot c_j \cdot \frac{\partial X_i}{\partial R_i} dv)
(embedding html code from https://www.codecogs.com/latex/eqneditor.php?lang=enus)
The density you define [=sum(c_i * X_i)] looks like an LCAO expansion for a single MO orbital?
Hi @hokru,
Yes, everything you say is correct (way to go with the embedded LaTeX!), with the exception that it should be X_j. The density does look like a LCAO expansion of a single MO.
I thought I’d post here to say I did end up finding a (much simpler) way to compute HellmannFeynman Forces from psi4. I just compute the electric field at the nuclei with something like,
e, wfn = gradient('scf',return_wfn=True)
oeprop(wfn, 'GRID_FIELD')
And then multiply the resulting electric field by the nuclear charge. This works and gives the correct HF forces on each atom.
What I still do not know how to do is give my own density matrix to oeprop
. That would be very cool.
oeprop
reads the density matrix from the wavefunction, so you just need to set the density matrix of the wavefunction object.