A Linear-Scaling Method for Noncovalent Interactions: An E ﬃ cient Combination of Absolutely Localized Molecular Orbitals and a Local Random Phase Approximation Approach

: A novel method for the accurate and e ﬃ cient calculation of interaction energies in weakly bound complexes composed of a large number of molecules is presented. The new ALMO+RPAd method circumvents the prohibitive scaling of coupled cluster singles and doubles while still providing similar accuracy across a diverse range of weakly bound chemical systems. Linear-scaling procedures for the Fock build are given utilizing absolutely localized molecular orbitals (ALMOs), resulting in the a priori exclusion of basis set superposition errors. A bespoke data structure and algorithm using density ﬁ tting are described, leading to linear scaling for the storage and computation of the two-electron integrals. Electron correlation is included through a new, linear-scaling pairwise local random phase approximation approach, including exchange interactions, and decomposed into purely dispersive excitations (RPAxd). Collectively, these allow meaningful decomposition of the interaction energy into physically distinct contributions: electrostatic, polarization, charge transfer, and dispersion. Comparison with symmetry-adapted perturbation theory shows good qualitative agreement. Tests on various dimers and the S66 benchmark set demonstrate results within 0.5 kcal mol − 1 of coupled cluster singles and doubles results. On a large cluster of water molecules, we achieve calculations involving over 3500 orbital and 12,000 auxiliary basis functions in under 10 min on a single CPU core.


INTRODUCTION
Chemistry relies on the idea that results obtained for molecules in one context are transferrable to other, similar situations. The idea of a molecule implies there is an inherent locality to the underlying description. While no such physical distinction truly exists, it is undeniably helpful to compartmentalize and categorize chemical systems. The latter is especially true in the case of noncovalent interactions, the ubiquitous forces underlying phenomena as diverse as geckos attaching themselves to walls, 1,2 to the cohesion of asteroids, 3 and more. 4,5 Such interactions are defined to be any stabilization weaker than would be expected for a chemical bonda definition so broad as to be insurmountably difficult to study without distinguishing further. Typically, such a distinction is found through decomposition of the interaction into terms such as electrostatic, polarization, and dispersion; 6 the complex is then typified by the predominant component.
Noncovalent interactions are difficult to study, both experimentally and theoretically. Isolating weak forces in the condensed phase is problematic, 7 although progress has been made in the use of molecular balances. 8 Similarly, attempts in the gas phase are hindered by the fact that the thermal energy may overcome the weak binding. Theory therefore is of vital importance to understand these weakly bound systems.
However, the presence of small energy differences and the importance of dispersionan effect entirely due to dynamical electron correlationmean that only the most accurate computational methods perform consistently well. 9−11 The lack of a reasonable description of dispersion, in particular, excludes the majority of density functionals, which have been shown to often give qualitatively and quantitatively inaccurate results, although progress has been made on this front. 12−16 Even the cheaper correlated methods, especially second-order perturbation theory, will consistently overbind systems with significant delocalization such as the benzene dimer, 9 while underbinding many others, for example, in saturated systems. 17 The gold standard is considered to be coupled cluster with singles, doubles, and perturbative triples [CCSD(T)], 18,19 which scales as N () 7 , where N is a measure of system size. Problematically, the correlation energy generally converges very slowly with basis set size, such that very large bases are needed to achieve high accuracy. 20 The result is that quantitative study is restricted to small molecular systems in the gas phase.
Considerable attention has been given to trying to reduce the cost of coupled cluster calculations. 21−24 The primary reason for slow convergence with basis set size is the poor description of the interelectronic cusp 25 when Gaussian-type orbitals are used. Recognition of this led to the development of explicitly correlated F12 methods, 26,27 where terms linear in the interelectronic separation are included as geminals in the wave function. Such methods dramatically improve the rate of convergence to the complete basis set limit, 28 allowing much smaller basis sets to be used without loss of accuracy. However, this comes at the cost of additional many-electron integrals; it was only with the advent of density fitting, 29,30 or resolution of the identity, approaches that it became feasible to carry out explicitly correlated calculations efficiently. Density fitting involves approximating the many-electron integral tensors by lower rank tensors using large auxiliary basis sets, a technique which has been found to be more generally applicable. 31,32 For example, it has been used in the Fock build 33,34 and integral transformation 35−37 stages of canonical wave function-based methods, yielding considerable savings in computational effort. 38−40 Unfortunately, in the case of coupled cluster, the use of density fitting does not improve the overall scaling (although it can considerably reduce the scaling prefactor). More recently, there have been attempts to use low-rank decompositions on the coupled cluster amplitude equations directly, with some promising early results. 41 The prohibitive scaling behavior of correlated wave function methods is unphysical. 42 It is an artifact of the use of canonical molecular orbitalsthe requirement for diagonalization results in orbitals extended over the entire system. As noted above, however, the properties of molecules are usually localized, and we would expect this to be reflected in the orbitals, at least in the occupied subspace. It is possible to use the invariance of the Hartree−Fock reference to rotations within the occupied and virtual subspaces to localize the molecular orbitals after the fact. 43 When used subsequently in the correlated part of the calculation, this has two advantages: the pair correlation between distant, localized orbitals is small and so can either be treated more cheaply or not at all; the steep and unphysical expansion of the virtual space available to each electron can be eliminated by restricting excitations to atomic orbitals spatially close to the occupied molecular orbitals. There are now a myriad of such "local correlation" methods, 23,44−49 which allow for coupled cluster calculations that scale essentially linearly with system size. None of the orbital localization procedures used in local correlation approaches are perfect, however, and great care needs to be taken in choosing thresholds for when terms are neglected; in particular, a lack of consistency across different geometries could result in discontinuous potential energy surfaces.
One method, related to coupled cluster, which has recently seen a lot of interest in the density-functional theory community is the random phase approximation (RPA). This method was first introduced in the 1950s by Bohm and Pines 50−52 as a way to determine exactly the correlation energy of the uniform electron gas. Previous perturbative attempts failed, as the perturbation series were divergent, while the renormalization inherent in the diagrammatic approach of Feynman, and later Goldstone, 53 circumvented these problems. Since then, it has been demonstrated that the RPA yields the exact long-range dispersive behavior, and it has been successfully applied to a wide range of molecular systems. 54,55 As we discuss later, it is also intimately connected with coupled cluster theory and as such has been used to bridge the gap between density functional-based and wave function-based approaches. 55 The problems in performing high-accuracy calculations are further compounded by the fact that in practically interesting systems, e.g., biological and supramolecular contexts, complexes are not normally found in their equilibrium geometries. The interaction of interest will usually be strongly influenced by the environment, such that the surroundings must be taken into account. In large systems with multiple components, new difficulties present themselves: (1) Many-body effects become important, meaning that pair potential approximations are inadequate. 56 (2) The number of local minima increases exponentially with system size. 57 (3) The intermolecular and intramolecular degrees of freedom become strongly coupled. 58 (4) The gradients become tainted with basis set superposition errors (BSSE). 59 The usual counterpoise scheme 60 for correcting BSSE would involve calculations on all fragments in the full complex basis, in addition to the supermolecular calculation itself. The cost therefore explodes as the number of fragments increases. In addition, BSSE is larger for smaller basis sets to the point where for very small bases it can be substantial enough as to result in a false minimum on a repulsive potential energy curve. 61 Therefore, to study intermolecular interactions in large systems, it is preferable to eliminate the superposition errors explicitly from the beginning.
Several different methods have been suggested for the a priori elimination of BSSE. These include symmetry-adapted perturbation theory (SAPT), 62−65 block-localized wave functions, 66 absolutely/extremely localized molecular orbitals, 67−69 local correlation methods, 42,49,58,70,71 dual basis set methods, 72−74 and the chemical Hamiltonian approach. 61,75,76 The first of these is perhaps the most popular and has the added advantage of providing a decomposition of the interaction energy into physical components, such as electrostatics and dispersion. However, the cost of high-order SAPT is not much better than traditional correlated methods, although recent density-fitted implementations are very efficient. 7,77 The chemical Hamiltonian approach, in contrast, attempts to directly identify the terms in the supermolecular Hamiltonian that introduce BSSE and removes them. 61,75,78 This has the advantage of being readily extendable to traditional correlated methods, in a way that SAPT cannot. The mean-field results obtained in this way have been seen to be virtually indistinguishable from counterpoise-corrected Hartree−Fock results, 75 without the need for monomer calculations in the full basis. However, it offers no computational savings beyond this and does not solve the underlying problem of localization.
An alternative method is to constrain the molecular orbitals (MOs) on a fragment to only use the atomic orbital (AO) basis local to that fragment. As the AOs are by construction localized, this results in absolutely localized molecular orbitals (ALMOs), with all BSSE necessarily eliminated as a result. This was first suggested by Stoll et al. 79 and then further developed by Cullen. 80 It has sometimes been attributed to the later work of Gianinetti et al. under the name self-consistent field for molecular interactions (SCF-MI); 81,82 this was then extended by Iwata and co-workers. 83−87 More recently, the Head-Gordon group have made extensive use of ALMOs, 69,88 introducing charge transfer corrections and an energy decomposition analysis. 89−92 Applications of the latter to a wide range of systems have yielded interesting results. 93,94 In the present work, we use absolutely localized molecular orbitals to develop an inherently fragment-based method for calculating interaction energies free from BSSE. These are used to eliminate the localization problems in a local treatment of the correlation energy and in this way naturally get a correlated energy decomposition analysis. Following development of the theory, we demonstrate how the new method can be made to be essentially linear scaling in the number of fragments and describe its amenability to massive parallelization. Finally, the applicability and efficiency of the method is demonstrated on a range of chemically interesting systems, showing that it gives essentially coupled cluster quality results at a fraction of the computational cost.

