The linear-scaling of the presented algorithm, the low memory requirements, its natural parallel implementation, the beneficial self-averaging of errors together with the possibility of employing a limited set of ab initio derived parameters opens up the possibility of evaluating many-body dispersion energies of large biomolecular systems within a few minutes'time. This is exemplified below for the solvated [Ala] 40 test system counting approximately 100000 atoms shown in Figure . Firstly atom types are assigned to each atom in the system: for simplicity we followed the AMOEBA's atom type assignment. Secondly, ab initio atomic parameters such as isotropic static polarizabilities are extracted from subsystems (isolated water molecules as well as [Ala] 3 for this specific test case) and then averaged to obtain atom type-based parameters. This approach bypasses the DFT computation of the electron density for the whole system and its subsequent atom-in-molecule decomposition, which, even for a linear-scaling DFT approach, is prohibitive for a system of the size as the one here considered. This strategy allows for evaluating the energy of the 100000 atoms system in approximately 2 minutes' time on 400 cores from Intel Skylake 8168 (2.7 GHz) fine nodes. We emphasize that these remarkable performances have to be attributed not only to the linear-scaling and nearly embarassingly parallel implementation of the algorithm but also to the self-averaging of errors which, for this 100000 atoms system, allows us to get an estimated relative standard deviation of 0.5% corresponding to an error of 0.003kcal/mol per atom by only employing 300 random samples. In general, systems composed by few thousands atoms (or above) do not necessitate of large numbers of random samples in virtue of the algorithm's pleasant self-averaging of errors. On the other hand, different types of systems may exhibit different degrees of locality: this can be taken care of by choosing a suitable cut-off radius in building the elements of the MBD potential matrix thus tuning its sparsity. In general, delocalized and periodic systems are expected to exhibit longer-range interactions therefore a suitable and large enough cutoff radius must be carefully chosen. The ab initio nature of the input quantities employed by the MBD model as well as the contained number of natural amino acids would make it possible to generate rigorously a complete set of general parameters to be used in the modelling of proteins in solution where the plasmonic nature of the solvated system can require a many-body treatment of dispersion interactions. This will be the object of forthcoming studies coming along with a further reduction of MBD's time-to-solution thanks to Tinker-HP's accelerated multi-GPUs (Graphics Processing Units) implementation. 64