A new method for the multipole evaluation of contracted Cartesian Gaussian-based electron repulsion integrals is described, and implemented in linear scaling methods for computation of the Hartree–Fock exchange matrix. The new method, which relies on a nonempirical multipole acceptability criterion [J. Chem. Phys. 109, 8764 (1998)], renders the work associated with integral evaluation independent of the basis set contraction length. Benchmark calculations on a series of three-dimensional water molecule clusters and graphitic sheets with highly contracted basis sets indicate that the new method is up to 4.6 times faster than a well optimized direct integral evaluation routine. For calculations involving lower levels of contraction a factor of 2 speedup is typically observed. Importantly, the method achieves these large gains in computational efficiency while maintaining numerical equivalence with standard direct self consistent field theory.