THEORY
The theory behind absolutely localized molecular orbitals has been developed, and redeveloped, by several groups. We present it here briefly in the context of second quantization, as this is helpful when including electron correlation later in the theory.
Consider a system of F fragments. Each fragment, X, has an atomic orbital (AO) basis, {χ Xμ }, associated with it, containing n X functions. From these, o X occupied and v X virtual molecular orbitals, ϕ Xp , will be formed. The total supermolecular basis then comprises the union of all such fragment bases, with N = ∑ X n X AOs, O = ∑ X o X occupied orbitals, and V = ∑ X v X virtual orbitals. Throughout the following, we define Tensor notation is used throughout, as is the Einstein summation convention with the exception of sums over fragments, which are always explicitly shown for clarity. As is usual, subscripted indices indicate covariant quantities, while superscripted indices represent contravariant quantities. Dots are used as placeholder indices to clarify the composite indices required to denote the fragmentation; for example, Xi is a single index.
2.1. Absolutely Localized Molecular Orbitals (ALMOs). The covariant ALMOs are expanded as a linear combination of atomic orbitals, with the coefficients constrained to be fragment localized. In tensor notation, this is written as where the coefficient tensors, T and V, are constrained such that This is in contrast to the canonical theory, where we would instead constrain all MOs to be orthogonal; the above removes the degrees of freedom that would allow us to do so, such that the orbitals so obtained are in general nonorthogonal. The wave function is more compactly represented in second quantization by operators, a pγ † and a pγ , which create or annihilate, respectively, an electron in orbital ϕ p with spin γ. The lack of orthogonality means that these satisfy the anticommutation relations where s pq = ⟨ϕ p |ϕ q ⟩. This completely encodes the antisymmetry of the wave function, but the nonvanishing nature of s pq complicates matters. Instead, we introduce a symmetrically orthogonalized set of operators for the occupied orbitals as where σ is the occupied−occupied overlap metric. The density operator, which projects onto the occupied subspace, is thus ρ= |ϕĩ⟩⟨ϕĩ|; we transform the virtual orbitals by projecting them out of the occupied subspace The occupied−occupied and occupied−virtual anticommutation relations thus become the usual Writing a single excitation operator in the form Ep q = ∑ γ ap γ † aq γ , the usual molecular Hamiltonian is given by pqrs pq rs ps rq (6) where hp q and gp qrs are the one-and two-electron integrals in the orthogonalized basis. For simplicity, we consider only a restricted closed-shell wave function, |Φ⟩ = ∏ i∈occ. a iα † a iβ † |⟩, but the extension to the spin-unrestricted case is simple, as has been detailed elsewhere. 95 The energy is thus given by where H and F are the AO representations of the core Hamiltonian and Fock operators, ĥand f, with F μν = H μν + ⟨μν∥λτ⟩P τλ . From this, we see that the density matrix is given by (8) Note that this couples fragments together, such that the energy is not simply a sum of disjoint monomer energies.
Taking the derivative of the energy with respect to a rotation, e −κ̂, of the orthogonalized orbitals, where κ̂= κ pq (Ẽp q − Ẽq p ), we get the stationary condition (9) This is equivalent to the requirement qfρ̂|ϕ̃X i ⟩ = 0. Thus, occupied−occupied and virtual−virtual rotations within each fragment's orbital space do not affect the energy, such that we may choose orbitals within each fragment to be orthonormal, diagonalizing the on-fragment blocks of the Fock matrix: ⟨ϕX j |f| ϕ̃X i ⟩ = ϵ Xi δ ij . We therefore define fragment-localized density and Fock operators as Note that ρ̂= ∑ X ρX, and that fX has been constructed to be Hermitian. Using the condition in eq 9, along with the facts that ρX |ϕ Yi ⟩ = δ XY |ϕ̃X i ⟩ and q̂|ϕ Xi ⟩ = 0, we see that Transferring to the AO representation, this yields a set of modified Roothaan−Hall equations for each fragment block

