In QM/MM studies with large MM regions, the calculation of electrostatic embedding integrals can become a computational bottleneck. To overcome this problem, an asymptotic expansion for nuclear attraction-type integrals is developed. As a result, the long-range interactions between the QM and MM atoms reduce to atom-centered multipole moment-like expansions. The algorithm uses a natural spatial division of the molecular structure. To further improve the computational performance, a cutoff radius for the multipole moment-like expansion is introduced. The new code was validated and benchmarked with deMon2k/CHARMM QM/MM calculations on an RNA polymerase II model with almost 350 000 atoms. It is shown that the computational time for the calculation of the embedding integrals in this system can be reduced below 200 s on a small parallel architecture (eight cores) without a loss of accuracy.