  Pulay terms arise in the Hellman-Feynman forces in electronic structure
calculations when one employs a basis set made of localized orbitals that move
with their host atoms. If the total energy of the system depends on a subspace
population defined in terms of the localized orbitals across multiple atoms,
then unconventional Pulay terms will emerge due to the variation of the orbital
nonorthogonality with ionic translation. Here, we derive the required exact
expressions for such terms, which cannot be eliminated by orbital
orthonormalization. We have implemented these corrected ionic forces within the
linear-scaling density functional theory (DFT) package ONETEP, and have used
constrained DFT to calculate the reorganization energy of a pentacene molecule
adsorbed on a graphene flake. The calculations are performed by including
ensemble DFT, corrections for periodic boundary conditions, and empirical Van
der Waals interactions. For this system we find that tensorially invariant
population analysis yields an adsorbate subspace population that is very close
to integer-valued when based upon nonorthogonal Wannier functions, and also but
less precisely when using pseudoatomic functions. Thus, orbitals can provide a
very effective population analysis for constrained DFT. Our calculations show
that the reorganization energy of the adsorbed pentacene is typically lower
than that of pentacene in the gas phase. We attribute this effect to steric
hindrance.
