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:

enter image description here

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

  • $\begingroup$ Related reading: books.google.com/… $\endgroup$
    – Tyberius
    Commented Aug 17, 2017 at 0:39
  • 2
    $\begingroup$ You may wish to see the paper: Bakken, V., Helgaker, T., Klopper, W., & Ruud, K. (1999). The calculation of molecular geometrical properties in the Hellmann—Feynman approximation. Molecular Physics, 96(4), 653-671. This makes it clear that a complete basis is needed for the Hellmann-Feynman theorem to hold. Also, one thing which I don't quite understand, but seems to be necessary is that the AOs are not supposed to remain centered on the atoms as the geometry is distorted. Thus, you are always calculating in reference to the original geometry, and convergence is necessarily slower. $\endgroup$
    – jheindel
    Commented Aug 17, 2017 at 4:43
  • $\begingroup$ @jheindel "AOs are not supposed to remain centered on the atoms as the geometry is distorted" Without doing out the math to prove this, it would ensure that any derivative of an AO w.r.t. a nuclear coordinate is zero. $\endgroup$ Commented Aug 18, 2017 at 13:08
  • $\begingroup$ @pentavalentcarbon I agree that the concept is weird and I don't really understand it. I just brought it up because it's unusual so it would have been something easy to do wrong assuming the paper is right. I don't think it's necessarily zero though. All that having the AOs uncentered does is cause you to evaluate the derivative further away from the peak of the gaussian. I don't see why this should be zero. $\endgroup$
    – jheindel
    Commented Aug 18, 2017 at 15:38
  • $\begingroup$ @ErikKjellgren While in general you would expect improvement as the basis set is enlarged, this isn't necessarily a guarantee so I don't think it is too surprising to see slightly higher errors. $\endgroup$
    – Tyberius
    Commented Sep 25, 2018 at 14:50

1 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.)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.