How does a complete basis set remove Pulay forces/stress?

Question

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?

Answer 1

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.

Answer 2

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.