When applied to electronic wavefunctions calculated with Gaussian-type basis functions, the Hiller-Sucher- Feinberg (HSF) identity improves the accuracy of the electron density at non-hydrogen nuclei by more than an order of magnitude, yielding approximate electron nuclear cusps. However, the HSF electron densities at hydrogen nuclei bound to heavy atoms are greatly overestimated. This phenomenon is associated with the asymptotic behaviour of the HSF density, which incorrectly decreases to a constant when the sum of Hellmann-Feynman forces acting on nuclei is finite. A method for constraining variational wavefunctions to yield vanishing Hellmann-Feynman forces is described. Hartree-Fock calculations of the constrained HSF (CHSF) electron densities with the 6-31G, 6-31G**, and 6-311++G** basis sets are reported at the nuclei of various diatomic molecules, and are compared with their corresponding conventional, HSF, and Hartree-Fock limit values. These calculations show that differences between HSF and CHSF densities are minor at non-hydrogen nuclei. Importantly, the calculated HF/6-311 + +G** CHSF densities are on average three times more accurate than the conventional densities at hydrogen nuclei.