One of the several problems that plague majority of density functional theory calculations is their inability to properly account for long-range correlations giving rise to dispersion forces. The recently proposed many-pair expansion (MPE) [T. Zhu *et al.*, Phys. Rev. B **93**, 201108(R) (2016)] is a hierarchy of approximations that systematically corrects any deficiencies of an approximate functional to finally converge to the exact energy. This is achieved by decomposing the total density into a sum of two-electron densities and accounting for successive two-, four-, six-,… electron interactions. Here, we show that already low orders of MPE expansion recover the dispersion energy accurately. To this end, we employ the Pariser-Parr-Pople Hamiltonian and study the behavior of long-range interactions in *trans*-polyacetylene as well as stacks of ethylene and benzene molecules. We also show how convergence of the expansion is affected by electron conjugation and the choice of the density partitioning.

## I. INTRODUCTION

Density functional theory (DFT) in the Kohn-Sham (KS) scheme appears today as the most prominent electronic structure method used throughout various disciplines.^{1–3} This results from an unprecedented combination of computational efficiency and accuracy. Despite its undeniable success, the need for an approximate exchange-correlation functional puts serious limitations on the reliability of the so obtained results. While there exists a hierarchy of approximate functionals that puts them on different rungs depending on their ingredients,^{4} their performance for a particular class of systems or properties is often hard to assess without careful benchmarking.^{5–8} There are several deficiencies that pertain to the vast majority of exchange-correlation functionals used in practice and are the root of systematic errors in DFT calculations. The major reason for many failures of the density functional approximation is the self-interaction error (SIE), which results from a spurious interaction of an electron with itself.^{9,10} The SIE is exposed most obviously in calculations of reaction barriers, charge-transfer, and ionization processes. Another systematic failure of common density functionals is their inability to describe strong electron correlations, which occur when the independent-particle picture breaks down. This often happens for molecules and materials containing transition metals with incompletely filled *d* or *f* shells or for bond-breaking processes.^{11,12} Finally, the semi-local correlation functionals are intrinsically not capable of accounting for long-range correlations, i.e., dispersion interactions.^{13,14} As dispersion is often a substantial or even dominating term of interaction energy, this puts severe limitations on applicability of DFT to intermolecular interactions or layered materials.

