Generalized Many-Body Dispersion Correction through Random-phase Approximation for Chemically Accurate Density Functional Theory

We extend our recently proposed Deep Learning-aided many-body dispersion (DNN-MBD) model to quadrupole polarizability (Q) terms using a generalized Random Phase Approximation (RPA) formalism, thus enabling the inclusion of van der Waals contributions beyond dipole. The resulting DNN-MBDQ model only relies on ab initio-derived quantities as the introduced quadrupole polarizabilities are recursively retrieved from dipole ones, in turn modelled via the Tkatchenko-Scheffler method. A transferable and efficient deep-neuronal network (DNN) provides atom in molecule volumes, while a single range-separation parameter is used to couple the model to Density Functional Theory (DFT). Since it can be computed at a negligible cost, the DNN-MBDQ approach can be coupled with DFT functionals such as PBE,PBE0 and B86bPBE (dispersionless). The DNN-MBQ-corrected functionals reach chemical accuracy while exhibiting lower errors compared to the DNN-MBD dipole-only counterparts as well as to other MBD-based dispersion correction models where the accuracy gain can reache up to 45%.


Introduction
The importance of modelling matter through computer simulations has risen tremendously in the past decades in virtue of both the theoretical achievements and the advent of mass-produced computers whose performances have increased exponentially.
Despite these tremendous achievements, the exact solution of the non-relativistic electronic Schrödinger equation remains out of reach for multi-electron systems and different approximations have been introduced in order to model and tackle systems of chemical relevance. In particular, Kohn-Sham Density Functional Theory (KS-DFT) has establish itself the most widely used electronic structure method being the cheapest way for introducing electronic correlation. KS-DFT is based on the idea of evaluating the kinetic energy from a single Slater determinant, thus assuming the electrons being non-interacting.
The difference between the true and KS kinetic energies as well as the difference between the true electronic and exchange interaction ones is essentially embedded and modelled by the exchange correlation functional, key quantity in KS-DFT, however, in practice unknown.
The different approaches taken in modelling the exchange-correlation functional define the plethora of KS-DFT variants which, however, are mainly capable of capturing local correlation effects. Dispersion interactions, on the other hand, are rooted in the longrange electron correlation that can clearly not be captured by the intrinsic locality of common exchange correlation functionals. This inadequacy of KS-DFT in modelling noncovalent dispersion interactions, ubiquitous in nature and materials, has risen the attention towards the development of dispersion correction models, 1 the most widely used being based on empirical pairwise terms, Eq.(1). 2 This pairwise approach coupled with KS-DFT has shown to provide very good accuracy despite its very simple nature that adds basically negligible computational time and, for this reason, it is also employed in the majority of force fields as attractive component of the Lennard-Jones potential.
Each of the above pairwise terms can be further expanded via a second-order perturbative approach 3 (in the limit of large inter-atomic separation) to include the higher-order contributions shown in Eq. (2) where α j k (iν) is the k-pole polarizability at imaginary frequency of atom j (only the single pair ij is considered for the sake of simplicity). 4 While providing the correct long-range asymptotic limit, C 6 terms alone are not enough to describe the short-and medium-range dispersion effects and these higher-order terms have proven to increase significantly interaction energies near equilibrium regions as well as condensed phase properties obtained with molecular dynamics simulations based on force field potentials. 5 Regardless of the higher-order terms inclusion discussed, this pairwise approach completely neglects the collective many-body dispersion (MBD) effects inherited from the intrinsic quantum mechanical nature of long-range electron correlation and their relevance has been proven in modelling extended systems, supramolecular complexes and proteins in solutions. [6][7][8][9] These non-additive dispersion effects have been modelled via a set of coupled fluctuating dipoles 10,11 as well as quantum Drude oscillators. 12

