We propose a new strategy to solve the Many-Body Dispersion (MBD) model by Tkatchenko, DiStasio Jr. and Ambrosetti. Our approach overcomes the original O(N**3) computational complexity that limits its applicability to large molecular systems within the context of O(N) Density Functional Theory (DFT). First, in order to generate the required frequency-dependent screened polarizabilities, we introduce an efficient solution to the Dyson-like self-consistent screening equations. The scheme reduces the number of variables and, coupled to a DIIS extrapolation, exhibits linear-scaling performances. Second, we apply a stochastic Lanczos trace estimator resolution to the equations evaluating the many-body interaction energy of coupled quantum harmonic oscillators. While scaling linearly, it also enables communication-free pleasingly-parallel implementations. As the resulting O(N) stochastic massively parallel MBD approach is found to exhibit minimal memory requirements, it opens up the possibility of computing accurate many-body van der Waals interactions of millions-atoms’ complex materials and solvated biosystems with computational times in the range of minutes.