Since the foundations of DFT have been laid, much of the progress is motivated by efforts to overcome these systematic errors of approximate functionals. The problem of self-interaction has been approached from two different angles. One is to apply a self-interaction correction (SIC) to a (semi-)local exchange functional that removes spurious one-electron self-interaction.^{9} The other approach to combat SIE is to use exact exchange which is one-electron self-interaction-free by construction. The challenge is to design a matching correlation functional, which would not only include the dynamic part but also account for static correlation.^{15–18} Attempts to include strong correlations in DFT have led to the development of such approaches as adding a Hubbard-like repulsion term (DFT+U),^{19} combining DFT with dynamic mean field theory,^{20} using the strong interaction limit of the Hohenberg functional^{21} or inverting the KS potential from accurate density matrix renormalization group calculations.^{22} The problem of missing dispersion interactions has also triggered the development of various methods to include them in DFT calculations. The conceptually simplest approach is to add a dispersion energy correction based on the known asymptotic limit of dispersion energy, in which the leading term behaves as −*C*_{6}/*R*^{6}, where *R* is the distance between interacting atoms and *C*_{6} coefficients are either tabulated or computed from density.^{23–30} A less empirical route to dispersion in DFT are van der Waals correlation functionals, which model dispersion energy through a non-local functional of the density.^{31–35} Finally, functionals using unoccupied Kohn-Sham orbitals are also capable of accounting for long-range correlations. They either take the form of double hybrids,^{36,37} which add a fraction of correlation energy calculated from perturbation theory, or make use of the adiabatic connection fluctuation-dissipation theorem^{38–40} to calculate the correlation functional.

As briefly overviewed, there are multiple ways to correct for a particular deficiency of an approximate exchange-correlation functional. Unfortunately, none of these is able to cure all problems at the same time. Different correction schemes often rely on empirical parameters that have been fitted to a particular data set and their reliability is difficult to predict. Also, improving one property may lead to the deterioration of another as it is often the case for self-interaction corrections.^{41,42} In other words, there is no established systematic and practical way to improve DFT results that starts from any approximation and performs uniformly better for the whole spectrum of properties and systems. To address this problem, we have recently proposed the many-electron expansion (MEE),^{43} a hierarchy of density functional approximations, which converges to the true ground-state functional. Applications to the Hubbard^{44} and Hubbard-Peierls^{45} Hamiltonians have proved its usefulness in dealing with strongly correlated electrons and the self-interaction error. The aim of this work is to elucidate how MEE deals with long-range Coulomb interactions. In particular, we are interested in answering the question how long-range correlations (dispersion) are recovered by successful MEE corrections. To this end, we use the Pariser-Parr-Pople (PPP)^{46–48} lattice Hamiltonians for polyacetylene as well as stacked ethylene and benzene molecules. We show that MEE corrections decay rapidly with distance and low orders of expansion are sufficient to obtain accurate results.

## II. THEORY

### A. Many-pair expansion (MPE)

The working equations of MEE were derived in Ref. 43. For a self-consistent presentation, we repeat the final results also in this work. As we deal with spin compensated systems only, we will use the many-pair expansion (MPE) formalism, which is based on the decomposition of the total electron density $ \rho ( \mathbf{r} ) $ into a sum of pair densities $ { \rho i ( \mathbf{r} ) } $ such that

where *N* is the number of electron pairs in the system. Note that by a pair density we mean one-particle density that integrates to two electrons and not a two-particle density $ \rho ( \mathbf{r} , \mathbf{r} \u2032 ) $. The idea of MPE is to start with an approximate density functional $ E a [ \rho ] $ and systematically improve its errors by accounting for corrections calculated for a few-electron densities at a time. For any *ρ*, the correction is defined as $ \Delta E [ \rho ] \u2261 E v [ \rho ] \u2212 E a [ \rho ] $, where $ E v [ \rho ] $ is the exact energy functional. Then we consider the following hierarchy of approximations to the true energy of the 2*N*-electron system:

In this expansion, only $ E 0 [ \rho ] $ is an explicit functional of the total density $ \rho ( \mathbf{r} ) $, while all higher order terms are functionals of the pair densities $ { \rho i ( \mathbf{r} ) } $ defined in Eq. (1). Since there are infinitely many ways to decompose the total density, MPE becomes an implicit functional of the total density once the decomposition is prescribed. This is analogous to orbital-dependent functionals,^{49} which are only implicit functionals of the density, through explicit dependence on orbitals defined by Kohn-Sham equations. In general, MPE energies *E*_{i} will depend on the particular choice of partitioning, except for the *N*th order energy *E*_{N}, which recovers the exact energy for a 2*N*-electron system. This means that MPE is a finite expansion that converges to the exact result for any functional chosen for the zeroth-order approximation. Here, by the exact result we mean the exact energy of the fully interacting system with the ground-state density $ \rho ( \mathbf{r} ) = \u2211 i \rho i ( \mathbf{r} ) $. This is equivalent to solving the exact diagonalization (full configuration interaction) problem constrained to give $ \rho ( \mathbf{r} ) $ as the ground-state solution.

For MPE to be a practical method we need to be able to compute approximate and exact ground-state energies (*E*_{a} and *E _{υ}*) for all pair densities as well as their sums, which are combined to form many-pair subsystem densities ($ \rho i + \rho j + \cdots $). We will assume that the approximate energy functional is defined within the Kohn-Sham scheme, which makes it an orbital-dependent functional. Therefore, to get the approximate energy, we do a non-interacting potential inversion

^{50–52}to find a ground-state Kohn-Sham (KS) determinant that gives the input density $ \rho q ( \mathbf{r} ) $. Once we have the orbitals, they can be plugged into the approximate functional $ E a [ { \varphi i } q ] \u2261 E a [ \rho q ] $, which now has an explicit dependence on the orbitals. Similarly, by doing interacting inversion,

^{53}we can find an interacting wavefunction that yields the same density $ \rho q ( \mathbf{r} ) $. Calculating the expectation value of the fully interacting Hamiltonian in this state gives the exact energy of the system $ E v [ \rho q ] $.

To solve a noninteracting potential inversion problem for a given density $ \rho q ( \mathbf{r} ) $, we search for a local potential υ_{s}(**r**) such that

i.e., the input density is the ground-state density of a non-interacting system in the sought external potential. The wavefunction of this system is a Slater determinant constructed out of orbitals $ \varphi k ( \mathbf{r} ) $. To find this wavefunction, we search for a stationary point of the Lagrangian

where $v$_{s}(**r**) is the Lagrange multiplier enforcing the density constraint. At the solution each orbital satisfies a one-electron Schrödinger equation

and the density can then be expressed as

Since the potential *υ*_{s}(**r**) that yields $ \rho q ( \mathbf{r} ) $ through Eqs. (5) and (6) is not known in advance, we need to find it numerically to satisfy Eq. (3). To this end, we start with an initial guess $ v s 0 ( \mathbf{r} ) $ and use Newton’s method^{54} to find an improved potential

where $ \delta \rho ( \mathbf{r} ) / \delta v s ( \mathbf{r} \u2032 ) $ denotes the response kernel, which can be computed analytically from perturbation theory

Newton’s search process in Eq. (7) is iterated until the procedure is converged. Solving the Kohn-Sham potential (Eq. (5)) with the converged potential $ v s q ( \mathbf{r} ) $ allows to evaluate $ E a [ \rho q ] $ through its explicit dependence on Kohn-Sham orbitals $ E a [ { \varphi i } q ] $.

For an interacting inversion problem, a similar approach can be applied except that now we are searching for the exact wavefunction of a fully interacting system that has exactly the same ground-state density $ \rho q ( \mathbf{r} ) $. This wavefunction is a stationary point of the following Lagrangian:

where *υ*_{ex}(**r**) enforces the density constraint. Now the interacting wavefunction $ \Psi ( \mathbf{r} ) $ is the ground state of the interacting Schrödinger equation

and the corresponding density is

The search for the local potential *υ*_{ex}(**r**) is carried out analogously to the non-interacting case except for the fact that we use numerical finite differences to compute the response kernel. Once the converged potential $ v e x q ( \mathbf{r} ) $ is found, the exact energy is calculated as the expectation value of the interacting Hamiltonian $ E v [ \rho q ] = \u27e8 \Psi | H ^ | \Psi \u27e9 $, where $\Psi $ is the ground-state solution of Eq. (10).

Both $ E a [ \rho q ] $ and $ E v [ \rho q ] $ are ground-state energies of systems that have $ \rho q $ as the ground-state density. Since in the inversion procedures we are searching over local potentials only, $ \rho q $ has to be both interacting and non-interacting *υ*-representable. Formally, this condition needs to be enforced by the density partitioning. At present it is not clear how to construct such a scheme since the sufficient conditions for *υ*-representability are not known.^{55} It is possible that a constrained-search extension of the formalism would side-step this issue at the formal level. Nevertheless, this is unlikely to be a problem in practice provided that pair densities are sufficiently smooth. Note that even non-v-representable densities can be approached arbitrarily close with smooth potentials^{56,57} and that within the finite basis set approximation even node-containing densities can be effectively v-representable.^{58} Since a practical procedure would enforce the density match only up to some numerical precision, it would be sufficient to converge a v-representable density that is sufficiently close to that resulting from partitioning.

### B. Relation to the Perdew-Zunger self-interaction correction (PZ-SIC)

The MPE (or MEE) hierarchy systematically removes any possible errors of an approximate density functional $ E a [ \rho ] $. In particular, it can be seen as a method to systematically remove the self-interaction error, which is one of the major sources of density functional approximation failures. For this reason, it is instructive to establish a connection between MPE and the Perdew-Zunger self-interaction correction (PZ-SIC).^{9} PZ-SIC is exact for any one-electron system, while MPE converges to the exact energy for an arbitrary number of electrons; therefore, the latter can be viewed as a generalization of the former. Nevertheless, several important differences need to be noted.

For an N-electron system, the PZ-SIC can be written in the following form:

where *E*_{a} is an approximate functional of the total energy, *E*_{EXX} is the functional using the exact exchange for the exchange-correlation part, and the sum runs over a set of one-electron densities. For a one-electron system (*N* = 1), PZ-SIC is equivalent to the first order of MEE, because in the absence of electron-electron interactions *E*_{EXX} is the exact functional. In a many-electron case, Eq. (12) and MEE1 are still formally equivalent, but they are defined on a different set of admissible one-electron densities $ { \rho i } $. In PZ-SIC $ \rho i = | \varphi i | 2 $ are orbital densities, where $ { \varphi i } $ are occupied Kohn-Sham orbitals or localized orbitals obtained through a unitary transformation. As orbitals are in general excited states of a one-electron Schrödinger equation, they exhibit pronounced nodal surfaces and, as a result, the corresponding densities are not smooth. Considering that $ E a [ \rho ] $ is an approximation to the ground-state functional, its performance for such excited-like, non-smooth densities is likely to be much worse than for the total density, which naturally is a ground state of some potential. This caveat has been linked to the often unsatisfactory performance of PZ-SIC.^{41} The complicated nodal structure of higher-lying orbitals is a consequence of their mutual orthogonality. Recently, it has been proposed to use complex orbitals^{59–61} in PZ-SIC as the nodes of the real and imaginary parts do not need to coincide, so the resulting densities are considerably smoother. A similar problem with orbital nodes appears for functionals using the ratio of von Weizsäcker and exact non-interacting kinetic energy densities to detect one-electron regions of the density.^{62,63} Recently, a density-based alternative to this detector has been shown to be free of the nodal problem.^{64–66} In MEE, one-electron densities $ { \rho i } $ are required to be v-representable (be ground states of some potential). In practice, this means that orbital densities are not allowed and that admissible densities should be smooth. This requirement makes the energies $ E a [ \rho i ] $ compatible with the approximate total energy $ E a [ \rho ] $, as in both cases a ground-state functional is applied to densities that are known to be ground states. V-representability of one-electron densities needs to be assured by the initial decomposition, which is not uniquely defined. This is a similar situation as for PZ-SIC, which is not invariant with respect to orbital rotations. In both cases, the partitioning of the total density may be optimized variationally.

Another possible reason why PZ-SIC can deteriorate some results is that the SIE often mimics to some extent static correlation,^{67} so removing it destroys the helpful compensation of errors in DFT. PZ-SIC is designed to remove only the one-electron SIE, which is easy to define for every orbital separately and can be eliminated by the use of the one-electron self-interaction-free exact exchange functional. Unfortunately, as PZ-SIC does not add any correlation beyond the initial approximation, it can destroy the balance between exchange and correlation functionals in modeling of the total exchange-correlation hole. Also, PZ-SIC has no means to correct for the many-electron self-interaction error,^{68} which is also pervasive in approximate density functionals. In contrast, MPE (MEE) is able to systematically remove the many-electron SIE and include the corresponding degree of many-electron correlation. In the MPE formalism pair densities are treated at the first order, so MPE1 in addition to correcting for the one-electron self-interaction corrects the correlation energy of two electrons within each pair. This feature makes MPE1 conceptually similar to wavefunction-based pair theories like the antisymmetrized product of strongly orthogonal geminals (APSG).^{69–71} Both theories can account for static correlation and at the same time include some amount of dynamic correlation. They would also have a similar problem where the division into electron pairs needs to be reorganized during the optimization. However, MPE provides a hierarchy of improvements, which APSG lacks, and avoids a complicated optimization problem in terms of partitioning to orbital sets. At higher levels of MPE, more electrons are correlated at a time, so many-electron self-interaction and many-electron correlations are treated in a balanced way.

### C. Pariser-Parr-Pople (PPP) model

The equations presented up to this point are general. The actual implementation depends on the type of the Hamiltonian and representation of the density. While work on the implementation of MPE for *ab initio* Hamiltonians is ongoing, here we focus on lattice models for unsaturated hydrocarbons. To model the electronic structure of *π* electrons, we employ the Pariser-Parr-Pople (PPP) model, defined by the following Hamiltonian:

This can be viewed as an extension of the Hubbard model, which in addition to on-site repulsion characterized by the parameter *U*, accounts for long-range electron-electron interactions determined by the inter-site potential *V*_{ij}. A common choice for *V*_{ij} is the Ohno potential,^{72} which is an interpolation between the on-site interaction U and asymptotic 1/*r* dependence of the Coulomb potential. The Ohno potential has the following form:

where *r*_{ij} is the distance between sites *i* and *j*, and *e* is the elementary charge. The parameters *t*_{ij} and *U* are usually fitted to reproduce experimental data. For two carbon atoms separated by $ 1.4 \xc5 $, the standard PPP parameters are *U* = 11.26 eV and *t*_{ij} = −2.4 eV. Further, we assume exponential dependence of the transfer integral on the bond distortion $ \varphi = ( 1.4 \u2212 r i j ) $, i.e., $ t i j = \u2212 2.4 e \alpha \varphi $ eV, for covalently bound sites. The parameter $ \alpha = 3.785 $ was fitted to reproduce the transfer integral *t* = −2.9 eV for ethylene ($ r = 1.35 \xc5 $).

Application of DFT to lattice Hamiltonians is somewhat ambiguous as there is no real-space density $ \rho ( \mathbf{r} ) $ and there is some freedom in the choice of the basic variable.^{73–75} We use a formulation where site occupations, i.e., the diagonal of the density matrix, are employed for this purpose. Therefore, the density operator is simply $ \rho i ^ = a ^ i \u2020 a ^ i $. Further in this work, we consider only hydrocarbons whose all site occupations are equivalent by symmetry. Additionally, each sp^{2} carbon atom contributes only one electron to the $\pi $ framework, so we will consider only half-filled PPP models. This means that for each site *k* the ground-state density is just $ \rho k = 1 $. We impose these constraints to avoid the need for optimization of the density. This restriction is not a limitation of MPE itself but makes the initial study simpler and faster without any substantial loss of generality.

The crucial step of MPE calculations is the decomposition into pair densities according to Eq. (1). As the way to carry out this procedure is arbitrary, it is clear that MPE energies are not invariant with respect to different partitioning except for MPE0, which depends on the total density and MPE*N*, which is exact for a 2*N*-electron spin compensated density. For a quick convergence of the expansion, the choice of pair densities has to be made in such a way that subsequent orders of MPE correct for gradually weaker correlations. In practice, it means that we want pair densities to be localized in space and possibly well separated from the others. For the models considered, the partitioning is particularly straightforward as it is always possible to choose a pattern where each pair occupies two neighboring sites and all pairs do not overlap with each other. For systems with alternating bond lengths, we allocate a pair to sites connected by a shorter bond, while longer bonds connect sites occupied by two different pairs (Fig. 1). In a previous work,^{43} we have shown that such decomposition leads to a quick convergence of MPE energies for the Hubbard-Peierls model. This also corresponds well with intuition which tells us that electron pairs would localize on double (short) carbon-carbon bonds rather than on single (long) ones.

When doing the potential inversion calculations (Eqs. (3)–(11)) in the site basis, the density and potential are represented as vectors **𝝆** and **𝒗**, while the coordinates **r** and r′ become the lattice site indices. As a result, $ \delta \rho ( \mathbf{r} ) / \delta v ( \mathbf{r} \u2032 ) $ is replaced by a Jacobian matrix $ d \bm{\rho} / d \bm{v} $. We start the calculations by building the noninteracting and interacting *N*_{s}-site PPP Hamiltonians with the initial guess for potential $ \bm{v} 0 ( 1,2 , \u2026 , N s ) $ such that for site $\alpha $, $ \bm{v} 0 ( \alpha ) = \u221e $ if $ \bm{\rho} q ( \alpha ) = 0 $. By diagonalizing the Hamiltonian $ \mathbf{H} [ \bm{v} 0 ] $, we obtain the eigenvalues $ { \u03f5 i } $ and eigenvectors $ { \varphi i } $ and the corresponding density $ \bm{\rho} [ \bm{v} 0 ] $. As the Jacobian $ d \bm{\rho} / d \bm{v} $ may be ill-conditioned or singular, to compute Newton steps in Eq. (7) we use the Tikhonov regularization^{54}

where $ \sigma i $ are singular values of the Jacobian and $ \bm{U} , \bm{V} $ are unitary matrices obtained through singular value decomposition (SVD).^{76,77} We set the regularization parameter $ \lambda = 10 \u2212 6 $.

Finally, for the approximate energy functional $ E a [ \rho ] $ we choose exact exchange (EXX) and to compute $ E v [ \rho ] $ we use the exact diagonalization (full configuration interaction) method.^{78,79} With this choice, the zeroth-order energy is one-electron self-interaction free. However, it is completely missing any Coulomb-like electron-electron correlation effects and is affected by the many-electron self-interaction error.

## III. RESULTS

### A. Rapid decay of many-pair corrections

Due to locality of interactions in the Hubbard model, MPE has *O*(*N*_{s}) scaling at any order, where *N*_{s} is the size of the system. In the PPP model, the presence of long-range interactions requires the calculation of energy corrections for every possible combination of pair densities. This immediately leads to $ O ( N s M ) $, where *M* is the order at which MPE is truncated. However, we can expect that energy corrections will decay rapidly with the distance between electron pairs that are correlated. If we could neglect contributions from combinations of pair densities separated by a distance greater than some threshold, MPE would still scale linearly with the size.

To check how rapidly MPE correction vanishes with the distance between pair densities, we first consider a 100-site model for *trans*-polyacetylene with imposed cyclic boundary conditions. The bond lengths are $ 1.36 \xc5 $ and $ 1.44 \xc5 $ for the double and single bond, respectively. We consider four types of pair density patterns (see Fig. 2) that have one gap in the partitioning. For MPE2, one pair density is fixed and the other one is moved along the chain (MPE2 (1+1)). Similarly, for MPE3 (2+1) and MPE4 (3+1), two and three adjacent pair densities are fixed, and the remaining one is translated by two sites at a time. Additionally, for MPE4 we consider a (2+2) pattern, where two adjacent pair densities are fixed and the remaining two are translated as contiguous blocks. Fig. 3 shows decay of the magnitude of MPE2-MPE4 energy corrections as a function of separation between localized pair densities. Energy corrections for contiguous blocks of pair densities (shift = 1) are suppressed as they are two orders of magnitude larger than that for a pattern with a gap between the pairs. It is evident that corrections decay rapidly and fall below $ 10 \u2212 6 $ eV beyond the separation of 5 localized pairs (10 sites). More insight can be gained from plotting the same data on a double logarithmic scale (Fig. 4), this time including also the corrections for contiguous blocks. For shift = 1, the corrections are $ \u2212 0.064 $ eV, $ \u2212 0.011 $ eV, and $ \u2212 0.0039 $ eV for MPE2, MPE3, and MPE4, respectively. This already shows that, if pairs are localized in space, corrections decay quickly with the level of the MPE, each being a few times smaller than the preceding one. Furthermore, the magnitudes of corrections fall off rapidly also with the separation between localized pairs. The decay is polynomial as the dependence on the separation is approximately a straight line on a double logarithmic scale. The blue broken line in Fig. 4 is tangent to the MPE2 curve at shift = 10 and it follows closely the MPE2 curve for shifts greater than 4. The slope is −5.96, which means that MPE2 corrections decay as 1/*R*^{6} in the asymptotic range, which corresponds to the fact that they account for unscreened long-range correlations due to dispersion interactions.

The rapid decay of MPE corrections, both with the level of the expansion and the distance between pairs being correlated, is very promising for the applications of MPE. First, we can truncate the expansion at low orders as there is no need to correlate more than a few electrons at a time. This behavior has been already shown for Hubbard and Hubbard-Peierls Hamiltonians.^{43} Second, we can suppress the corrections for pairs separated by large distances, which would reduce the polynomial $ O ( N s M ) $ scaling at the cost of some small loss in accuracy. Both these features hinge on the decomposition of the total density to well localized pair densities, which is an essential component of MPE calculations.

### B. Dispersion interactions

Dispersion forces result from the interaction of instantaneous dipole moments, which in turn are generated by electron correlations. This is a weak but long-range effect, which lies at the heart of supramolecular chemistry and is decisive for the existence and structure of many systems of biological or technological importance. Unfortunately, being an effect of electron correlations, dispersion interactions are completely absent in the Hartree-Fock theory. Also, while the exact exchange-correlation functional takes account for these interactions, common semi-local approximations to the correlation functional are intrinsically unable to capture this effect due to its fundamentally non-local character.

To study how MPE accounts for dispersion interactions, we first start with a model for the ethylene dimer. Dispersion interactions in this system are dominated by correlations of *π* electrons, which are the most polarizable. In reality, different orientations of the dimer strongly affect the interaction energy.^{80} This includes also the dispersion contribution due to the directional character of the *π* cloud. The anisotropy cannot be accounted for in the model Hamiltonian that we use, which is represented in the basis of four sites in the positions of carbon atoms only. Therefore, we stipulate that our model represents a configuration where two ethylene molecules are stacked perfectly on top of each other and interact mostly through *π*–*π* stacking.

Due to the use of the site basis in the PPP model, there is no overlap between basis functions on non-bonded sites. Since the short-range Pauli repulsion depends on the orbital overlap, this effect is missing in our calculations, even though the EXX approximation can account for it explicitly. Allowing electrons on both molecules to correlate generates a force that is attractive at all separations. As a result, MPE interaction energies do not behave like in real $\pi $–$\pi $ interacting complexes. In particular, they cannot predict the existence of an equilibrium distance. To account for this effect, we add a pairwise correction to the energy that depends exponentially on distance between interacting sites

where $ \u03f5 = 377.2 eV $, $ D = 0.3455 \xc5 $, and *r*_{ij} is the distance between sites *i* and *j* belonging to two different molecules. The parameters of this potential have been chosen so that the total MPE2 interaction energy closely follows the Buckingham potential from the MM3 force field^{81} at all distances. Note that when comparing energy differences at fixed distances, the repulsive part will cancel out, so its particular form is somewhat arbitrary in this case.

Fig. 5 shows the dissociation curve calculated at the MPE2 level, which is exact for the four-electron PPP Hamiltonian used to model the interaction. Naturally, the curve has the qualitative character of the potential energy surfaces of a van der Waals complex. The shapes of MPE2 and Buckingham potentials are very similar at all separations, differing mostly by a slight shift at the short and intermediate range, which justifies the parameterization of the repulsion correction we use. The attractive (MPE2) and repulsive components of the total interaction energy are also plotted in Fig. 5. To analyze the behavior of the attractive part of the total MPE2 (exact) interaction energy, we plot its absolute value on the double logarithmic scale (Fig. 6). Since the attractive component of the MPE2 interaction energy is pure dispersion, it should scale like *R*^{−6} at large separations between monomers. The broken blue straight line in Fig. 6 is a monomial fitted to be tangent to the dispersion energy curve at $ R = 50 \xc5 $. Its slope is −5.98, which matches the asymptotic limit almost perfectly. The exact dispersion energy approaches this limit rather slowly and at $ 10 \xc5 $ the slope is still noticeably smaller. This is expected since the Ohno potential that has been used in the PPP Hamiltonian reaches the 1/*r* dependence of the true Coulomb potential only at infinity. For this type of two-electron interactions, the correct long-range scaling of the dispersion energy is $ ( ( e 4 / U 2 ) + R 2 ) \u2212 3 $ (see Eq. (14)). The red curve in Fig. 6 is a similarly fitted monomial of $ ( e 4 / U 2 ) + R 2 $ plotted as a function of *R*. It follows the exact dispersion energy much closer and around $ 7 \xc5 $ these two curves become visually indistinguishable.

As a second example, we study the model for a parallel stacked benzene dimer. The same parameters have been used as for the ethylene dimer, both in the PPP Hamiltonian and the repulsive energy correction. Assuming that the density is decomposed into maximally localized pairs, two types of partitioning are possible (Fig. 7). Unless otherwise stated, we will assume the partitioning of the *D*_{3h} symmetry, which makes both benzene molecules symmetric. Fig. 8 shows MPE potential energy surfaces up to the 6th order, which is equivalent to the exact diagonalization for the entire system.

As the MPE1 energy adds only a portion of intramolecular correlation energy, the interaction energy is exactly the same as for the EXX, so both curves are on top of each other. MPE2 adds pairwise correlations between electron pairs on two separate benzene molecules and leads to a minimum at $ 4.1 \u2009\xc5 $. Compared to the exact result, the equilibrium distance is about $ 0.1 \u2009\xc5 $ too long and the interaction energy is underestimated by almost 37%; also the long-range decay is too quick. MPE3 corrects this result substantially, yielding the correct equilibrium distance. The binding energy is overestimated by approximately 8% at equilibrium, but the interaction energies are very accurate beyond $ 6 \u2009\xc5 $. Interestingly, MPE4 gives a substantial overcorrection, shifting the minimum again to $ 4.1 \u2009\xc5 $ and underbinding by over 15%. MPE5 is only a slight improvement, while MPE6 is a relatively significant correction leading to the exact result. Overall, already MPE2 captures the qualitative behavior of the PES and locates the equilibrium distance fairly accurately. The deterioration of MPE4 energies compared to MPE3 might be related to the fact that three pairs are conjugated within a benzene ring. We can break the conjugation completely by setting every second transfer integral to zero (Fig. 9); thus generating effectively a Kekule structure of the benzene molecule. The geometry of the dimer is not changed, so all two-electron interactions in the PPP Hamiltonian are the same as in the original dimer model. In this case MPE2 leads to approximately 45% overbinding. MPE3 energies are already close to the reference, indicating that three-pair correlations can be important even in the model with broken conjugation. MPE4 and MPE5 results are visually indistinguishable from the reference values. This example confirms that in the presence of conjugated electron pairs, the convergence of MPE is likely to be more challenging.

The results presented so far were obtained with pair densities on one benzene that are a mirror reflection of the pair densities on the other molecule (left-hand side of Fig. 7). However, MPE results clearly depend on the way how the total density is decomposed. To investigate this effect, we consider also a pattern where pair densities on one benzene are rotated by $ 60 \xb0 $ with respect to the other one (right-hand side of Fig. 7). Comparing Figs. 8 and 10 reveals that the choice of partitioning is decisive for the convergence of MPE results. For the *C*_{3}_{υ} partitioning (Fig. 10), MPE2 produces only a shallow minimum at $ 4.8 \u2009\xc5 $ and the binding energy is underestimated by approximately 80%. MPE3 corrects this result to a major extent, exhibiting a minimum that is only 12% too low, but MPE4, instead of getting closer to the reference, leads to even greater overbinding. MPE5 is again very close to the exact result.

The example of benzene dimer shows that the convergence of MPE interaction energies may be non-monotonic. This points to the fact that MPE may behave less well if density decomposition breaks a system of strongly delocalized electrons. In such a case, the results are also sensitive to the partitioning pattern. Good results obtained with MPE3 suggest that the order of the expansion needs to be somehow balanced with the number of electron pairs that are delocalized.

Interaction energy is non-additive and so is the dispersion energy component.^{82,83} In the model of fluctuating multipole moments, the non-additivity of dispersion results from dynamical screening of polarizabilities by the surrounding subsystems and from the truly many-body nature of interactions between them, which cannot be reduced to a sum of pairwise contributions. Since MPE includes electron correlation through the use of many-determinantal wavefunctions, these effects are naturally accounted for to the degree determined by the number of electrons treated explicitly at any given order of the expansion. Due to its long-range character, dispersion is a collective effect, whereby the total stabilization energy grows superlinearly with the size of the system. These effects contribute to the overall cooperativity of binding of van der Waals complexes. To study how the MPE interaction energy scales with the system size, we consider stacks of ethylene and benzene molecules. To ensure equivalence of sites and thus uniform site densities, we impose cyclic boundary conditions. The distances between monomers are fixed at $ 4.1 \u2009\xc5 $ for ethylene and $ 4.0 \u2009\xc5 $ for benzene stacks. Figs. 11 and 12 show additional stabilization per monomer due to intermolecular dispersion interactions as a function of the number of monomers in the system. To facilitate comparison, plots show relative changes with respect to the interaction energy of a trimer for every MPE level. For the ethylene stack, MPE3 and MPE4 energies are visually indistinguishable, so MPE3 is already the converged result. Cooperativity of dispersion interactions is noticeable as additional molecules in the stack lower the stabilization energy per monomer. This monotonic decrease flats out at 7 monomers for all considered levels of MPE. This indicates that while MPE2 underestimates the additional stabilization by about 25%, it properly scales with the system size. The converged total interaction energy per monomer is 0.76 meV, so the cooperativity effect accounts for about 3% of the total value.

The results for benzene stacks are qualitatively similar. The relative effect of adding MPE3 corrections is more pronounced as MPE2 recovers less than 50% of the additional stabilization. Adding MPE4 correction stabilizes the system further; however, the effect is comparatively small. For 7 and 8 monomers in the stack, some loss of accuracy becomes apparent due to accumulation of numerical errors. Contrary to the total interaction energy in the benzene dimer (Figs. 8–10), additional stabilization due to many-body interactions between different monomers seems to converge monotonically with the level of MPE. As the total MPE4 interaction energy per monomer stabilizes at 4.7 meV, the additional stabilization due to cooperativity accounts for almost 5% of the effect.

These applications of MPE to the PPP Hamiltonian clearly show that already low levels of the expansion successfully recover the dispersion interactions, which are absent in the mean field description. As all the terms are defined as density functionals, the implication is that MPE provides an efficient way to formulate a dispersion correction scheme for any class of approximate exchange-correlation functionals. The procedure is completely non-empirical and naturally accounts for both intra- and inter-molecular dispersion. As it does not assume the asymptotic 1/*R*^{6} dependence, it works equally well at any separation between interacting subsystems and does not require any damping at short ranges. The effect of cooperativity of many-body interactions is properly accounted for and the scheme goes beyond the pair-wise additive correction and effectively includes the dynamical dielectric screening effects. These attractive characteristics warrant pursuing further development of the proposed scheme and incorporating it in practical DFT computations.

## IV. CONCLUSIONS

In this work we have applied MPE to PPP Hamiltonians for several $\pi $-electron systems. The aim was to study the convergence and accuracy of MPE calculations including weak long-range electron correlations. The analysis revealed that energy corrections decay rapidly with the distance between pair densities in the fragment for which correlated calculations are performed. For a large system, this observation allows to neglect a vast number of corrections, leading to reduced scaling at the cost of introducing small error. Additionally, we have shown that the convergence of the expansion depends on the partitioning of the total density. In particular, an optimal partitioning leads to accurate results already at low orders, making MPE calculations practical for large systems. It has been shown that already MPE2 recovers the qualitative and semi-quantitative features of dissociation curves of van der Waals complexes and recovers the 1/*R*^{6} long-range asymptotic behavior of dispersion energy. The importance of three-pair contributions has been shown on the example of stacked benzene molecules. For ethylene and benzene stacks, MPE3 accurately accounts for the cooperativity of dispersion interactions.

This and previous work on the applications of MPE to model Hamiltonians show that relatively low orders of expansion are able to account for both strong and weak long-range correlations, for which common density functional approximations often fail miserably. Missing dispersion is a long-standing problem in DFT and the present work illustrates how these interactions could be systematically incorporated within the DFT framework without any empiricism. This approach goes beyond the typical atomic pair-wise correction schemes and implicitly accounts for non-additive many-body and dynamical screening effects. It is also valid at any separation between subsystems, in particular damping at short distances is not needed. While in this work we focus on model systems, the MPE dispersion correction can find numerous applications to studies of materials and biological systems, for which dispersion interactions are of vital importance for their functionality.

The future work on MPE should focus on the implementation for *ab initio* Hamiltonians, for which the decomposition of the density into smooth pair densities is the starting point. An efficient implementation will rely on the observation that contributions from distant pairs can be often neglected, which would substantially reduce the scaling of the method. Given the smooth long range decay of the dispersion term, the computational cost could be reduced further by combining MPE with a more efficient non-local correlation functional like VV10^{34} or vdwDF.^{38,39} Finally, this work points out the potential problems of MPE, which is sensitivity on the density partitioning. Finding a systematic way to break up the density in an optimal way would be highly desirable.

## ACKNOWLEDGMENTS

This work was funded by a grant from the NSF (No. CHE-1464804), T.V.V. is a David and Lucile Packard Foundation Fellow.