Journal of Chemical Theory and Computation
T h i si st h es a m er e s u l to b t a i n e db yS t o l le ta l . 79 By construction, there can be no BSSE present in the solution to these. Moreover, given that the diagonalization of an M × M matrix is an M () 3 process, the cost for diagonalizing F fragment Fock matrices is (Fn ) 3 ,a so p p o s e dt o = NF n () ( ) 33 3 for the full matrix. 2.2. Perturbative Correction to the ALMO SCF Energy. By restricting the orbitals to be absolutely localized, the energy obtained must, in general, be higher than the full HF solution. The missing density can be attributed to the lack of orbital rotations between fragments, i.e., the lack of density transfer between fragments. In fact, a Mulliken analysis 96 trivially shows that the charge on each fragment remains static under the ALMO approximation; we note, however, that this is not true using a Loẅdin analysis, 97 implying there is still a degree of uncertainty.
Equivalently, the difference can be thought of as the monomers each being restricted to a smaller, local AO basis instead of being allowed to relax in the full supermolecular basis. Dual-basis methods offer an efficient way to expand the converged solution into the larger basis. 98−101 In this approach, we diagonalize the Fock matrix and construct a new density, P ∞ , from which the energy correction is determined as However, diagonalization of the full Fock matrix is a costly N () 3 procedure. Instead, it can be circumvented by iteratively finding an orbital rotation, u, that diagonalizes the Fock operator, as outlined by Liang and Head-Gordon. 102 Defining , where the o and v subscripts denote occupied and virtual blocks, running over all fragments, the energy correction is simply given by (12) The rotation parameters are found by iterative solution of vv vo vo oo vo vo vo (13) This formulation is useful for two main reasons. First, its most expensive part computationally is the quadratic term, involving two matrix−matrix multiplications. It therefore scales an order of magnitude better than diagonalization and despite being iterative involves far fewer floating point operations. It should be noted that the quadratic term can be dispensed with, leading to the usual second-order peturbative correction, removing the need for iterations. However, the second advantage of the present approach is seen later in its utility for decomposing the energy correction into contributions between fragments.
2.3. Including Electron Correlation. As discussed in the Introduction, dispersion is often very important in intermolecular interactions. 7 Dispersion can be identified with dynamical electron correlation, of which only exchange terms are included in the mean-field treatment above. As such, it is vital to go beyond the mean-field limit if we are to treat supermolecular systems accurately. It is well known that dispersion can be formulated exactly in terms of the frequencydependent electric polarizabilities of a system, using a technique called the "adiabatic connection". 54,103−106 In simple terms, the Hamiltonian in eq 6 is decomposed into a zerothorder part (the Fock operator) and a fluctuation potential, ŵ, controlled by a parameter λ: Ĥ= f̂+ λŵ. By slowly, or adiabatically, "switching on" the perturbation from λ =0toλ = 1 and integrating, the correlation is recovered. If we make the assumption that only the part of the response of the electrons to a light wave that is in phase with the wave makes a nonnegligible contribution, we arrive at the random phase approximation (RPA). [50][51][52]107 The correlation energy is then found by integrating the fluctuation over the response in the electric polarizabilities, themselves determined by solution of a Dyson equation. This is the approach taken both in timedependent and range-separated density functional theory. 54,105,108−113 In this way, it can be shown that the correct R −6 long-range dispersive behavior in the interaction energy is recovered. 6 More recent interest in RPA has arisen from its intimate connection with coupled cluster methods. 114−119 In fact, direct RPA (dRPA), 103 where the fluctuation w iajb is simply the twoelectron integral (ia|jb), is almost equivalent to coupled cluster doubles (CCD) with only the diagrams with ring topologies included. 117 This is termed ring coupled cluster, or rCCD; there is a slight difference, however, in that the cluster amplitudes are not required to be symmetric with respect to orbital permutation. Starting from the orthogonalized basis of eq 4, we further orthogonalize the virtual orbitals within themselves in the same way where π is the virtual−virtual overlap metric. Note that the projection out of the occupied space necessarily delocalizes the virtual orbitals somewhat, but this can be shown to be in the sense that leaves them closest to the original, absolutely localized orbitals. 120 They are therefore expected to remain well localized. The double excitation operator T2 can then be written as tĩ j ab aã † ab † ajaĩ where tĩ j ab are the excitation amplitudes. Writing |Φ ij ab ⟩ = aã † ab † ajaĩ|Φ⟩, the CCD energy and amplitude equations are given by It can easily be shown 117 diagrammatically that the above amplitude equation with only ring diagrams retained leads to the dRPA residual =̃+̃̃+̃̃+̃̃̃= R K At tA tKt 0 (16) Here, K̃i ajb =(iã|jb), and A = ϵ̃+ K, where ϵĩ ajb = δ ij fã b − fĩ j δ ab . This is an algebraic Riccatti equation, for which many simple and efficient methods of solution exist, 121 allowing iterative determination of the RPA energy, given from eq 15 as However, we would like to make use of the inherent locality of the orbitals, achieving a local correlation method with the orbitals localized a priori. To do so, we back-transform eqs 16 and 17 using the relations in eqs 4 and 14 to get the following  (18) Equation 18 is written entirely in terms of the ALMOs, allowing us to only consider pairs of fragments and thus reach linear scaling due to the rapid decay of interactions with distance between fragments, as is discussed in the pairwise RPA approximation later. However, written as above, we can evaluate the amplitudes as t ik ac π cd σ kl , retaining all of the computational advantages of traditional RPA, minimizing the on-fragment scaling.

Journal of Chemical Theory and Computation
2.3.1. Including Exchange. Direct RPA performs very well in the long range but provides a much poorer description of electron correlation in the short range. 119,122,123 For example, the pair-correlation function in the high-density electron gas becomes negative at small separations. 55,123 This is due to the missing exchange terms in the Coulomb interaction, i.e., the fact that antisymmetrized integrals are not used. At short separations, exchange-repulsion dominates as the distributions of electrons of the same spin attempt to overlap. Several different schemes have been suggested to ameliorate this problem, 55,119,122 one of the simplest of which is to simply calculate the energy by contracting the amplitudes with the antisymmetrized integrals This is called the second-order screened exchange (SOSEX) correction. 124 Such a correction only affects the energy, however, whereas a full solution requires inclusion of exchange effects in the residuals. By replacing all instances of K with B in eq 16 (and the equivalent local forms), we arrive at the RPA plus exchange (RPAx) version. 119 The energy becomes Tr RPAx (19) where the extra factor of 1/2 is necessary to avoid double counting. SOSEX is equivalent to including all contractions of the amplitudes with the interaction that does not form "loops", i.e., where exactly one particle and one hole line have been exchanged, whereas RPAx also includes exchange in the Coulombic screening when determining the amplitudes. We note that the second-order contribution in RPAx equates to the MP2 energy.

