Hellmann–Feynman forces & the Pulay term

The bridge from electronic structure to the lattice: how the force on a nucleus is computed once the electronic problem is solved, and the one correction — the Pulay term — that decides whether those forces, and the phonons built from them, are right. The force is the gradient of the Born–Oppenheimer PES, so everything downstream — geometry optimization, molecular dynamics, and the force constants of finite-displacement phonons — rests on getting this gradient exactly.

The Hellmann–Feynman theorem

Let $|{\Psi }_{\lambda }⟩$ be an eigenstate of a Hamiltonian ${\hat{H}}_{\lambda }$ depending on a parameter $\lambda$ (here a nuclear coordinate $R_A$), with energy $E_{\lambda }=⟨{\Psi }_{\lambda }|{\hat{H}}_{\lambda }|{\Psi }_{\lambda }⟩$ and $⟨{\Psi }_{\lambda }|{\Psi }_{\lambda }⟩=1$. Differentiating,

$$\frac{dE}{d\lambda }=\underset{\text{wavefunction-derivative terms}}{\underbrace{⟨\frac{\partial \Psi }{\partial \lambda }\left|\hat{H}\right|\Psi ⟩+⟨\Psi \left|\hat{H}\right|\frac{\partial \Psi }{\partial \lambda }⟩}}+⟨\Psi \left|\frac{\partial \hat{H}}{\partial \lambda }\right|\Psi ⟩.$$

For an exact eigenstate the first two terms collapse: $⟨{\partial }_{\lambda }\Psi |\hat{H}|\Psi ⟩+c.c.=E{\partial }_{\lambda }⟨\Psi |\Psi ⟩=E{\partial }_{\lambda }(1)=0$. More generally they vanish whenever the energy is stationary with respect to all variational parameters of $|\Psi ⟩$ — the key point for approximate methods. What survives is the Hellmann–Feynman theorem

$$\boxed{\frac{dE}{d\lambda }=⟨\Psi \left|\frac{\partial \hat{H}}{\partial \lambda }\right|\Psi ⟩}$$

so the force on nucleus $A$ is the expectation value of an operator, no derivative of the wavefunction required,

$${\mathbf{F}}_A=-\frac{\partial E}{\partial {\mathbf{R}}_A}=-⟨\Psi \left|\frac{\partial \hat{H}}{\partial {\mathbf{R}}_A}\right|\Psi ⟩.$$

The statement is exact at a variational stationary point and is the reason forces are cheap: solve the electronic problem once, then evaluate a single expectation value rather than finite-differencing the energy or propagating ${\partial }_{\lambda }\Psi$.

The electrostatic (Feynman) force theorem

For the electronic Hamiltonian, the only ${\mathbf{R}}_A$-dependence is in the electron–nucleus and nucleus–nucleus Coulomb terms — the kinetic and electron–electron operators are independent of nuclear positions. Their derivative is purely electrostatic, so the Hellmann–Feynman force evaluates to a classical expression, the electrostatic force theorem (Feynman; later Politzer & Murray):

$${\mathbf{F}}_A=Z_A\left[\int n(\mathbf{r})\frac{\mathbf{r}-{\mathbf{R}}_A}{|\mathbf{r}-{\mathbf{R}}_A{|}^3}d\mathbf{r}-\sum _{B\ne A}{Z}_B\frac{{\mathbf{R}}_B-{\mathbf{R}}_A}{|{\mathbf{R}}_B-{\mathbf{R}}_A{|}^3}\right].$$

The exact force on a nucleus is literally the classical Coulomb attraction of the ground-state electron density $n(\mathbf{r})$ plus the bare repulsion of the other nuclei — no kinetic, exchange, or correlation piece appears explicitly. All quantum mechanics enters only through the shape of $n(\mathbf{r})$. The bonding force, chemically, is the electron density piled up in the bond pulling the nuclei together. This is conceptually clean but numerically delicate: the force is dominated by the density very near the nucleus, so it amplifies small density errors there.

The Pulay (wavefunction) force — the practical core

The Hellmann–Feynman theorem holds only when the energy is fully stationary with respect to every parameter in the wavefunction. With atom-centred basis sets that condition is silently broken, and this is the single most important caveat in computing forces.

Expand the orbitals in basis functions ${\chi }_{\mu }(\mathbf{r}-{\mathbf{R}}_{A(\mu )})$ that are attached to the nuclei — Gaussian-type orbitals in CRYSTAL, numerical atomic orbitals in SIESTA. When nucleus $A$ moves, the basis functions centred on it move too, so $|\Psi ⟩$ acquires an explicit ${\mathbf{R}}_A$-dependence through the basis, not only through the variationally optimized LCAO coefficients. The total derivative then carries an extra term,

$$\frac{dE}{d{\mathbf{R}}_A}=\underset{\text{Hellmann–Feynman}}{\underbrace{⟨\Psi \left|\frac{\partial \hat{H}}{\partial {\mathbf{R}}_A}\right|\Psi ⟩}}+\underset{\text{Pulay (wavefunction) force}}{\underbrace{\sum _{\mu }\frac{\partial E}{\partial {\chi }_{\mu }}\frac{\partial {\chi }_{\mu }}{\partial {\mathbf{R}}_A}}},$$

the Pulay force (a.k.a. incomplete-basis-set correction, or wavefunction force), arising because the variational condition $\partial E/\partial {\chi }_{\mu }=0$ fails for a finite basis — the energy is stationary with respect to the coefficients but not with respect to the basis-function positions. In practice it is evaluated through derivatives of the overlap and Hamiltonian integrals (an “energy-weighted density matrix” times $\partial \mathbf{S}/\partial {\mathbf{R}}_A$ term).

Neglecting the Pulay term gives wrong forces, and the damage propagates. Wrong forces mean wrong gradients of the PES; finite-differencing those to build force constants then corrupts the dynamical matrix and the phonon spectrum — the failure shows up as spurious imaginary modes or shifted frequencies that are artifacts of the missing term, not physics. Any atom-centred-basis phonon workflow must include Pulay forces.

When the Pulay term vanishes

Two clean cases recover the exact Hellmann–Feynman force:

• Position-independent bases. Plane waves are not attached to nuclei — $\partial {\chi }_{\mathbf{G}}/\partial {\mathbf{R}}_A=0$ — so the Pulay force is identically zero and the bare Hellmann–Feynman expression is exact. (Pseudopotential and PAW formalisms add their own non-local-projector contributions, but those are part of $\partial \hat{H}/\partial {\mathbf{R}}_A$, not Pulay terms.) This is a real advantage of plane-wave codes for forces.
• Complete bases. In the complete-basis limit the energy is stationary against any basis variation, so the Pulay term goes to zero with basis-set incompleteness. The Pulay force is thus a measure of how far the finite basis is from completeness.

The underlying distinction is variational vs non-variational: at a true variational minimum the Hellmann–Feynman theorem applies and the only extra force is the explicit basis-motion (Pulay) term; for genuinely non-variational energies (e.g. some correlated post-HF methods) even the Hellmann–Feynman term must be supplemented by a response (Z-vector / relaxed-density) contribution. For ground-state SCF in a position-independent basis, neither correction is needed and the electrostatic force theorem holds exactly.

Prerequisites

Born–Oppenheimer PES

Builds toward: finite-displacement phonons (forces from displaced supercells) · DFPT (the same force constants via analytic linear response) · MLIPs (forces as analytic-gradient training labels)

Key references

• The theorem — H. Hellmann, Einführung in die Quantenchemie (Deuticke, Leipzig, 1937); R. P. Feynman, Phys. Rev. 56, 340 (1939).
• The basis-set (Pulay) force — P. Pulay, Mol. Phys. 17, 197 (1969).
• The electrostatic force theorem in modern context — P. Politzer & J. S. Murray, Theor. Chem. Acc. 137, 14 (2018).