  The Hellmann-Feynman (HF) theorem provides a way to compute forces directly
from the electron density, enabling efficient force calculations for large
systems through machine learning (ML) models for the electron density. The main
issue holding back the general acceptance of the HF approach for atom-centered
basis sets is the well-known Pulay force which, if naively discarded, typically
constitutes an error upwards of 10 eV/Ang in forces. In this work, we
demonstrate that if a suitably augmented Gaussian basis set is used for density
functional calculations, the Pulay force can be suppressed and HF forces can be
computed as accurately as analytical forces with state-of-the-art basis sets,
allowing geometry optimization and molecular dynamics to be reliably performed
with HF forces. Our results pave a clear path forwards for the accurate and
efficient simulation of large systems using ML densities and the HF theorem.