Energy Decomposition Analysis.
The motivation for the choice of development of the method outlined above has been 2-fold: to achieve linear scaling of computation cost with respect to the number of fragments, as is discussed in the next section, and to yield a physically motivated energy decomposition. Building on previous work using ALMOs, 89,90 we describe four main components.
The first of these is the "frozen" energy, describing the electrostatic interaction between the unrelaxed, frozen monomer orbitals in their distorted geometries. After the full SCF calculations have been performed on the monomers in their interacting geometries, each fragment has a set of unperturbed occupied orbitals, T X,frz. . Inserting these in eq 8 gives the frozen density, which due to the presence of the inverse metric is, in general, not the superposition of noninteracting fragment densities, as the orbitals will not be orthogonal between fragments. Using this density in eq 7 gives the total energy of the unperturbed system; the difference between this and the sum of the monomer energies is then defined to be the frozen energy. This therefore represents the interactions between unperturbed fragment electron densities, represented by the occupied orbitals, which in turn corresponds to multipole−multipole interactions. A further term, the relaxation energy, could be included at this point as the difference between the monomer energies in their relaxed and interacting geometries.
The next two terms are polarization and charge transfer. During the ALMO SCF procedure, the frozen orbitals above are allowed to relax, corresponding to the distortion of the monomer densities due to the presence of other fragments, i.e., polarization. The polarization energy can thus be defined as the difference between the energy calculated with the frozen and full densities. However, as was noted earlier, the absolute localization constraint prevents transfer of electron density between fragments. Allowing the orbitals to relax in the full basis, as described by eq 12, can therefore be robustly defined as representing charge transfer. Charge transfer is by definition a directional quantity, however, and it is desirable to be able to determine the contribution from "charge" transferred between fragments X and Y. We do this by introducing Mulliken partition operators Xa noting that the sum over all fragments and orbitals of these gives the ALMO occupied (density) and virtual subspace projectors, respectively. The energy due to charge transfer between fragments X and Y is then found by inserting these into eq 12 This highlights the usefulness of taking the iterative approach of eq 13. We note that whereas the frozen and polarization energies are constructed to be free from BSSE, the charge transfer term necessarily reintroduces some error. However, due to the perturbative nature of the correction, this error is expected to be small, as is shown later. By transforming the density in the same way, we can also define (Mulliken) partial charges: However, we stress that this does not represent a physical transfer of electrons from one entity to another; it is rather the nonlocal response of the electron density on each fragment when brought into proximity with the other fragments. That is, we have split polarization into local and nonlocal portions; ΔE pol. reflects the redistribution of electron density constrained to a finite volume around the fragment, defined by the extent of the local atomic orbitals, while ΔE CT encompasses the remaining relaxation.
The dynamical electron correlation in this method is determined using RPA, which encompasses to some extent infinite-order many-body effects on the single particle level. The a priori local nature of the orbitals allows for considerable savings in computational cost but also for the explicit removal of BSSE from the correlation energy. Solving eqs 18 (or the RPAx equivalent) leads to excitation amplitudes in terms of orbitals strictly localized to specific fragments. As only double

Journal of Chemical Theory and Computation
Article excitations are included in rCCD, such excitations can involve at most four fragments (one for each index). Extending an idea first introduced by Schuẗz et al. in the context of local MP2, 58 we can classify these excitations as described in Figure 1.W e note that this approach is similar to the diatomics/triatomics in molecules approach of Head-Gordon and co-workers 125,126 but across molecules rather than atoms.
We identify five categories of excitation, based on what they physically represent. When all of the orbitals are located on a single fragment, not shown in Figure 1, a term is clearly intramolecular. As these are only slightly perturbed by the presence of the rest of the complex, it is reasonable to approximate that these do not form part of the interaction energy. Dispersion is taken to be the response of on-fragment excitations to simultaneous excitation within a separate fragment, i.e., the dispersive coupling between the separated fragment electron densities. This corresponds to class (a) in Figure 1, whereas class (b) is exchange-dispersion. This is again a number-conserving coupling of excitations involving exchange between fragments. This is not well separated from the dispersion term, as the reference wave function is fully antisymmetrized and so inherently includes exchange-repulsion in all terms. As such, this is not directly comparable with the same term in SAPT. The dispersive energy contributions are therefore taken from eq 18 as The two remaining classes in Figure 1 are ionic and BSSE, all of which involve excitations from X to Y, without a corresponding excitation into X. It should be stressed that no electrons are transferred; these are rather couplings between different states in the underlying Fock space. The term ionic is denoting that the excitations are not particle numberconserving within a fragment, representing another form of "charge transfer". In particular, it is supposed that diagrams (e) and (i), which are double excitations from the same fragment into distant virtual orbitals on separate fragments, represent an attempt to improve the description of intramolecular correlation using the extended basis. In principle, any terms involving excitations between fragments could contain some measure of superposition error, but the aforementioned terms are likely to be the main contributors due to the sharp distance dependence of the underlying physical interactions.
As the RPA is performed on the ALMO density directly no charge-transfer correction is made to the wave function, only the energythe ionic terms will likely contain not only some BSSE but also a significant amount of the charge transfer effects that would be encompassed at the mean-field level by ΔE CT . In addition, the virtual orbitals are projected out of the occupied subspace, as per eqs 5 and 14, such that they are not completely localized. As a result, the ionic and BSSE terms are entangled, albeit only to a small extent; coupled with the fact that superposition errors are known to be much larger at the correlated as opposed to mean-field level, we choose to define the charge-transfer contribution as in eq 12 and neglect all but the dispersion and exchange-dispersion terms in Figure 1.Itis important to note that in doing this we are removing some of the higher-order many-body effects from the dispersion energy. These divide into two categories: many-body excitations and many-body interactions between fragments. The former are still included through the renormalization inherent in RPA and are not reduced any more than they would be in, say, CCSD. The latter, however, are necessarily only included indirectly through the mean-field density upon which the dimer RPA calculations are performed. Initial tests on benzene trimers have shown that the three-body effects are not removed completely (unlike in MP2 calculations) but are reduced compared to CCSD. In principle, higher-order effects could be included by doing a tripletwise RPA and so on, and this will need careful investigation in the future, as such many-body interactions can be very important in extended systems.

