Hellmann-Feynman Forces with Hartree-Fock

Question

The energy in the Hartree-Fock approximation is given as:

$$E_{HF}=\left<\psi_{HF} \left| \hat{H} \right| \psi_{HF} \right>=\sum_{i,j}P_{i,j}H_{i,j}^{core}+\frac{1}{2}\sum_{i,j,k,l}P_{i,j}P_{l,k}(ij||kl)+V_{NN} \tag{1}$$

The geometrical derivative of the Hartree-Fock energy can be shown to be [1]:

$$\frac{\partial E_{HF}}{\partial X_A}=\sum_{i,j}P_{i,j}\frac{\partial H_{i,j}^{core}}{\partial X_A}+\frac{1}{2}\sum_{i,j,k,l}P_{i,j}P_{l,k}\frac{\partial (ij||kl)}{\partial X_A}-\sum_{i,j}Q_{i,j}\frac{\partial S_{i,j}}{\partial X_A}+\frac{\partial V_{NN}}{\partial X_A} \tag{2}$$

The Hellmann-Feynman Theorem states:

$$\frac{dE}{d X_A}=\left<\psi \left| \frac{d\hat{H}}{d X_A} \right| \psi \right>$$

The Hellmann-Feynman Theorem implies that the energy derivative only depends on the parts of the Hamiltonian that have a dependency on the derivative. This leads to the geometrical derivative being equal to:

$$\frac{dE}{d X_A}=\left<\psi \left| \frac{d\hat{V}_{NN}}{d X_A} + \frac{d\hat{V}_{Ne}}{d X_A} \right| \psi \right> = \left<\psi \left| \frac{d\hat{V}_{NN}}{d X_A} \right| \psi \right> + \left<\psi \left| \frac{d\hat{V}_{Ne}}{d X_A} \right| \psi \right>=\left<\psi \left| Z_A\sum_B\frac{Z_B(X_B-X_A)}{R^3_{AB}} \right| \psi \right> + \left<\psi \left| -Z_A\sum_i\frac{X_i-X_A}{r^3_{iA}} \right| \psi \right>$$

Now the Hellmann-Feynman theorem can be applied to find the derivative geometrical of the Hartree-Fock energy:

$$\frac{\partial E_{HF,\ce{Hellmann-Feynman}}}{\partial X_A}=\left<\psi_{HF} \left| \frac{\partial\hat{H}}{\partial X_A} \right| \psi_{HF} \right>=\sum_{i,j}P_{i,j}\frac{\partial V_{Ne,ij}}{\partial X_A}+\frac{\partial V_{NN}}{\partial X_A} \tag{3}$$

I tried to implement equation $(2)$ and equation $(3)$, and used them to calculate, the force between H and Li for different distances, with cc-pVDZ and cc-pVTZ. and got the following:

Hartree-Fock refers to equation $(2)$ and Hellmann-Feynman refers to equation $(3)$. I have read that the Hellmann-Feynman theorem is only valid for exact solutions, but should become a better and better approximation, when going towards the Hartree-Fock limit. As can be seen in the picture, equation $(3)$ performs worse with the larger basisset cc-pVTZ. This now leads me to my question, were am I wrong in my understanding of the application of the Hellmann-Feynman theorem to the Hartree-Fock approximation?

[1] Attila Szabo, Neil S. Ostlund; Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; equation C. 12

Answer 1

It is not at all surprising that the forces you get from application of the Hellmann-Feynman theorem do not match the true force. The reason is the appearance of Pulay forces: the atomic orbital basis set depends on the nuclear coordinates, and changing the position of the nuclei also changes the basis set. This yields the contribution you're missing in the plots.

cc-pVDZ and cc-pVTZ are still pretty small basis sets. I would expect that if you add diffuse functions, the result improves quite a bit, since the basis set covers more degrees of freedom. Furthermore, to get something meaningful you'd have to compare e.g. cc-pV5Z vs cc-pV6Z, which are closer to the complete basis set limit.

However, already in your plots you can see that the agreement is better from cc-pVDZ to cc-pVTZ when you look at the far-range behavior. A better way to plot these would be to plot the difference, i.e. the Pulay force.

I would also imagine that decontracting the basis sets helps, because if you imagine just moving the electron density but not the basis functions, then you will see that the description of the core orbitals is impoverished quite a bit as the contraction coefficients are wrong.

Actually, the result in cc-pV6Z is probably still not very close to the Hellman-Feynman limit. When you move the nucleus, 1s core orbitals start to gather p, f, etc character. Thus, when you move the nuclei but not the basis functions, your original basis set should have e.g. very tight p functions that are able to describe this effect.

If you want to toy around with this, you should look at calculations in larger and larger primitive basis sets. Even-tempered exponents should work fine. Start with a small set of s functions, and study how the result changes when you add steeper or more diffuse functions. Next, look at the p shell: which exponent changes the result the most? And then again start adding steeper and more diffuse functions. Repeat ad nauseam. (This is called completeness-optimization, and I've published an automatic procedure for this.)