O (N) Stochastic Evaluation of Many-Body van der 1 Waals Energies in Large Complex Systems 2 Pier Paolo Poier †,* , Louis Lagard` ere †,‡ , and Jean-Philip Piquemal †,§,¶,* 3 * corresponding authors: Pier Paolo Poier (pier.poier@sorbonne-universite.fr) and Jean-Philip Piquemal 4 (jean-philip.piquemal@sorbonne-universite.fr) 5 † Sorbonne Universit ` e, LCT, UMR 7616 CNRS, Paris, France 6 ‡ Sorbonne Universit ` e, IP2CT, FR 2622 CNRS, Paris, France 7 § Institut Universitaire de France, Paris, France 8 ¶ The University of Texas at Austin, Department of Biomedical Engineering, TX, USA 9 ABSTRACT 10 We propose a new strategy to solve the Tkatchenko-Scheffler Many-Body Dispersion (MBD) model’s equations. 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. 11 Introduction 12 Van der Waals dispersion interactions play a crucial role in nature, determining the structure and functionality of molecular 13 systems and materials, 1, 2 two examples being the deviation from the ideal gas law observed in real gases and the ability of 14 geckos to stick to walls. 15 These inherently quantum-mechanical interactions originate from the non-local electronic correlation among distant densities 16 and their accurate description is crucial for a reliable computational modelling of matter. As a consequence of the underlying non- 17 local electronic correlation, Kohn-Sham (KS) Density Functional Theory (DFT) is unable to model dispersion interactions 3–6 18 and expensive wavefunction-based methodologies are usually required. In this connection, three-body Symmetry Adapted 19 Perturbation Theory (SAPT) approaches have been shown to provide very accurate evaluations of dispersion energies. However, 20 these methods are computationally challenging for large systems with scalings ranging from O (N 7 ) for HF-SAPT 7, 8 to O (N 5 ) 21 for SAPT(DFT) 9 . 22 To retain the attractive performances of DFT methods, several dispersion-corrected functionals have been proposed where 23 dispersion is added via empirical pairwise London contributions. 10 24 This simple and computationally appealing approach has proven successful in including dispersion effects within the DFT 25 framework although the semiempirical nature of exchange-correlation functionals may have effects on their overall accuracy. 11 26 One limitation of the meanfield pairwise approach in modelling dispersion interactions is the impossibility of capturing 27 many-body dispersion (MBD) effects: many evidences have in fact demonstrated the importance of a non-additive treatment 28 of dispersion interactions in modelling supramolecular complexes 12 , clusters of atoms 13 , one-dimensional wires 14 , extended 29 systems 15 as well as in molecular crystals. 16–18 30 The non-additive long-range character of dispersion interactions has been modelled, among the many different approaches 19 , via 31 a set of coupled fluctuating dipoles (CFD) 20, 21 or alternatively by quantum Drude oscillators 22–25 , the former being preferred 32 over the latter for computational reasons. 26 33 In recent years, Tkatchenko et al. have formulated the MBD@rsSCS model based on the CFD and introducing the range- 34 separation of the self-consistent screening (rsSCS) of polarizabilities which are related to each of the quantum harmonic 35 oscillator in the system. 27 In particular, the model relies on ab initio derived input quantities thus almost removing the 36 presence of empirical parameters. 28 Moreover, the self-consistent screening procedure improves the physical description of the 37 polarizability anisotropy which strongly depends on the atomic environment. 38