O(N) Stochastic Evaluation of Many-Body van der Waals Energies in Large Complex Systems
We propose a new strategy to solve the key equations of 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 frequencydependent 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.
Introduction<!>Reviewing the MBD@rsSCS model<!>2/27<!>3/27<!>Linear-scaling iterative self-consistent screening equations<!>2<!>Stochastic Lanczos trace estimator's expressions for E int<!>6/27<!>Computational Details<!>Results and discussion<!>Iterative solution of the reduced rs-SCS Dyson-like equation<!>5:<!>Stochastic Lanczos based estimation of the MBD energy<!>Conclusion<!>Supporting Information
<p>Van der Waals dispersion interactions play a crucial role in nature, determining the structure and functionality of molecular systems and materials, 1,2 two examples being the deviation from the ideal gas law observed in real gases and the ability of geckos to stick to walls. These inherently quantum-mechanical interactions originate from the non-local electronic correlation among distant densities and their accurate description is crucial for a reliable computational modelling of matter. As a consequence of the underlying nonlocal electronic correlation, Kohn-Sham (KS) Density Functional Theory (DFT) is unable to model dispersion interactions [3][4][5][6] and expensive wavefunction-based methodologies are usually required. In this connection, three-body Symmetry Adapted Perturbation Theory (SAPT) approaches have been shown to provide very accurate evaluations of dispersion energies. However, these methods are computationally challenging for large systems with scalings ranging from O(N 7 ) for HF-SAPT 7,8 to O(N 5 ) for SAPT(DFT) 9 . To retain the attractive performances of DFT methods, several dispersion-corrected functionals have been proposed where dispersion is added via empirical pairwise London contributions. [10][11][12][13] This simple and computationally appealing approach has proven successful in including dispersion effects within the DFT framework although the semiempirical nature of exchange-correlation functionals may have effects on their overall accuracy. 14 One limitation of the meanfield pairwise approach in modelling dispersion interactions is the impossibility of capturing many-body dispersion (MBD) effects: many evidences have in fact demonstrated the importance of a non-additive treatment of dispersion interactions in modelling supramolecular complexes 15 , clusters of atoms 16 , one-dimensional wires 17 , extended systems 18 as well as in molecular crystals. [19][20][21] The non-additive long-range character of dispersion interactions has been modelled, among the many different approaches 22 , via a set of coupled fluctuating dipoles (CFD) 23,24 or alternatively by quantum Drude oscillators [25][26][27][28] , the former being preferred over the latter for computational reasons. 29 In recent years, Tkatchenko et al. have formulated the MBD@rsSCS model based on the CFD and introducing the rangeseparation of the self-consistent screening (rsSCS) of polarizabilities which are related to each of the quantum harmonic oscillator in the system. 30 In particular, the model relies on ab initio derived input quantities thus almost removing the presence of empirical parameters. 31 Moreover, the self-consistent screening procedure improves the physical description of the polarizability anisotropy which strongly depends on the atomic environment.</p><p>The MBD@rsSCS model can be summarized into three main steps. Firstly, a set of atom-in-molecule static dipole polarizabilities are obtained from the partitioning of the molecular electron density. Secondly, the polarizabilies are made frequency-dependent via Padé approximation and subsequently a Dyson-like self-consistent screening linear equation is solved for a selected set of frequencies. Lastly, the set of screened frequency-dependent polarizabilities are employed in the solution of Casimir-Polder integrals and the outcoming screened atomic charateristic excitation frequencies are then employed to compute the target MBD energy.</p><p>The MBD@rsSCS model coupled to DFT methods involves a small additional computational cost compared to the solution of the KS equations, however, for increasingly large systems where linear-scaling DFT methods are employed, the O(N 3 ) scaling of the diagonalization procedure as well as the solution of the self-consistent screening equations becomes no longer negligible. In the present work, we considerably extend the applicability of the MBD@rsSCS model via a novel computational resolution of its equations. A robust and efficient linear-scaling scheme for the iterative solution of the self-consistent screening equations is introduced together with a stochastic approach allowing for a linear scaling estimation of the MBD energy which bypasses matrix diagonalization. The article is organized as follows: first the MBD@rsSCS model is briefly reviewed at the beginning of the Theory section. Later in the same section the theoretical aspects of the efficient solution of rsSCS equations are presented as well as the theoretical foundations of the linear-scaling stochastic estimator of the MBD energy. A detailed discussion of the algorithms and of their massively parallel implementations are presented in the Numerical Result section that is followed by the Conclusion.</p><!><p>The MBD energy E MBD corresponds to the correlation energy for a system of coupled fluctuating dipoles, i.e. the difference between the system's full interaction energy E int and the zero-point (non-interacting) energy E 0 , Eq.( 1).</p><p>We start focusing on the zero-point energy which, for a N atoms system, is defined in Eq.( 2) where ωi represents the screened resonance frequency of the harmonic oscillators associated to each atom in the system.</p><p>To identify the set of { ωi } we start by introducing a corresponding set of bare isotropic and static atomic polarizabilities {α i (0)} which, in the original Tkatchenko-Scheffler scheme 31 , is obtained by partitioning of the total electron density via the Hirshfeld 32 decomposition or, alternatively, via an iterative Hirshfeld scheme. 33,34 The static polarizabilities are turned into frequency-dependent ones by employing the [0/2] Pade approximants' form in Eq. (3). ω 0 j is the free-atom resonant frequency which is defined in terms of the free-atom's first dispersion coefficient and static polarizability, C 0 6, j and α 0 j (0) respectively. We note in passing that both quantities can be obtained from accurate ab initio calculations.</p><p>The above frequency-dependent bare polarizabilities are then used to define the (3N, 3N) block diagonal bare super-polarizability matrix A(iν) where the j-th diagonal block is shown in Eq. (5).</p><p>The next step is to obtain the range-separated screened super-polarizability matrix Ā(iν): this is accomplished by solving the Dyson-like self-consistent screening equation below.</p><!><p>The T SR (3N, 3N) matrix represents a damped dipole-dipole interaction operator applied to the Coulombic interaction between two spherical Gaussian charge distributions. Eq.( 7) shows a general element of the i j block of T SR for the atom pair placed at distance</p><p>The attenuation length σ i j (iω) is defined below in terms of the width of the two Gaussian distributions which depend on the bare frequency-dependent dipole polarizability introduced in Eq. (3).</p><p>The explicit form of Eq.( 7) is shown in Eq.( 9), T σ τ i j being the frequency-independent classical interaction term between two point dipoles. For the sake of notation, the frequency-dependency is dropped.</p><p>The damping function f (R i j , S i j ) defined in Eq.( 10) is responsible for the range-separation of the dipole-dipole interaction which is modulated via the β parameter.</p><p>In the above equation, R vdw i represents the atom-in-molecule van der Waals radius for atom i which is obtained from the atom-in-molecule static polarizability α i (0), the free-atom static polarizability α 0 i (0) and the free-atom van der Waals radius R 0,vdw i , the latter defined for noble gases as half of the distance between two atoms where exchange-repulsion balances London's dispersion attraction.</p><p>The rs-SCS equation in Eq.( 6) can be solved by multiplying from the left both of its sides by A −1 and by subsequently inverting the resulting matrix, Eq. (12).</p><p>Screened atomic polarizabilities are then obtained from a partial contraction of Ā(iν) as shown in Eq. (13).</p><p>We note that the resulting set of screened polarizabilities { ᾱ j (iν)} are not isotropic due to the effect of self-consistent screening. Finally, the screened resonant frequencies used in the zero-point energy E 0 , are defined similarly to Eq.( 4): this time instead, screened quantities are used and the screened dispersion coefficients must be explicitly computed via the Casimir-Polder integral, Eq. (14).</p><!><p>The Casimir-Polder integral is efficiently solved via Gauss-Legendre (GL) quadrature where a suitable change of variable is performed in order to work with a finite interval of integration. [35][36][37] This requires Eq.( 12) to be solved for a N ν number of frequencies.</p><p>Moving now to the system's full interaction energy E int , we introduce the MBD range-separated Hamiltonian in Eq.( 15), d i = √ m i ξ i being the mass-weighted dipole moment displaced by the vector ξ i from its equilibrium position. 24,30 ĤMBD = − 1 2</p><p>The long-range dipole-dipole interaction is taken care of by the range-separated T LR matrix which i j block is shown in Eq.( 16), T i j being is the classical dipole-dipole interaction tensor introduced earlier.</p><p>The damping function f (R i j , Si j ) differs from the one introduced in Eq.( 10) for the use of screened van der Waals radii in expressing Si j , Eq. (17).</p><p>From the Hamiltonian in Eq.( 15) we identify the (3N, 3N) MBD interaction matrix V MBD in Eq.( 18) and its full spectrum of eigenvalues {λ i } is used to finally define the interaction energy E int .</p><p>By putting Eqs.(1),( 2) and ( 19) together we can finally define the explicit MBD correlation energy expression, Eq. (20).</p><p>The overall computational cost involved in the evaluation of the MBD correlation energy arises from the N ν -times solution of Dyson-like self-consistent screening equations (Eq.( 6)) as well as from the eigendecomposition of V MBD both of which being as O(N 3 ) in scaling. In the following part of this section we first describe a robust and linear-scaling iterative self-consistent screening algorithm as well as a linear-scaling stochastic estimation of E int based on a state-of-the-art linear algebra numerical algorithm.</p><!><p>The set of screened polarizabilities ᾱ j (iν), key quantities in the aforedescribed MBD@rsSCS model, are typically obtained by partial-contraction of the screened super-polarizability Ā(iν). We propose, instead, to first contract Dyson-like rs-SCS equation prior its iterative solution which directly provides the set of screened polarizabilities. This approach not only involves a reduced number of variables with the consequent reduction of memory requirements but also allows us to obtain screened polarizabilties in a linear-scaling fashion as it will be discussed profusely.</p><p>We begin by introducing the full (3N, 3) contraction block-vector operator C which j-th block elements consist of the (3, 3) identity matrix I 3 , Eq. (21).</p><p>By projecting C on both the left sides of Eq.( 6), we obtain a partially contracted rs-SCS set of equations, Eq.( 22), where P and P represent the (3N, 3) arrays gathering all bare and screened atomic polarizabilities respectively. ĀC = AC − AT SR ĀC P = P − AT SR P ( 22)</p><p>We note that the partial contraction has, for a given frequency, reduced the number of variables from 3N(3N+1)</p><!><p>(nonredundant elements of Ā) to 9N (elements of P). Eq.( 22) is simplified by splitting P (and P) column-wise, Eq.( 24): the resulting three independent linear sets of equations are then solved iteratively as shown in Eq.( 25), [k] denoting the iteration.</p><p>p[k+1]</p><p>The overall computational cost of Eq.( 25) arises from the AT SR p[k] σ matrix-vector products involved in the iterative procedure, where the diagonality of A allows the computation of AT SR on the fly at nearly no cost. The matrix-vector products described above can be performed very efficiently by exploiting the short-range nature of T SR , i.e. its sparsity. This is achieved by building specific short-range neighbor lists based on a cut-off radius which choice will be discussed and motivated in the numerical results section. A further aspect to consider in improving the matrix-vector products' efficiency is the way the T SR elements are built. We remind that Eq.( 25) should be solved as many times as the number of frequencies (N ν ) needed in the GL quadrature employed to solve Eq. (14). In this regard it is possible to compute the purely geometry-dependent parts of T SR once and this allows their re-use in the on-the-fly building of frequency-dependent elements. This, as it will be analyzed later, decreases considerably the computational cost. Efficient matrix-vector products would be of little help if the overall convergence required a large number of iterations: it is therefore important to make use of a robust and efficient iterative scheme. The simple self-consistent solution of Eq.( 25) can be considered as a fixed-point iterative procedure 38 where the previous iteration's result is used as input for the next one. The rate of convergence of fixed-point iterative algorithms can be generally improved by constructing, at each iteration, a guess vector being a linear combination of previous iterations' resulting vectors where the coefficients of the linear combination are determined to minimize in a least squares sense the error. This extrapolation scheme known as diret inversion of the iterative subspace (DIIS) was introduced with the aim of accelerating and stabilizing the convergence of the Roothaan-Hall self-consistent field equations. 39 The application of the DIIS together with fixed-point iterative methods in solving linear systems of equations makes the overall procedure equivalent to the Generalized Minimal Residual method (GMRES) which, belonging to the Krylov subspace class of iterative methods, is known for its fast and robust convergence properties. 40 We mention that the coupled fixed-point/DIIS iterative schemes have proven successful in the context of modelling many-body polarization interactions [41][42][43] and this motivates the use of this mixed procedure in the solution of Eq. (25). Finally, when the px , py and pz arrays are converged for all the N ν frequencies, the frequency-dependent screened atomic polarizabilities are isotropized (second equality in Eq.( 13)) and the set of { ωi } constituting the zero-point energy E 0 are obtained via Eqs.( 14) and (2) respectively.</p><!><p>The E int defined in Eq.( 19) requires the full eigendecomposition of V MBD , Eq.( 26), which, especially for large systems, becomes computationally expensive.</p><p>We propose, in the following, an alternative approach bypassing the diagonalization step. We first express E int in terms of the trace of √ Λ: by recognizing the trace invariance under any change of basis, the interaction energy is written as shown in Eq.( 27) which, at this point is still exact. For simplicity, the MBD lower script is dropped from now on.</p><p>The E int expression is now turned into the evaluation of the trace of a square-root matrix function: this very problem has been recently tackled by Ubaru, Chen and Saad by means of a stochastic Lanczos (SL) trace estimation. 44 For the sake of generality, the algorithm is discussed for Tr [ f (V)], where f denotes a general analytic function inside a closed interval.</p><p>To estimate Tr [ f (V)], the proposed algorithm inkoves Hutchinson's stochastic trace estimator 45 , Eq.( 28), where v l is one of the R normalized random vectors of dimension D (in our case D = 3N), which elements follow a Rademacher distribution, i.e. they can assume values of either 1 or −1 with equal probability.</p><p>Each of the R scalar expectation values in Eq.( 28) can be expressed in terms of f (Λ) and the unitary transformation W as reported in Eq.( 30) where we introduced</p><p>The last equality in Eq.( 30) is found to be the Riemann-Stieltjes integral 46 defined in Eq. (31).</p><p>The estimation of each v † l f (V)v l scalar thus corresponds to the solution of a Riemann-Stieltjes integral which can be approximated by means of quadrature techniques, {τ k } and {θ k } representing the unknown weights and nodes respectively of the general (M + 1)-point quadrature.</p><p>By inserting Eq.( 32) in Eq.( 28), one can identify the complete expression for the stochastic trace estimation, Eq. (33).</p><!><p>The stochastic Lanczos algorithm, as the name suggests, makes use of Lanczos quadrature in approximating the Riemann-Stieltjes integrals 46 . Below we provide and motivate the explicit expressions for nodes and weights within the Lanczos quadrature framework.</p><p>The starting idea in the discussion of the quadrature is the generation, via the Lanczos algorithm reviewed in the Supporting Information, of an orthonormal basis for the (M + 1)-dimensional Krylov subspace, where v l is taken as starting vector.</p><p>The orthonormal basis vectors are compactly gathered as columns of the Y matrix of dimension (D, M + 1) which allows for transforming V into the tridiagonal (M + 1, M + 1) ∆ (l) matrix, Eq. (35).</p><p>The diagonalization of ∆ (l) via the unitary transformation U (l) generates Λ (l)</p><p>, Eq.( 33), which corresponds to the full spectrum of V projected onto K (l) M+1 (V, v l ). In the equation below, the equality</p><p>We are now in the position of introducing the essential point of Lanczos quadrature, namely the approximation of v † l f (V)v l by means of the just introduced quantities Λ (l) and W (l) respectively.</p><p>)</p><p>Eq.( 37) simply represents the (M + 1)-point quadrature approximation of Eq.( 30). Moreover, we note that Eq.( 37) involves scalar products between v l and all the basis vectors gathered in Y spanning K</p><p>M+1 : being the basis orthonormal, the only non-vanishing element of the resulting array is the one corresponding to its first entry, Eq. (38).</p><p>The insertion of Eq.( 38) into Eq.( 37) finally allows us to identify the nodes and weights being the eigenvalues λ</p><p>k and the squares of the first entries of the eigenvectors [U (l)</p><p>An alternative derivation of the nodes and weights expressions based on the properties of Lanczos polynomials can be found in reference. 46 We now replace the general nodes and weights in Eq.( 32) with the just derived ones to express the stochastic Lanczos estimation of Tr [ f (V)].</p><p>The overall computational cost of the SL algorithm is to be mainly addressed to the construction of the tridiagonal matrices ∆ (l) via the Lanczos algorithm which, however, only relies on matrix-vector products and can thus be performed in a linear-scaling fashion. On the other hand, the computational cost required by the diagonalization of ∆ (l) is negligible as the (M + 1) number of quadrature points (and therefore the dimension of the tridiagonal matrix to be diagonalized) is very small, typically not exceeding 50, regardless of the system size. Moreover, the exact diagonalization of small tridiagonal matrices can be performed efficiently by means of several O(N) efficient libraries. We finally note that the SL algorithm is embarrassingly parallelizable with respect to the R random vectors as they involve independent processes thus requiring minimal communication.</p><p>7/27</p><!><p>The presented algorithms have been implemented in the massively parallel Tinker-HP molecular dynamics package 47 . Atomsin-molecule quantities were computed with the Multiwfn software 48 by partitioning DFT electron densities generated with the Gaussian quantum chemistry software 49 at the PBE/cc-pVQZ 50,51 level. Free-atom parameters are taken from reference 52,53 except for noble gases' van der Waals radii which have been taken as the half of the equilibrium distance obtained from accurate reference potential energy surfaces. [54][55][56] .</p><!><p>Before the discussion and analysis of the newly proposed algorithms, we validated the implementation of the standard model which will then represents our target reference. This was done by selecting, from the S66 57 dataset, 9 structures for which dispersion plays a major role in inter-molecular interaction energies. For these selected dimers the implemented model generates MBD energies being in agreement (average deviation of 0.4 kcal/mol, Supporting Information) with the ones computed with the MBD correction implementation included in the Q-Chem package. 58 We note however, that the Hirshfeld partitioning can give as output slightly different parameters according to the software implementation employed thus leading to potentially small discrepancies.</p><!><p>We start by analyzing the iterative self-consistent solution of the reduced Dyson-like screening equation from which the set of frequency-dependent screened polarizabilities are obtained and subsequently used in the solution of the Casimir-Polder integral in Eq.( 14). Before proceeding with the numerical analysis it is first important to understand the adequate number of quadrature points (i.e. frequencies) to be used in the GL quadrature. While a too low number of quadrature points may lead to numerical errors, a too large one, on the other hand, would lead to an unnecessary increase of the computational cost as well as memory requirements and a balance must be therefore found. The analysis reported in the Supporting Information shows that a ten-points GL quadrature leads to errors in the order of 10 −6 and this number of choice will be embraced throughout the discussion.</p><p>We are now in position to analyze numerically the proposed algorithms, starting from the convergence rate for the iterative solution of Eq.( 25) which is schematically summarized below. In this and in the following section, argon clusters are taken as test systems as this only involves free-atom parameters and decouples the MBD model from the Hirshfeld partitioning, thus allowing to easily work with increasingly large system sizes without facing the DFT computational bottleneck.</p><p>Algorithm 1 Schematic representation of the iterative solution of Dyson-like rsSCS equation. G i j and F i j (iν n ) are used to represent the geometry-and frequency-dependent parts of the T SR i j (iν n ) block respectively .</p><p>for ( j = neighbor(i)) do</p><!><p>Build G i j 6:</p><p>end if 10:</p><p>Assemble</p><p>end for</p><p>The rate of convergence is studied on a cluster of 30 argon atoms where the inter-atomic distance is, on average, less than the half of the argon-argon equilibrium distance. This high-density system allows for the testing the algorithm's robustness under high screening regimes. Figure ?? shows the convergence of the screened frequency-dependent polarizabilities in terms of the logarithm (base 10) of the residual defined in Eq.( 41) for a single frequency ν n , as function of the number of iterations in the case of a simple fixed-point solution of Eq.( 25) where line 18 in Algorithm 1 is skipped.</p><p>First we note that the rate of convergence is different for different frequencies even though a possible pattern has not been identified. This confirms the importance of monitoring the convergence in a frequency-dependent manner in order to avoid unnecessary iterations for those frequencies already below a given convergence threshold. Moreover the convergence of the simple fixed-point algorithm suggests that the rs-SCS problem is well-conditioned even for the dense cluster of atoms taken as test system. It is also of interest to probe the effects of the DIIS extrapolation coupled to the fixed-point iterative solution of the screening equations and this is shown in Figure 2. As expected, the DIIS extrapolation accelerates the convergence, in particular for the slow converging frequencies observed previously in Figure ??. The total number of iterations is reduced, for the full set of frequencies, by 55%, underlining the benefits of the DIIS inclusion. As mentioned in the theoretical section, an additional aspect considered to improve the performance of the algorithm is the re-use of purely geometric-dependent terms (referred to as G i j in Algorithm 1) in the on-the-fly building of T SR i j (iν n ): this allows for further decreasing the computational cost by up to 40% for a 10-points quadrature as shown in Figure 3. A final aspect concerning the efficient solution of the rs-SCS sets of equations is the choice of the proper cutoff radius from which the neighbor lists can be built. The short-range nature of T SR responsible for its sparsity can be seen in Figure 4 where the isotropized rs-SCS polarizability of one of the argon atoms in a dimer is plotted as function of the inter-atomic distance. At inter-atomic distances above 3 Ångstroms the effects of the self-consistent screening becomes negligible suggesting that a cutoff radius of 4 Ångstroms is more than enough to include all the anisotropic effects arising from the neighboring atoms. We wish add to the discussion of the proposed iterative rs-SCS algorithm described in Algorithm 1 the comparison between its computational scaling and the one relative to the direct inversion approach in Eq.( 12), Figure 5. The remarkable computational gain introduced by the iterative solution of the reduced rs-SCS Dyson-like equations shown in Figure 5 is, however, only one of the advantages introduced: a further aspect to take into account is in fact the memory efficiency. The solution of Eq.( 6) requires to store the upper-triangular part of Ā(iν) for all the set of frequencies employed, thus involving quadratic-scaling memory requirements with the system size. The variable reduction involved by Eq.( 25), on the other hand, only requires 9N elements to be stored for each of the frequencies employed, making the needed memory allocation linear scaling with the system size as shown in Figure 6. The advantages introduced by the here proposed strategy thus concerns both the computational performance and memory allocation requirements.</p><!><p>We now turn our attention to the proposed Stochastic Lanczos trace estimator, summarized in Algorithm 2, which is employed to estimate the many-body interaction energy (E int ). Lanczos iteration: y</p><p>for (k = 1, M + 1) do 10:</p><p>end for 12: end for 13: Calculate the average over samples: Tr [</p><p>The detailed description of Lanczos algorithm with explicit recursive expressions for y k , ∆</p><p>kk and ∆ (l) (k−1)k are given in the Supporting Information. Being the algorithm stochastic in its nature, it is of interest to analyze the distribution of the energy estimations in terms of the number of random samples (R) as well as the dimension of the Krylov subspace identifying the number of Lanczos quadrature points. This is shown in Figure 8 where three distributions relative to three different values of R are reported for a Krylov subspace of dimension 20. The SL distributions are centered close to zero and this denotes a negligible bias toward the target energy value obtaineded via 11/27 Figure 3. Cpu-time as function of the number of quadrature points in case of a complete rebuild of T SR i j (iν n ) (red plot) and with the re-use of purely geometrical-dependent terms (black plot) a direct diagonalization. Moreover, the distributions narrow with the increment of R as one would expect, thus allowing for a systematic improvement of the estimation. We note that by decreasing the dimension of the employed Krylov subspace to 10, the distributions remain well centered around zero with a bias smaller than 1% on the total MBD energy even for low values of R. A more complete analysis of the standard deviation dependence on the number of random samples R is shown in Figure 9 for two different Krylov subspace dimensions, 10 and 20 respectively: this reveals a typical Monte Carlo 1/ √ R scaling where the dimension of the Krylov subspace (i.e. the number of Lanczos quadrature points) has very little effects unless unresonably small subspaces are employed. We note at this point that Lanczos quadrature is extremely efficient in approximating the terms reported in Eq.( 32) as it requires half of the points, for given target accuracy, compared to other quadrature strategies such as the ones based on Chebyshev polynomials. This, together with a rigorous mathematical error analysis of the SL trace estimator, is discussed in the original work of Ubaru, Chen and Saad. 44 Figure 10 shows the evolution of the stochastic trace estimation with the number of random samples for a single point estimation where the trace approaches the target value with a standard deviation magnitude becoming smaller with the increase of the number of random samples, in line to what shown in Figure 9. A further aspect to probe is the standard deviation's dependency on the system size. This is shown in Figure 11 for increasingly larger argon clusters. The trend shows how the error decreases with the increase of the system size thus denoting a desirable statistical self-averaging of errors, in line with the fact that, in the thermodynamic limit, fluctuations die out. 59 We note in passing that a similar behavior was observed in stochastic approaches to compute electronic correlation energies. 60 Further words on this aspect will be spent later on in this section. It is now of interest to analyze the scaling with system size of the proposed SL trace estimation algorithm and to compare it with the standard approach requiring the full eigendecomposition of V. The building of the Krylov subspace only requires matrix-vector products which are performed in a linear-scaling fashion in contrast to the standard full-diagonalization approach which is cubic scaling with the system size as shown in Figure 12. We stress that the plot was generated by keeping the number of random samples constant, however, in virtue of the self- averaging property of the algorithm the number of random samples could have been notably reduced (for a given target accuracy) along with the system size and this would have contributed to an even better scaling of the SL algorithm's CPU-time with system size.</p><p>The stochastic algorithm presented here is by nature embarassingly parallel by dividing the random samples among the available processes while the only communication required is a reduction occurring when each of the independent processes has computed the contribution to the trace associated with its subset of random vectors. Though not being the computational bottleneck, the resolution of the Dyson equations for a set of frequencies is also parallelized with the spatial decomposition strategy embraced within the Tinker-HP package. Similarly to the resolution of the polarisation equations characterizing the AMOEBA force field, each process computes the contribution of the matrix-vector product due to its local p component. At each iteration, in order to compute this matrix vector product, a communication of p at the previous iteration has to be made before the computation and a communication of the result of the product has to be made after. Besides, reductions have to be made to compute the residual (Eq.( 41)) and the DIIS matrix at each iteration. Because the parallelization scheme is not the same between the iterative solution of the screening equations and the stochastic Lanczos trace estimation, a global broadcast of the screened polarizabilities has to be done between these two steps. Overall, the moderate amount of communication explains the small deviation from the ideal scaling that is observed in Figure 13 after 200 cores.</p><p>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 14. Firstly atom types are assigned to each atom in the system: for simplicity we followed the AMOEBA's atom type assignment. 61 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. 62,63 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 14/27</p><!><p>In this work we proposed two strategies to improve the scalability of the range-separated self consistent screening Many-Body Dispersion model. The first strategy focuses on the efficient solution of the Dyson-like self-consistent screening equations from which a set of frequency-dependent screened polarizabilities are obtained. In particular, the proposed algorithm drastically reduces the number of variables involved and relies on matrix-vector products which can be performed efficiently and in a linear-scaling fashion with minimal memory requirements. The algorithm is coupled to the DIIS extrapolation which not only reduces by 55% the number of iterations but also guarantees convergence thus making the algorithm efficient and robust at the same time. Secondly, we introduced an alternative stochastic strategy to estimate the many-body interaction energy of coupled quantum harmonic oscillators characterizing each atom in the system. This approach is based on the state-of-art Stochastic Lanczos trace estimation and exhibits a linear-scaling with the system size as well as an embarassingly parallel implementation arising from the minimal communication required among processes. Compared to the exact diagonalization approach, this strategy introduces negligible errors also in virtue of their the self-averaging characterizing the algorithm itself. Our approach, in fact, does not modify the physics and the underlying accuracy of the original Tkatchenko-Scheffler MBD model, rather it reformulates the resolution of the key equations in a O(N) fashion. We wish to highlight that the Stochastic Lanczos trace estimation algorithm, here profusely analyzed and discussed in connection to the MBD model, can be readily applied in the context of electronic structure theory and specifically in the Random Phase Approach (RPA) where the correlation energy is expressed in terms of a trace of some matrix-products or matrix-functions. The contained number of random samples required by the SL approach, compared to other stochastic Monte Carlo techniques employed to solve the RPA equations 60 , is likely to speed up by orders of magnitude the evaluation of electronic correlation energies. The algorithm's scalability opens up the computation of accurate many-body dispersion energies within a few minutes' time for up to millions-atoms' systems. This finally enables to investigate the role of physically grounded dispersion effects in complex materials and biologically relevant processes, which, due to the cubic scaling of standard approaches, are still mostly unexplored at realistic scales. The present implementation will be available in the next release of the Tinker-HP molecular dynamics package 47 and it will also be included in the Libmbd 65 external open-source library. Efficient O(N) DFT models represent the natural context of application for the newly presented algorithm, however, in virtue of its remarkable performances and making use of ab initio-derived quantities, the algorithm can give the basis for the rising of next-generation ab initio-based accurate force fields where both polarization and dispersion are consistently treated at the many-body level. 66,67 In that connection, a forthcoming study dedicated to the derivation of the MBD energy gradients will follow the present work.</p><!><p>An analysis on the Gauss-Legendre quadrature is provided together with a table summarizing the energy values for the structures used in the validation of the implementation. Finally the Lanczos algorithm is reviewed and explicit recursive expressions are reported.</p>