In situoptimization of a set of localized orbitals with respect to a systematically improvable basis set independent of the position of the atoms, such as psinc functions, would theoretically eliminate the correction due to Pulay forces from the total ionic forces. We demonstrate that for strict localization constraints, especially with small localization regions, there can be non-negligible Pulay forces that must be calculated as a correction to the Hellmann-Feynman forces in the ground state. Geometry optimization calculations, which rely heavily upon accurate evaluation of the total ionic forces, show much better convergence when Pulay forces are included. The more conventional case, where the local orbitals remain fixed to pseudo-atomic orbital multiple-ζ basis sets, also benefits from this implementation. We have validated the method on several test cases, including a DNA fragment with 1045 atoms