Theory
The correlation energy E c for a system of interacting electrons can be rigorously expressed with the adiabatic connection fluctuation-dissipation formula in Eq. (3), v being the Coulomb potential |r − r | −1 whose coupling strength is governed by λ while χ λ and χ 0 are the interacting and non-interacting response functions where the latter can be evaluated from a set of single particle orbitals φ i (r) with their corresponding energies i and occupation numbers f i via the Adler-Wiser formalism, Eq.(4). 20,21 By defining χ λ (r, r , iν) via the self-consistent screening Dyson equation, Eq.
In virtue of its generality, the framework here discussed for the case of electronic correlation, can be used in different contexts and in particular in connection with the many-body dispersion model where the target energy arises from the correlation of coupled fluctuating dipoles. 23 In this case the response function assumes the form of the frequency-dependent atom-in-molecule (AIM) isotropic dipole polarizability α µ i (iν) localized at the atomic position R i while the Coulomb potential v(r, r ) is replaced by the Eq.(6) (last equality) is thus written in terms of the MBD model's quantities as shown in Eq. (7), A(iν) being the (3N, 3N) diagonal superpolarizability matrix having atomic polarizabilities as elements and T LR being the full dipole-dipole interaction tensor excluding self interactions (T LR ii = 0). We note that Eq. (7) does not include the χ 0 v term appearing in Eq.(6) due to the traceless property of the A(iν)T LR matrix product.
For MBD models based on coupled fluctuating dipoles, the RPA formula in Eq.(7) does not represent the most efficient strategy to compute correlation energies solution as the problem can be equivalently solved via exact diagonalization of the MBD potential matrix for a system of coupled quantum harmonic oscillators. 23 We will now focus our attention on A and, in particular, on its dipole and quadrupole dynamic polarizability entries being the model's key parameters.
As mentioned earlier, a pleasant feature of the MBD model is the fact that it only relies on ab initio-derived parameters. In particular, AIM dipole polarizabilities are typically obtained via the volume rescaling approach shown in Eq.(9), α µ i and V i being the AIM static dipole polarizability and volume respectively of the i-th atom while the zero superscript denotes free-atom reference quantities.
This volume/polarizability proportionality was discussed by Brinck et al. 24  polarizability based on single-particle systems' analysis. 28 We will, in the present work, stick to the common volume/polarizability proportionality in Eq.(9) as this will allow us to make ready comparisons with our dipole-only DNN-MBD model.
The AIM volume V i is commonly accessed by Hirshfeld partitioning the explicit electron density obtained via the solution of KS equations. Instead, we will adopt the recently proposed efficient and accurate 5-hidden layers DNN model trained with the minimal basis iterative Stockholder atom 29 (MBISA) volumes of approximately 4.6 millions molecules that provides AIM volumes bypassing the electron density partitioning. 18 Compared to Hirshfeld partitioning as well as its Iterative HI variant more suitable for ionic systems, 15 the MBISA scheme used to generate the training set does not suffer from asymmetric AIM densities arising from the use of free-atom reference densities that, especially for hydrogen atoms, lead to an overestimation of high radial moments, i.e. AIM volumes. 29 A further disadvantage related to the use of free-atom reference densities affecting the HI partitioning arises from the density interpolation for negatively charged atoms. In fact this procedure is, for some anionic species and for all doubly charged anions typically encountered in inorganic oxide clusters, ill-defined.
The next step is to find a suitable expression for α Q i (0), possibly without the introduction of empirical parameters to preserve the original MBD spirit. At this point, it is possible to relate two consecutive dispersion coefficients by the scaled Starkchal-Gordon relation 30 shown in Eq.(10) for the specific case of homonuclear C ii 6 and C ii 8 coefficients, where Z i is the atomic number, r n are expectation values or multipole-type moments derived from atomic densities while γ scaling factor will be introduced shortly. 30 We note, in passing, that Eq.(10) represents a fundamental recursion relation for Grimme's successful DFT-D3 dispersion correction model. 31 In order to retrieve α Q i (0) in terms of known quantities, we can now assume to model atom i via a quantum Drude particle characterized by the frequency ω i as this allows us to express C ii 6 and C ii 8 in terms of dipole and quadrupole polarizabilities, Eqs. (11) and (12). 32 By inserting Eqs. (11) and (12) in Eq. (10), it is finally possible to isolate the static AIM quadrupole polarizability in Eq. (13).
The combination of the scaled Starkchal-Gordon relation with a quantum Drude particle expression for the C ii 8 dispersion coefficient was originally adopted by Carter-Fenk et al.
in their extended symmetry-adapted perturbation theory with MBD (XSAPT+MBD) approach including pairwise dipole-quadrupole dispersion effects. 33 In that occasion, the γ scaling factor was successfully modelled according to Eq. (14), where γ 0 was chosen such that noble-gas quadrupole polarizabilities are reproduced and the same strategy is embraced in this work.
A further discussion about the combination of the Starkchal-Gordon rule with a quantum Drude particle expression is found in the "Annex" section at the end of the manuscript.
The above discussed recursive relation represents an efficient strategy to derive quadrupole polarizabilities that can be readily applied to any AIM volume partitioning scheme as the atomic volume is the only density-related quantity which, however, in this work is retrieved from an accurate DNN model.
In particular, we employ a common ω M,0 k parameter for both dipole and quadrupole dynamical polarizabilities and we express it as a function of free atom C 6 and static dipole polarizability, Eq. (16).
This assumption is, however, well defined as it guarantees to model dynamical quadrupole polarizabilities as lower bound quantities while the introduction of ab initio-derived scaling factors could be used to improve over this assumption. 34 We note that the dipole polarizabilities used in Eq.(13), can be either taken directly from Eq.(9) or can be further screened by solving Dyson-like self-consistent screening equations. 13,35 In this work the latter is chosen and, as a consequence, these screening effects are thus transferred to quadrupole polarizabilities via Eq. (13).
The explicit expressions of the multipole interaction tensor T LR ik , is given in Eq.(17) for the ik pair placed at distance R ik and greek letters are used to denote x, y, z Cartesian components. The These l MM terms could in principle be set as parameters and optimized together with β to reduce a target energy error, however, to avoid overfitting these terms will, in this work, be chosen once for all regardless of the chosen functional and not optimised.
We solve the RPA formula in Eq.(7) by means of Gauss-Legendre quadrature, nevertheless, more sophisticated techniques such as the Clenshaw-Curtis quadrature 37,38 could also be potentially employed.
The DNN-MBDQ model is implemented in the Tinker-HP package 39 where the extremely efficient linear-scaling stochatic Lanczos-based DNN-MBD model is also implemented. 40 In order to design an accurate and reliable dispersion-corrected KS-DFT model, both the functional and the dispersion correction must fulfill specific criteria. These requirements are essential for avoiding error cancellation effects in favor of more physically grounded achievement of small errors, as recently discussed by Price et al.. 41 While the dispersion correction should possibly include higher-order terms as well as many-body effects, the DFT functional should be dispersionless since by adding a dispersion correction to an exchange-correlation functional that somewhat partially includes dispersion, will lead to double counting at shorter ranges, with a consequent loss of accuracy. In fact, in order to reduce this artifact, dispersion corrections are often excessively damped with an inevitable deterioration of the mid-range interactions' description. Among the possible dispersionless functionals, the B86bPBE has proven numerically stable and, coupled with the exchange-hole dipole moment dispersion model, 42,43 has performed well for molecular and material applications. 44