IMPLEMENTATION
The main aim of the present work is to present an accurate method for noncovalent interactions that is efficient and scalable to large systems. As we have described the system in terms of fragments, there are two scaling regimes to consider: scaling with fragment size and with the number of fragments. The latter magnifies the former considerably and therefore is the larger source of computation cost; effectively, it is a measure of the entire system size. As such, we focus here on achieving linear scaling with regard to the number of fragments. There are three main expensive parts of any SCF procedure: evaluation (and storage) of the two-electron integrals, forming the Fock matrix, and diagonalizing the Fock matrix. The latter exhibits the correct scaling by virtue of the use of ALMOs, as was mentioned earlier. In the RPA portion of the calculation, the main difficulties are in the transformation of the integrals to the molecular orbital basis and the solution of the amplitude equations. If the full integrals are then required, they can be formed in the AO or MO bases as bb T or bb̃T, respectively. The evaluation of B clearly requires NM () 2 computations, where N and M are the number of functions in the orbital and auxiliary bases, respectively. The Cholesky decomposition is M () 3 but with a very small prefactor, while the transformation is (ON ) 2 for the inner step and (OVN) for the outer step. Overall the algorithm is therefore cubic scaling, as is the memory cost. While this is a considerable improvement on the quartic and quintic scalings of the conventional integral and transformation algorithms, it is still problematic for large systems.
However, such steep scaling with system size is unphysical. If we consider two one-particle densities separated by a distance R, represented by Gaussians of combined angular momentum l and exponents α p and α q , it can be shown 128 that the electron repulsion integral, V, between the two asymptotically follows the relation Hence, the worst case is for spherical Gaussians, where the integral decays at the same rate as the Coulombic potential itself. Nonetheless, the integral tensors should eventually become highly sparse, with each fragment having an effectively fixed domain size over which integrals do not fall below a given threshold.
One possibility for exploiting this would be to directly use conventional sparse matrix formats, such as the popular compressed row storage data structure. 129 However, this is not useful if elementwise access is needed, as this would require expensive binary searches. To solve this problem, we have devised an adapted sparse data structure and accompanying algorithm, which fully utilizes the underlying fragmented structure of the system. Consider a slice of the three-index DF tensor B along the ABS axis, P. The matrix can be represented as shown in Figure 2, where a typical sparsity pattern is also indicated, and the blocks are fragment by fragment. This matrix can be represented by two arrays with length N Blocks = F(F + 1)/2, after symmetry is taken into account. The first is an array of zero-sized matrices, while the second is a list of boolean values reflecting whether the corresponding integral block is zero. The approach can trivially be expanded to include the ABS blocks by "z-marching" along the index P indexing as shown within each P slicethen incrementing the ABS index and repeating. In this way, only the nonvanishing blocks need be stored, but these can still be rapidly accessed individually due to the redundant, zero-sized matrices. As such, sparsity is not accurately reflected within each block. However, the number of significant blocks per fragment will eventually reach a fixed size, such that the storage and evaluation of the integrals reaches linear scaling. The data redundancy is a small compromise that allows for a considerably faster Fock build routine, as is discussed shortly.
An estimate for the maximum value within a block, weighted by the initial density P, can be found using the Cauchy− Schwarz inequality 128 (25) where R (XY)Z is the distance from fragment Z to the midpoint of fragments X and Y. This requires an accurate initial density and evaluation of the diagonal two-electron integrals, g XYXY . The latter can very efficiently be computed and stored once at the beginning of the calculation, and the density is formed from the converged monomer densities, such that it should be a good initial guess. Using this coarse, fragment-blocked For example, all diagonal blocks will necessarily be "dark" and need to be included, whereas most off-diagonal blocks (such as blocks 3 and 5 on the top row) are negligible.

Journal of Chemical Theory and Computation
Article screening, the integral evaluation and storage becomes linear scaling in system size. However, the Cholesky decomposition of the two-index quantities is still necessary, and eventually for very large systems, this will dominate the cost.

