The Hellmann-Feynman (HF) theorem provides a way to compute forces directly from the electron density, affording an approach to calculating forces of large systems with machine learning (ML) models that predict electron density. The primary 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 construct specialized atom-centered Gaussian basis sets to reduce the Pulay force, and demonstrate the basis sets' effectiveness in computing accurate HF forces. We find that HF forces computed using the σNZHF (N = Single, Double, Triple) basis sets developed in this work yield comparable accuracy to forces computed with the Pulay term using size matched cc-pVNZ [1] and pcseg-N [2] basis sets for water clusters and pcseg-N and aug-pcseg-N basis sets for DNA fragments. Our results illustrate that the σNZHF basis sets yield HF forces with state-of-the-art accuracy, paving a clear path forwards for accurate and efficient calculations of forces for large systems using the HF theorem and ML densities.