Results
The S66x8 dataset for noncovalent biologically relevant interactions 46 Table 1 together with the ones for dipole-only based DNN-MBD. l µµ , l µQ and l QQ factors are set to 1, 2.8 and 3.0 respectively as these have proven to be good and transferable values for different functionals. We observe that the PBE0/DNN-MBDQ method Table 1: Range-separation β parameter values optimized by minimizing the mean absolute relative error for the DNN-MBD and DNN-MBDQ models coupled with the PBE and PBE0 functionals for the S66x8 data set.     The dipole contribution to dispersion interactions dominates the long-range asymptotic limit while inclusion of quadrupole contributions improves the description at medium ranges near equilibrium distances. The S66x8 is an optimal dataset to optimize the rangeseparation parameters as it is not biased towards equilibrium distances, however, in virtue or its wide range of inter-molecular distances, it does not represent a fully optimal set to investigate the effects of quadrupole terms. In that regard, results on the S22 dataset, which is composed of model complexes at equilibrium distances, are a better judge of the accuracy of our model since S22 includes more interactions in the full quadrupole (and mixed dipole-quadrupole) operational regime. The S22 validation set thus enables us to better understand the effects of the added quadrupole terms as well as to investigate the transferability of the DNN-MBDQ model by employing the same range-separation parameters reported in Table (1). Table 4 shows the performances for some of the dispersion correction models reported in Table 2 for which results relative to the S22 set are avail- have been discussed profusely in references and previously recalled in this work. 18,29 Among MBD-based models, the DNN-MBDQ coupled with the functionals here consid-ered provides one of the lowest errors for the S22 set without an ad hoc parameter optimization. This is rather significant considering that the model fully relies on ab initio quantities. The accuracy performances of the DNN-MBDQ model could be potentially increased further by adding an extra parameter in the Fermi damping function (Eq. (18)), as well as by performing a complete optimization for the choice of the l MM terms and future efforts will be spent in this direction.