Local Fock Build.
The two-electron integrals and density matrix are used in each iteration of the SCF procedure to form the Fock matrix. The core Hamiltonian part is inexpensive: it is the Coulomb and exchange portions that conventionally show quartic dependence on the number of basis functions. Of these, the latter is, in principle, much more complicated than the former. Within the density fitting framework, the exchange matrix can be written as The contraction with the density occurs across the two threeindex integral tensors, resulting in an apparent quintic cost. However, if we write the density in the symmetrically orthogonalized basis, we get which is cubic scaling. However, a cubic scaling term every iteration is still a problem. Exchange is inherently a short-range effect, such that we can use the inherent localization of the ALMOs to reduce the unphysical cost. A "local exchange" approach was first proposed by Polly and co-workers 33 and has more recently been improved by Koppl and Werner; 34 we use their method as a starting point. The basic idea is that each localized orbital only has significant overlap with other orbitals within a fixed extent. Therefore, for each molecular orbital, i, an orbital domain can be assigned, [i], comprising the atomic orbitals that will give nonvanishing contributions to the exchange. As the system size increases, the domain sizes eventually become fixed due to the localization, such that the calculation becomes linear scaling.
A clear disadvantage of this approach, in general, is that the results so obtained are likely to be highly sensitive to the choice of domains. This is where the fragment localization of the present method is particularly useful. Instead of considering the significance of individual orbitals, the contribution of a fragment as a whole can be considered. This will of course result in the inclusion of orbitals that are not necessary, but the well-separated nature of the fragments precludes the exclusion of important contributions, somewhat resolving the strong dependence on choice of domain.
We use the following criteria for the domains. The AOs of fragment X are included in [i] MO if they make a significant contribution to the density on i; that is, if ii MO (28) where ϵ MO is a user-defined threshold. Typically a value of ϵ MO ≈ 10 −6 suffices. It is then presumed that, almost by definition, if the density of i has a significant contribution from {Xμ}, then this will also be true of k. Thus, [i] MO ⊂ [i] AO . However, as the Coulomb force decays slowly with separation, the AO domain will, in general, be considerably larger than the MO domain. Let μ denote the exponent of the most diffuse basis function on X and assume that the extent of this functionthe cutoff beyond which its norm changes less than a given threshold, ϵ AO is much greater than 1. Then, for i located on fragment Y, we include all the atomic orbitals from fragment X if where R XY is the center-of-mass separation of the two fragments. Again, a reasonable threshold is ϵ AO =10 −6 . Finally, the ABS domains are found by a population analysis, as suggested by Polly et al. 33 For i on Y, the ABS of X is included in [i] ABS if the partial charge associated with charge transfer from i to X is greater than a threshold, ϵ ABS . However, for ALMOs, a Mulliken analysis gives zero charge transfer, such that a Loẅdin analysis must be used instead. The condition is thus In the original local exchange method, 33,34 the domain would be further extended by including all atoms connected to the atom of the current orbital. However, by assumption, there are no covalent bonds between fragments here, such that this is not necessary.
In general, the domains should need to change with each iteration. However, given a high quality initial density, such as is provided by the converged monomer densities, the domains can be fixed in the first iteration. An additional complication, though, is that the restriction of the fitting basis means that the fitting metric, G, is no longer valid in the restricted space. As such, the full metric needs to be restricted by retaining only the blocks included in the ABS domain, but this needs to be done separately for each occupied orbital. In practice, this is a minor technical issue that does not affect the efficiency of the algorithm.
At face value, the Coulomb contribution to the Fock matrix should be much simpler. Using density fitting, it is given by Jb bP QT Q (31) meaning it can be formed by two NM () 2 matrix−matrix products. It is tempting to try to reduce this cubic scaling in the same manner as for the exchange, but as has already been mentioned, the Coulomb interaction falls off much more slowly with distance. As such, prohibitively large domains would be needed. Moreover, the local exchange algorithm requires storage of the raw three-center integrals, rather than those symmetrized with the inverse fitting metric. Therefore, we instead form the Coulomb matrix as (μ X ν Y |P Z )B̃Z P , where BJ ournal of Chemical Theory and Computation Article = G −1 P XY (XY|Z). In this way, we need only loop over nonzero integral blocks, (XY|Z), something which is also true for eq 27 with the raw integrals. Thus, the integral blocking not only severely reduces the memory imprint but also more rapidly brings the computational cost of the Fock build to linear scaling. The only remaining nonlinear term is the multiplication by the inverse metric in the Coulomb part, which is quadratic. Note that, while G will become sparse in the macroscopic limit, this does not imply its inverse will be sparse; in fact, the opposite is usually true. 129 However, as this is a single matrix−matrix multiplication, it is a negligible cost.
3.3. Pairwise RPAd. The calculation of the RPA dispersion energy requires integrals transformed to the molecular orbital basis, which can be achieved using density fitting as described earlier. Note that the requirements on the auxiliary basis to minimize DF errors are different for the integral transformation compared to the Fock build, so in general, a separate ABS must be used. For direct RPA and RPA+SOSEX, we have that K = bb̃T, allowing us to avoid ever forming the full four-index integral tensor. Moreover, we can use this to rewrite eq 16 as ϵ̃̃+ε̃=− tt u u T (32) where u =( 1+t)b. This is not possible for RPAx due to the exchange terms in the residual. The most expensive computation in each iteration will thus be the evaluation of the quadratic term on the right-hand side, which is OVM () 22 for dRPA and OV () 33 for RPAx. Thus, these are formally fifthand sixth-order scalings with system size, although the number of occupied orbitals, O, is generally much smaller than N.
As several iterations are required, this scaling is prohibitive. However, the fragment-based localization again allows for reduction to linear scaling. In the energy decomposition, we partitioned the RPA excitations as per Figure 1, dispensing with all but the dispersive terms, where such terms are inherently limited to pairs of fragments as rCCD (and therefore RPA) is restricted to double excitations. This does not mean that only two-body effects are included, though, as RPA implicitly includes contributions to infinite order at the orbital level. However, it does allow us to adopt a natural local correlation approach, where the orbital domains are restricted by performing separate RPA calculations on pairs of fragments. If the ALMO orbitals were used directly, this would be exact, but the projection of the virtuals makes the pairwise approach an approximation. On the other hand, one only needs to project the virtuals out of the smaller occupied space for the pair, such that localization is better maintained. Moreover, as dispersion has an asymptotic R −6 dependence, only pairs within a given distance of each other (we measure the centerof-mass separation) need be included, such that the number of RPA calculations per fragment eventually becomes fixed, leading to linear scaling. We term this approach "pairwise RPAd", where the "d" denotes the restriction to dispersion and exchange-dispersion excitations.
3.4. Parallelization. The fragment-based partitioning of the monomer calculations, integrals, Fock build, diagonalizations, and pairwise RPA calculations make this method what is termed "pleasingly" parallelizable, meaning that it can be divided into tasks that are essentially independent of one another. The lack of significant communication required in all but the formation of the density ideally suits a distributed memory parallelism, and several portions of the method would benefit from the use of accelerators. The technical hurdles involved with such massive parallelism are beyond the scope of this paper. Instead, we consider a simple multithreaded approach that can later be combined with a larger scale parallelization.
First, all integrals can trivially be threaded by balancing over shell pairs, triplets, or quartets, depending on the integral class. This can be done separately within each fragment integral block, such that the blocks can themselves later be parallelized over. Assuming that the orbital domains are of roughly similar sizes on each fragmenta coarse but reasonable assumption each fragment block of the Fock build can be threaded by chunking up the largest domain (always the AO domain). This is particularly simple for the exchange portion, but for the Coulomb portion, the first half transformation of the densityfitted integrals must be completed on all threads before the second half transformation is performed, somewhat limiting the effectiveness of the multithreading. Instead for this portion, it would perhaps be better to multithread the matrix multiplications. In the core ALMO iterations, the formation of the density matrix and inverse metric necessitates shared memory, such that these parts can simply be threaded fragment-by-fragment. The individual fragment diagonalizations can also make use of standard multithreaded linear algebra routines, as can the contractions in the iterative solution of the RPA amplitude equations.
3.5. Summary. The scalings of the different components of the calculation with respect to both fragment and system size are given in Table 1. It shows that all of the major steps are essentially linear scaling, with the exception of the Coulomb part of the Fock matrix. This is, however, only formally quadratic due to an unavoidable matrix−matrix multiplication; properly accounting for sparsity would eventually (but slowly) lead to linearity. The density fitting, however, requires a cubic Cholesky decomposition, which will dominate for very large systems, despite the efficiency of standard algorithms. We are currently looking at how this can be circumvented. Finally, the on-fragment scaling is still dominated by the RPA calculation. In principle, this could be alleviated by further use of local correlation methods, but applying such an approach would require a very careful consideration of the additional errors it would introduce.

RESULTS
The ALMO+RPA method has been implemented in our inhouse electronic structure package, 130 which utilizes the LIBINT a Upper and lower case quantities refer to the system and fragment sizes, respectively.

