Only point densities are sampled when matrix elements of a delta function operator are evaluated. Thus, the Fermi contact contribution to nuclear spin—spin coupling tensors is very sensitive to local inaccuracies of the wavefunctions. An alternative expression for the delta function which uses a global representation of δ( r) improves the agreement between calculated and experimental coupling constants, especially for small basis sets. The method is illustrated for both small and large basis sets for the HD molecule in a second-order polarization propagator approach.