Conclusions
We propose the DNN-MBDQ model where the density-free/deep learning-aided manybody dispersion model is extended to quadrupolar polarizability terms thanks to a Random Phase Approximation formalism. Quadrupole polarizabilities are recursively retrieved from free-atom dipole polarizabilities and atom-in-molecule volumes modelled via the accurate and transferable neuronal network recently proposed. 18 The described stategy can be readily applied to MBD model based on the volume reascaling polarizability partitioning as it does not require any additional electron density-derived quantity.
The density-free DNN-MBDQ model implemented in the Tinker-HP package exhibits improved accuracy compared its dipole-only DNN-MBD counterpart as well as reference MBD-based models. In particular, for the widely used S22 data set of dimers placed at equilibrum inter-molecular distances, the DNN-MBDQ shows a remarkable improvement (up to 45%) compared to dipole-only reference MBD models.
In addition to the common and successful PBE and PBE semi-local functionals, the DNN-MBDQ model is also coupled with the B86bPBE one, to explore the possibility of employing dispersionless functionals, preventing double counting dispersion, especially at shorter ranges. For this case, the errors are abundantly below the chemical accuracy threshold and similar (in some cases even lower) to the ones obtained for the PBE/PBE0 case.
We believe that the high accuracy exhibited by the DNN-MBDQ model coupled with semi-local and dispersionless functionals will be beneficial in pushing Kohn-Sham Density Functional Theory a step closer to post-Hartree Fock reference standards which, due to their high computational cost can not be employed to explore larger systems of relevance in bio-and material modelling. Further work will be therefore dedicated to the

Supporting Information Available
SI can be download directly via the Zenodo repository located at the address in reference. 51 SI contains raw PBE, PB0 and B86bPBE interaction energy data for the S66x8 and/or S22 data sets as well as their relative DNN-MBDQ dispersion corrections for the different β values discussed.

Annex
An analytical connection between the Starkschall-Gordon rule and the quantum Drude oscillator (QDO) expression is non-trivial to obtain, however, to verify the consistency of the two expressions, it is possible to evaluate the C 8 /C 6 proportionality factor for the hydrogen atom, for which α Q , α µ , r 2 and r 4 are analytically available. 30,52 We start by taking the QDO expression for the C 8 coefficient where we introduce ω = It is now possible to evaluate f QDO and f SG from the analytical α Q , α µ and r 4 , r 4 expres-