Journal of Chemical Theory and Computation
Article package 131 of Valeev and co-workers for the molecular integrals, and the EIGEN library 132 for all linear algebra, including tensor manipulations. Additional DF-MP2, 35 CCSD, 133 CCSD(T), 133,134 and M06-2X 135 calculations were performed in the MOLPRO package, 136 while DF-SAPT 77,137 calculations were carried out using Psi4. 138 Unless noted otherwise, all calculations use the aug-cc-pVDZ (aVDZ) orbital basis 139−141 with the associated JKFit 127,142,143 and MP2Fit 144 density fitting sets for the Fock-and MO-basis integrals, respectively. All timing and memory benchmarks are reported for calculations on a single processor on the same compute node, with two 10-core Intel Xeon CPUs and 256 GB of RAM. Scalings with system size are reported for the RPAx variant of the method, as this is the most computationally complex variant, as discussed earlier. In the following, we apply the ALMO+RPA method to a variety of different systems, including benchmarking against the S66 database. 10,145 In this way, we demonstrate and assess the accuracy of the method, the errors associated with the various approximations, and the scaling with both fragment and system size. Finally, we compare the energy decomposition with that from symmetryadapted perturbation theory. 4.1. Scaling. The scaling with number of basis functions per fragment is demonstrated for the water dimer in Figure 3 by using the progression of orbitals basis sets: cc-pVnZ and aug-cc-pVnZ, with n = D,T,Q. 139,140 The ALMO+RPAxd curve has a scaling exponent of 2.7, or approximately cubic, as should be expected from Table 1. Moreover, this represents savings of several orders of magnitude over conventional counterpoise-corrected CCSD calculations, as can be seen from Figure 3e.
To study the scaling with overall system size, we consider two extreme cases: a linear chain of hydrogen fluoride molecules and a series of water clusters. The former corresponds to the best possible scenario, with a single hydrogen bond added per new fragment and with interactions along a single spatial axis. The clusters, on the other hand, are known to be particularly difficult, 146 as they are tightly bound, with multiple interactions per fragment. First, the efficacy of the fragment-blocked density fitting scheme is shown in Figure  4. In both cases, the memory requirements become linear, but this happens much more rapidly for the chain. The reduction in memory means that 108 hydrogen fluoride molecules (roughly 3500 and 12,500 orbital and auxiliary basis functions) or 50 water molecules (2000 and 7500 functions) can be considered entirely in-core with slightly over 16 GB of memory. The different rates at which these savings are reached are due to the different sparsities of the underlying integral tensors, The time spent on each portion of an ALMO+RPAxd calculation as a function of total system size is shown in Figure  5. From this, the linear scaling is clear to see, with the exception of the integrals. In fact, the integral routines themselves are linear. This is a significant contributor to the cheapness of the Coulomb and exchange terms; instead, it is the Cholesky decomposition of the fitting metric hidden in this term that is leading to unfavorable scaling. In fact, for the 108 hydrogen fluoride fragments, just the Cholesky routine takes up roughly 40% of the total computation time. For larger fragments, however, it will take longer for this to become the bottleneck. Nevertheless, it would be best to try to avoid this problem, and attempts have been made in this direction elsewhere. 147 We also note that the Coulomb term should not formally be linear but appears to be so here, most likely by virtue of how inexpensive it is overall.
The significant reduction in cost is due to a number of approximations, each of which has an associated error. Figure 6 reports the average error in the interaction energy per fragment due to the density and local and pairwise RPA approximations for the first few fluoride chains. This shows that the error due to both, when the thresholds stated earlier are used, is less than 0.1 kcal mol −1 in the total interaction energy and stays remarkably static as the system size increases. We do not observe any of the convergence problems noted elsewhere, 148,149 most likely due to the coarseness of our approximation. In particular, the change in the density due to using the local approximation has essentially no effect on the   Finally, the performance of the multithreading described earlier is shown for both a large chain and a large water cluster in Figure 7. This shows excellent speedups up to around four threads before starting to tail off. The effect is much larger in the water cluster, as this is dominated by the parallelized bits of code, in particular, the Fock build. The fluoride chain, on the other hand, is dominated by the Cholesky decomposition, which was not threaded here. Overall, using only eight threads, the calculation on a cluster of 32 water molecules can be carried out in under 3 min; when combined with a distributed, fragment-based parallelism, there is clearly potential for very large systems to be considered with minimal time requirements. Moreover, it is important to go beyond a single interaction energy; the method should also accurately reproduce equilibrium geometries and molecular properties. We demonstrate this by considering potential energy curves for four different bimolecular complexes, as a function of the intermolecular separation. These systems are the water dimer, a widely studied example of hydrogen bonding; FCl··· OH 2 , an unusually strong halogen bond; 150 the helium−neon complex, which is purely dispersively bound; and the complex of benzene and HCN, where the cyanide is perpendicular to the benzene ring. Geometries are given in the Supporting Information. Figure 8 compares the ALMO+RPAxd results, where the exchange (x) is SOSEX, termed RPA(SOSEX)d, with other double excitation theories, namely, counterpoise-corrected CCSD and MP2. In addition, the CCSD result without counterpoise correction is given, demonstrating that the BSSE has been successfully eliminated. In all cases, the new method is in excellent agreement with the coupled cluster results, while MP2 performs considerably worse. Reassuringly, the locations of the energy minima and the overall shapes of the curves are reproduced by the new method.
4.3. Performance on the S66 Benchmark Set. We have performed more extensive benchmarks on single-point calculations at equilibrium geometries by doing so for every complex in the S66 database of biologically relevant noncovalent interactions. 145 These are of particular interest because one of the most appealing uses of the new method is to investigate how solvent affects the interactions in such systems. Calculations were performed using a variety of different methods with the aVDZ basis sets, with geometries (see SI) optimized at the CCSD/aVDZ level. The database is split primarily into three different "types" of interactions: hydrogen bonded, dispersively bound, and other (usually a mixture of the two). It is also possible to instead split by the dominant component in the SAPT decomposition: electrostatics, dispersion, or an even mix of the two; we will compare the classifications from the new method's energy decomposition to those from SAPT in the next section.
In Figure 9, we see that the ALMO+RPA(SOSEX)d method gives far better agreement with the equivalent CCSD results than any of the other methods, in particular, in the case of the

