Derivative integrals

[Emieeel]

To estimate Hellman-Feynman forces (with the SAOOVQE for example) we use

$$ \frac{\partial E}{\partial R} = \langle \Psi | \frac{\partial H}{\partial R} | \Psi \rangle $$

Suppose we have $N$ atoms and $M$ orbitals, meaning $R$ is a vector of size 3N. To calculate $\frac{\partial H}{\partial R}$, we need integral derivatives of four things:

• the overlap matrix $\langle \phi_p| \partial_R \phi_q \rangle$ (array of M x M x N x 3)
• the nuclear repulsion $\frac{\partial E_{\rm nuc}}{\partial R}$ (array of N x 3)
• the one-body integrals $\frac{\partial h_{pq}}{\partial R}$ (array of M x M x N x 3)
• the two-body integrals $\frac{\partial g_{pqrs}}{\partial R}$ (array of M x M x M x M x N x 3)

It would be great if these could be incorporated as possible data in our hdf5 scheme.

I already have code to generate these for both Psi4 and PySCF in AO basis (the latter being quite a hassle, see this issue). This is still quite inefficient so any improvements are welcome. I can make a branch here where I can incorporate them with the help of @Cmurilochem.

[Cmurilochem]

Thanks @Emieeel. I suspected that this idea would be handy in replacing the need to deal with the "orbital optimization" problem by making the derivatives of the hamiltonian directly instead of MO coeffs. Please, correct me if I am wrong.

Although I think that this could be nice to implement, specially if this is related to #24, I still think that for now and with the current infrastructure of the program, it would be hard to extend the QCSchema, specially because (for validation reasons) we are currently borrowing the implemented QCSchema dataclass (see here and here) implemented in qiskit-nature.

But if you have a nice idea on how to do this in connection with #24, please, feel free to make a branch from the main and try them out.