For testing, all integrals in Eqn. (109) were evaluated numerically. The centers of the two basis functions were drawn randomly from a ball of radius 5 Å, typical of the size of QM regions in QM/MM calculations. The polarizable atom was always placed at the origin. The exponents of the basis functions were also drawn randomly from the interval [0.01, 2.01]. In total, parameters for 50 pairs of shells are generated. For a pair of centers, all combinations of angular momenta up to d-functions were enumerated. Figure shows a scatter plot of the upper bounds versus the exact polarization integrals. The diagonal straight black line indicates a hypothetical perfect bound. All points in Figure lie above the diagonal black line, so the bound is rigorous. Furthermore, for most integrals, the bound is within seven orders of magnitude, which is typical for CS bounds involving Cartesian Gaussian basis functions. If the threshold for neglecting one-electron integrals is set to 10 -8 Hartree (red horizontal line in Fig. ), the CS bound shows that most integrals can be neglected. We note, however, that for some small integrals the bound exceeds the true integral by more than seven orders of magnitude, shown within the gold circle in Fig. . These can be understood as arising from symmetry. For some combinations of angular components, the integrals are numerically zero by symmetry, but they are still included since the screening is applied to blocks of basis functions with the same total angular momentum. The number of these integrals is relatively small, however, so including them does not incur a large computational effort.