Using a previously introduced weight scheme, microbatching, and grid compression [R. E. Stratmann, G. E. Scuseria and M. J. Frisch, Chem. Phys. Lett. 257, 213 (1996)], we significantly speed up the numerical integration of the exchange-correlation contribution to the Coupled-Perturbed Kohn–Sham equations. In addition, we find that the nature of the integrand is such that it is possible to employ substantially fewer grid points in the quadrature and to use the Gaussian very Fast Multipole Method (GvFMM) with very short multipole expansions for the Coulomb contribution, with negligible loss in accuracy. As a representative example, the computational demand for the exchange-correlation portion of a coronene (C24H12) frequency calculation with a 3-21G basis is reduced by more than one order of magnitude. The overall speed up achieved in this calculation is between a factor of 4 to 6, depending on the specific functional. We also present sample calculations using polarized bases, gradient-corrected functionals, and on even larger systems (C54H18 and C96H24), to illustrate the various effects and improvements that we have accomplished.