When finding nuclear forces for the true electronic wavefunction, thanks to the Hellmann-Feynman theorem we only need to consider explicit derivatives with respect to the nuclear coordinates $\mathbf{R}$: \begin{equation} \frac{dE}{d\mathbf{R}}=\left\langle\Psi\left|\frac{dH}{d\mathbf{R}}\right|\Psi\right\rangle+E\frac{d}{d\mathbf{R}}\left\langle\Psi\left|\right.\Psi\right\rangle=\left\langle\Psi\left|\frac{dH}{d\mathbf{R}}\right|\Psi\right\rangle \end{equation}

The 2nd term cancels as the overlap of the wavefunction with itself is constant (1 if normalized) and so the derivative is 0.

For approximate methods, this second term (the Pulay force or stress) does not generally go away, as the wavefunction is not an eigenfunction of the true Hamiltonian. For variational methods, these forces are said to go away in the limit of a complete basis. What is it about the complete basis limit that removes the Pulay forces for a variational method?

  • $\begingroup$ I am interested in this. Are there some references you could give which say that these forces go away in the CBS limit? I suppose that even in the CBS limit, these forces don't necessarily go away, if the wavefunction still for some reason doesn't have a norm of 1, so it will depend not only on the basis set but on the level of theory (coupled-cluster vs Hartree-Fock vs MP2 for example... perturbation theory possibly giving unphysical results at low order, and hence not giving perfectly physical wavefunctions even at the CBS limit). I would have to see the papers where they make this claim. $\endgroup$ Commented May 11, 2020 at 17:31
  • $\begingroup$ @NikeDattani This is one example that people seem to cite: journals.aps.org/prb/abstract/10.1103/PhysRevB.61.16207. The "alternate proof" on the Wikipedia page for the Hellmann-Feynman theorem suggests these terms should cancel at the CBS limit provided the method being used is variational, so coupled-cluster and MPn should still have Pulay forces even at the CBS limit. $\endgroup$
    – Tyberius
    Commented May 11, 2020 at 17:53
  • $\begingroup$ I still have to look at the Wikipedia page, but I would think that a terrible variational ansatz (like CISD for a sextuple bond like in the chromium dimer) would be a bad wavefunction even in the CBS limit: There is a reason why we have the denominator $\langle \psi|\psi\rangle$ in the variational formula, it's because $\langle \psi|\psi\rangle$ is not always 1 (especially in stochastic CI methods like FCIQMC). I will look into it when I have more time! Right now of course, the recruiting for more members is more important to me! $\endgroup$ Commented May 11, 2020 at 18:02
  • $\begingroup$ @NikeDattani it would definitely be a bad wavefunction. All this theorem is saying is that the forces would just come from the explicit dependence of the energy on that parameter, rather than having these Pulay terms. The calculated force could still be terrible, but it's the "correct" force for that method without having to compute any implicit dependence. $\endgroup$
    – Tyberius
    Commented May 11, 2020 at 18:16
  • 1
    $\begingroup$ I see that even if the wavefunction is not normalized, it can still be independent of R. I still haven't read the paper or Wikipedia page, but I see now that the theorem is probably saying that in the CBS limit, $\langle \psi | \psi \rangle$ won't have any type of dependence on R even if it's not normalized. Will look more into it when things calm down a bit! $\endgroup$ Commented May 11, 2020 at 18:21

2 Answers 2


I think a good way to see this is to simplify the stage a little. Imagine we want to solve the Schrödinger equation for a Hamiltonian $H$. For this we take a very simple basis, namely just two real-valued functions $\{f_1, f_2\}$. If we variationally optimise the trial wavefunction from this basis, we are presented with an approximation to the ground state in the linear combination $ \Psi = c_1 f_1 + c_2 f_2$ and some corresponding approximate energy $E$. Since we are variational, we have $$ 0 = \frac{dE}{dc_1} = 2 \left\langle \Psi \middle| H \frac{d\Psi}{dc_1} \right\rangle = 2 \left\langle \Psi \middle| H f_1 \right\rangle$$ and similarly $0 = \left\langle \Psi \middle| H f_2 \right\rangle$. Now the derivative of the energy wrt. positions can be written as \begin{equation} \frac{dE}{d\mathbf{R}}=\left\langle\Psi \middle|\frac{dH}{d\mathbf{R}}\Psi\right\rangle + 2 \left\langle \Psi \middle| H \frac{d\Psi}{d\mathbf{R}} \right\rangle \end{equation} with the second term being the Pulay forces of interest for us. Consider as an example the derivative wrt. $R_1$. Its Pulay term is $$ \left\langle \Psi \middle| H \left( c_1 \frac{df_1}{dR_1} + c_2 \frac{df_2}{dR_1} \right)\right\rangle. $$

When is this term zero? Either when both the derivative of $f_1$ and $f_2$ wrt. $R_1$ are zero, i.e. if the basis functions by itself is independent of atomic positions. This is the case for example for plane waves or generally all basis sets which are not atom-centred. The other option is if the derivatives $\frac{df_1}{dR_1}$ and $\frac{df_2}{dR_1}$ are itself basis functions or can be exactly represented by the basis. To see what happens then, assume this was the case. We could write $$ \frac{df_1}{dR_1} = k_{11} f_1 + k_{12} f_2 \quad\text{and}\quad \frac{df_2}{dR_1} = k_{21} f_1 + k_{22} f_2$$ for appropriate constants and obtain \begin{align} &\hspace{-30pt}\left\langle \Psi \middle| H \left[ (c_1k_{11} + c_2k_{21}) f_1 + (c_1k_{12} + c_2k_{22}) f_2 \right]\right\rangle \\ &= (c_1k_{11} + c_2k_{21}) \left\langle \Psi \middle| H f_1 \right\rangle + (c_1k_{12} + c_2k_{22}) \left\langle \Psi \middle| H f_2 \right\rangle \\ &= 0 \end{align} due to the first expressions we derived.

A complete basis set is just a special case of our argument where by definition all derivatives of basis functions can be represented by the basis itself, thus giving net zero Pulay forces.


To complement Michael's answer, Pulay forces can actually be dampened by a suitable design of the basis set, as noted by Pulay in 1977. We have recently shown in arXiv:2207.03587 that accurate Hellmann-Feynman forces can be evaluated with density functional theory using Gaussian basis sets that are not much bigger than established state-of-the-art triple-zeta basis sets.

  • $\begingroup$ This is interesting! It may be worthwhile to see if these basis sets can also reduce the errors of the Hellmann-Feynman non-adiabatic coupling matrix elements (see pubs.acs.org/doi/10.1021/acs.accounts.1c00312?ref=PDF for my recent benchmark results on this), especially the excited state-excited state NACMEs for which the Pulay force term involves CP-SCF equations and is a real computational burden. NACME calculations generally require less precision than gradients, which may be an advantage. $\endgroup$
    – wzkchem5
    Commented Oct 15, 2022 at 17:32

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .