Pulay stress

The Pulay stress or Pulay forces (named for Peter Pulay) is an error that occurs in the stress tensor (or Jacobian matrix) obtained from self-consistent field calculations (Hartree–Fock or density functional theory) due to the incompleteness of the basis set.[1][2][3]

A plane wave basis set is created for the hexagonal lattice (left), using the reciprocal lattice vectors inside the red circle. Then the lattice relaxes into a cubic symmetry (right). Keeping the red circle basis constant results in lattice vectors taken from an ellipsoid instead of a spherical area (compare to the blue circle).

A plane-wave density functional calculation on a crystal with specified lattice vectors will typically include in the basis set all plane waves with energies below the specified energy cutoff. This corresponds to all points on the reciprocal lattice that lie within a sphere whose radius is related to the energy cutoff. Consider what happens when the lattice vectors are varied, resulting in a change in the reciprocal lattice vectors. The points on the reciprocal lattice which represent the basis set will no longer correspond to a sphere, but an ellipsoid. This change in the basis set will result in errors in the calculated ground state energy change.

The Pulay stress is often nearly isotropic, and tends to result in an underestimate of the equilibrium volume.[2] Pulay stress can be reduced by increasing the energy cutoff. Another way to mitigate the effect of Pulay stress on the equilibrium cell shape is to calculate the energy at different lattice vectors with a fixed energy cutoff.[2]

Similarly, the error occurs in any calculation where the basis set explicitly depends on the position of atomic nuclei (which are to change during the geometry optimization). In this case, the Hellmann–Feynman theorem – which is used to avoid derivation of many-parameter wave function (expanded in a basis set) – is only valid for the complete basis set.[3] Otherwise, the terms in theorem's expression containing derivatives of the wavefunction persist, giving rise to additional forces – the Pulay forces:[4]

${\displaystyle {\begin{aligned}{\frac {\mathrm {d} E}{\mathrm {d} \mathbf {R} }}&={\frac {\mathrm {d} }{\mathrm {d} \mathbf {R} }}\langle \psi |{\hat {H}}|\psi \rangle \\&={\bigg \langle }{\frac {\mathrm {d} \psi }{\mathrm {d} \mathbf {R} }}{\bigg |}{\hat {H}}{\bigg |}\psi {\bigg \rangle }+{\bigg \langle }\psi {\bigg |}{\hat {H}}{\bigg |}{\frac {\mathrm {d} \psi }{\mathrm {d} \mathbf {R} }}{\bigg \rangle }+{\bigg \langle }\psi {\bigg |}{\frac {\mathrm {d} {\hat {H}}}{\mathrm {d} \mathbf {R} }}{\bigg |}\psi {\bigg \rangle }\\&=\underbrace {E{\bigg \langle }{\frac {\mathrm {d} \psi }{\mathrm {d} \mathbf {R} }}{\bigg |}\psi {\bigg \rangle }+E{\bigg \langle }\psi {\bigg |}{\frac {\mathrm {d} \psi }{\mathrm {d} \mathbf {R} }}{\bigg \rangle }} _{\ 0\ for\ complete\ basis\ set}+{\bigg \langle }\psi {\bigg |}{\frac {\mathrm {d} {\hat {H}}}{\mathrm {d} \mathbf {R} }}{\bigg |}\psi {\bigg \rangle }.\end{aligned}}}$

The presence of Pulay forces makes the optimized geometry parameters converge slower with increasing basis set.[3] The way to eliminate the erroneous forces is to use nuclear-position-independent basis functions,[4] to explicitly calculate and then subtract them from the conventionally obtained forces, or to self-consistently optimize the center of localization of the orbitals.[3]