Our interest in polarization integrals stems from the direct reaction field method.[3] DRF divides the system into an active region (solute, QM), where electrons are treated using quantum mechanics, and an inert region (solvent, MM) of polarizable atoms that model the polarization response of the solvent to the solute density. Assuming that the basis set for the QM part contains N AO atom-centered orbitals and is surrounded by N pol polarizable sites, the number of unique polarization integrals scales as O(N 2 AO N pol ). Although the number of combinations of basis functions increases quadratically, asymptotically the number of non-zero one-electron matrix elements scales roughly linearly with the size of the system. This is because atomic orbitals decay exponentially, so that two basis functions have to be relatively close for their product to give a non-negligible contribution to the charge density. Therefore, local operators, like Eq. ( ), are sparse in the atomic orbital basis. To take advantage of this sparsity, we seek a simple upper bound to the integrals that allows us to efficiently estimate their magnitude and prescreen out negligible integrals.