Journal of Chemical Theory and Computation
Article dispersion-dominated complexes. The dRPA and RPAX (full exchange) variants show similar but slightly worse errors, as they underpredict and overpredict the exchange contribution, respectively; all raw results are given in the Supporting Information. It should be noted that as MP2 tends to overbind the electrostatically dominated systems it can often appear closer to CCSD(T) results than CCSD for small basis sets. However, this is a fortunate error correction, rather than a physical feature of MP2; it does not contain the physics described by a triples correction. Moreover, the improvement is not consistent across all systems, so cannot be relied upon. This is why we in the first instance have compared it to a CCSD benchmark.
The perturbative triples contributions of CCSD(T) are known to be very important in accurately describing the interaction energies, however. 145 Figure 10 shows the error distributions instead compared to the CCSD(T) result. In this case, ALMO+RPA performs worse than MP2 overall but better on the hydrogen-bonded complexes. The improvement in the MP2 result is largely due to the dispersion-dominated systems, where the triples contribution is particularly large, such that the significant overbinding of MP2 relative to CCSD is compensated for. This demonstrates the necessity for the inclusion of triples effects in future iterations of the present method; we are currently investigating various approaches for doing this. In principle, as we have a CCSD-like formalism, we should be able to include a (T)-like correction, at relatively little extra cost.
It is useful to try to disentangle the error due to the use of ALMOs from that of retaining only the ring coupled cluster diagrams. This can be done by comparing the ALMO energies to counterpoise-corrected Hartree−Fock energies. Without the charge transfer correction, the mean-absolute error across the whole set is 0.95 kcal mol −1 , reducing to 0.33 kcal mol −1 when the full, iterative correction in eq 12 is applied. The latter error stays fairly consistent across all classes of interactions at 0.25, 0.42, and 0.32 kcal mol −1 for the hydrogen bond, dispersion, and other categories, respectively. The average BSSE present in the HF and CCSD correlation energies is −0.65 and −1.39 kcal mol −1 , showing that this has been effectively eliminated.
There is a technical issue here with the use of the local exchange procedure, however. As noted by Koppl and Werner,34 while the density is fairly insensitive to the local approximation, the energy is often quite sensitive to the choice of domains. As such, they recommend performing a single, full exchange build from the converged density; the difference between this and the local result will here be termed the Xcorrection. The coarseness of the domain building, i.e., its fragment-based nature, means however that the energy is far less sensitive in most cases. In fact, for all but the dispersiontype systems, it accounts for less than 3% of the total interaction energy. However, unless the distance threshold in the domain selection is adjusted to reflect the longer separations in the dispersive complexes, the X-correction can account in magnitude for up to 50% of the total interaction, and on average 9%, and so must be included. The larger percentage is also partly due to the fact that these systems tend to have smaller interaction energies. This may become a problem when more fragments are present, as the full exchange build will scale cubically with system size; in these cases, it is then best to select a distance cutoff where at least nearestneighbor fragments are included in the domains for a fragment.
The average timings for DF-MP2, CCSD, CCSD(T), M06-2X, and DF-SAPT2 relative to ALMO+RPAxd both with and without the X-correction are given in Table 2. These show that even with the correction term the new method is at least twice as fast as density-fitted MP2 and over an order of magnitude faster than CCSD, despite giving very similar results to the latter. Moreover, these speedups will increase rapidly when either the number of fragments or size of the basis is increased. To consider convergence of the energies with respect to the latter, calculations were performed on the benchmark geometries from the original S66 paper 145 using the aVTZ basis sets; the mean-absolute error compared to the complete basis set limit CCSD results was 0.60 kcal mol −1 .
4.4. Comparison with SAPT. In the above, DF-SAPT2 was chosen for comparison as it is probably the most similar method currently available to the ALMO+RPA approach: (1)

Journal of Chemical Theory and Computation
Article It eliminates BSSE from the beginning. (2) It is designed specifically for interaction energies. (3) It contains a coupled cluster-style correlation energy. (4) It provides an energy decomposition analysis. The components from both the present method and SAPT for the entire S66 data set using the aVTZ basis are given in the Supporting Information. Here, we just compare how the new method classifies the different interactions and thus whether this agrees qualitatively with SAPT. The original classification came from the ratio of dispersion to electrostatic energy terms. 145 However, there is no direct comparison of the SAPT electrostatic term with, for example, the frozen energy term in the ALMO EDA, because the latter already includes exchange. Instead, we note that by calculating the SAPT charge transfer term of Stone and Misquitta 151 we can extract an equivalent to the ALMO polarization term by removing charge transfer from the SAPT induction. Then, considering the ratio, r D/P , of dispersion to polarization, we can recover exactly the original classifications: electrostatic, r D/P < 2; mixed, r D/P < 8; and dispersive, r D/P ≥ 8. A comparison is shown in Figure 11. In general, it can be seen that qualitative agreement between the two is good, with only four complexes resulting in different categories being assigned. Two of these are borderline cases where one method gives r D/P ≈ 7.9 and the other slightly over the cutoff (T-shaped benzene dimer and benzene···MeNH 2 complexes). The other two (MeNH 2 with peptide and pyridine) the ALMO EDA classifies as electrostatic (r D/P ≈ 1.6, 1.8), whereas SAPT classifies them as mixed (r D/P ≈ 3.0, 5.6). In both instances, this is due to SAPT predicting a larger dispersion contribution than ALMO+RPA. This is borne out in the total interaction energies, where SAPT and ALMO+RPA overbind and underbind by about 0.5 kcal mol −1 , respectively, compared to CCSD/CBS results.

CONCLUSIONS
We have presented a new method capable of accurately determining BSSE-free interaction energies for supramolecular systems composed of many fragments. By building on absolutely localized molecular orbitals and applying a variety of cutting-edge computational chemistry techniques, it has been possible to arrive at a procedure that gives CCSD quality results for a range of chemically interesting systems at a fraction of the cost of full, counterpoise-corrected calculations. In particular, a chain of 108 hydrogen fluoride molecules Figure 11. Ratio of dispersion to polarization contributions, r D/P , in the interaction energies of the complexes in the S66 database. Results are compared for the SAPT (red, left bar) and ALMO (blue, right bar) energy decompositions, with the cutoffs for classification shown as dashed lines. The molecule number is as given in the original database. 145 All numerical data can be found in the Supporting Information.

Journal of Chemical Theory and Computation
Article corresponding to over 1000 electrons, 3500 orbital basis functions, and 12,500 auxiliary basis functionscould be treated on a single processor in 486 s, or slightly over 8 min, using less than 10 gigabytes of memory. It was found that by specializing to systems composed of a large number of smallto medium-sized monomers, it is possible to make full use of the inherent localization of the ALMOs both in the mean-field and correlated regimes. In particular, problems with large energy errors in the local Fock build are alleviated by using a coarser definition of the orbital domains, and large superposition errors are mostly removed from the correlation energy by decomposing the excitations into distinct terms. Finally, the new method allows the interaction energy to be decomposed into physically distinct terms in a quantitatively accurate way. This coupled with the essentially linear-scaling and inherently parallel nature of the calculation opens up the possibility of gaining insight into the interactions in extended systems, such as polymers, supramolecular complexes, and even crystals.
There are caveats, however, mainly in the scaling with individual fragment size, which is prohibited by the use of RPA. Possible ways to improve this could include developing method-specific basis sets to minimize the problem or translating the problem to a density-functional theory setting, where RPA correlation is routinely used. This would, however, remove the ability to decompose the excitations into different terms. In addition, the Cholesky decomposition required for the density fitting procedure becomes a bottleneck for very large systems. Again, this could be alleviated by utilizing massive parallelization. In this way, it would become possible to routinely treat systems involving many hundreds of molecules in an accurate manner and to gain physical insight into the main underlying physical effects through energy decomposition. The ability to partition this further into interactions between specific pairs of molecules will be particularly useful in studying the effects of cooperativity between multiple noncovalent interactions, which is especially important in the condensed phases. Finally, we reiterate that the inclusion of higher-order excitations, in particular triples, is vital to the accurate description of intermolecular interactions. The rCCD-based approach here is set up in such a way that such effects can readily be included, for example, as is done in the DLPNO−CCSD(T) method, 152 and we are currently working on implementing this.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.9b00615.
Interaction energies and energy decompositions, both tabulated and in raw formats; errors and calibration graphs for local thresholds; and geometries of the complexes in Figure 8.( PDF) Cartesian coordinates of the CCSD/aVDZ optimized geometries for the S66 set. (ZIP)