###### Abstract

We present a formally exact van der Waals inclusive electronic structure theory, called FDE-vdW, based on the Frozen Density Embedding formulation of subsystem Density-Functional Theory. In subsystem DFT, the energy functional is composed of subsystem additive and non-additive terms. We show that an appropriate definition of the long–range correlation energy is given by the value of the non-additive correlation functional. This functional is evaluated using the Fluctuation–Dissipation Theorem aided by a formally exact decomposition of the response functions into subsystem contributions. FDE-vdW is derived in detail and several approximate schemes are proposed, which lead to practical implementations of the method. We show that FDE-vdW is Casimir–Polder consistent, i.e. it reduces to the generalized Casimir–Polder formula for asymptotic inter-subsystems separations. Pilot calculations of binding energies of 13 weakly bound complexes singled out from the S22 set show a dramatic improvement upon semilocal subsystem DFT, provided that an appropriate exchange functional is employed. The convergence of FDE-vdW with basis set size is discussed, as well as its dependence on the choice of associated density functional approximant.

1.4

FDE-vdW: A van der Waals Inclusive Subsystem Density-Functional Theory

[3ex]

Ruslan Kevorkyants, Henk Eshuis, and Michele Pavanello

Department of Chemistry, Rutgers University, Newark, NJ 07102, USA

Department of Chemistry and Biochemistry, Montclair State University, Montclair, NJ 07043, USA

[2ex]

Date: | July 27, 2019 |

Status: | Submitted to J. Chem. Phys. |

## 1 Introduction

In subsystem Density-Functional Theory [1] (subsystem DFT), the total electronic density of a molecular system is partitioned into subsystem contributions:

(1) |

where is the total number of subsystems. This is an appealing strategy, as it invokes a kind of “divide and conquer” principle which makes the electronic problem more manageable by partitioning it into several coupled sub-problems. This strategy is general, and can be formulated in such a way as to involve no approximations [2, 3]. In practice, this divide and conquer approach is more easily implemented if the electron system is composed of weakly interacting subsystems [4, 5]. In this limiting case, the total electron density in Eq. (1) is already well approximated by using the electron densities of isolated subsystems (i.e. the molecular fragments composing the full system).

Subsystem DFT is designed (see Section 2) to reproduce Kohn–Sham DFT [6] (KS-DFT), and thus is formally exact [4]. Because of this, subsystem DFT inherits all the strengths and weaknesses of KS-DFT, e.g. the fact that the exact exchange–correlation (XC) functional, , is unknown. As a consequence, practical implementations of subsystem and KS-DFT involve employing XC functional approximants. The most common class of approximants rely on local and semilocal forms of the functional, i.e. it depends on the electron density and possibly also on its gradient at a point in space. Due to their character, these semilocal approximations are unable to account for interactions that are non-local in nature, such as dispersion interactions [7, 8] and generally all long–ranged interactions originating from the correlated part of the energy functional.

Several solutions have been proposed to address this deficiency of local and semilocal XC functionals [9]. One avenue is the design of fully non-local functionals, as proposed by Dion et al. [10] and by Vydrov and Van Voorhis [11]. A very popular and efficient alternative has been an additive pairwise correction that uses the so-called coefficients [12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. This type of correction is inspired by the Casimir-Polder formula,

(2) |

where is the distance between fragment and in a molecular system. Appropriate decompositions of the molecular system are used to yield pairwise corrections [22, 23, 24].

A better description of dispersion interactions is given by the generalized Casimir-Polder (GCP) formula [7, 8]

(3) |

where and are the electronic frequency dependent linear-response functions of the two isolated fragments. Here, the fragment’s electron densities are allowed to overlap and the restriction to large interfragment distance is also lifted. The GCP formula is derived to account for dispersion interactions to second order in perturbation theory, where the perturbation is the Coulomb interaction between the fragments.

A formally exact framework to obtain the full correlation energy (as opposed to only the dispersion part) is the adiabatic–connection fluctuation-dissipation theorem applied to DFT (ACFDT-DFT) [25, 26, 27] which relates the total correlation energy to the frequency dependent linear-response functions (response functions, hereafter). Since response functions can be efficiently calculated using semilocal DFT [28, 29, 30, 31, 25], ACFDT-DFT is a very appealing way to introduce dispersion interactions, and generally all interactions arising from the correlated part of the energy.

The simplest and most used method based on the ACFDT is the random phase approximation (RPA) [32, 33, 34, 35, 36] and its extensions [37, 38]. RPA is parameter free, and was shown to be an order of magnitude more accurate than semi-local DFT for mid- and long-range correlation [39] and of comparable accuracy as dispersion-corrected functionals, such as B3LYP-D3, or highly parameterized functionals, such as M06-2X [40]. For short-range correlation, RPA does not improve over semi-local DFT, thus range-separated schemes were proposed which combine a short-range DFT description of the system with a long-range RPA description. One advantage of these methods is a much faster convergence with basis set size at the cost of having to include a parameter in the formalism [41, 42, 43].

The main idea behind this work is to devise a theory that exploits the fragment picture provided by subsystem DFT and borrows ideas from ACFDT-DFT to obtain the (theoretically exact) correlation energy of interaction between fragments. This has, so far, not been attempted for overlapping electronic fragment densities (for non-overlapping densities see Ref. 27). As we exploit the Frozen Density Embedding (FDE) formulation of subsystem DFT, we call the resulting theory FDE-vdW.

FDE-vdW has several appealing properties: (1) full correlation energies are recovered, (2) the theory is formally exact and independent of the definition of the fragments (or subsystems), (3) the long-range correlation energy is rigorously incorporated in the subsystem DFT energy functional without resorting to an ad-hoc function that switches between the short and the long range.

Two recent publications provide a major motivation for this work. First, Visscher and coworkers [44], demonstrated that the FDE coupled with semiempirical dispersion corrections yields accurate interaction energies. Second, Chałasiński and coworkers [45], used a projector-based partitioning of HF and KS-DFT wavefunctions into fragment wavefunctions to define an inter-fragment dispersion correction calculated with the GCP formula which yielded excellent interaction energies. In both of the mentioned publications as well as in this work, the common double-counting of van der Waals interactions occurring when correcting the KS-DFT interaction energies is completely bypassed and the resulting theories are more balanced.

The paper is organized as follows. Section 2 recalls some basic ideas of subsystem DFT and its time-dependent extension. Section 3 develops the FDE-vdW theory by applying ACFDT-DFT to subsystem DFT. Section 4 gives an outline of possible approximations for practical calculations. Section 5 presents details of the FDE-vdW implementation. Section 6 displays pilot calculations carried out with one of the approximations proposed in Section 4. Finally, Section 7 delineates conclusions and future directions for the presented formalism.

## 2 Subsystem DFT and TD-DFT

The density partitioning in Eq. (1) is in principle arbitrary. One can use chemical intuition to identify subsystems (using for example a real-space partitioning), or any other criterion. Unlike the density, the energy functional is not additive, and an exact form of it involves non-additive corrections that arise at the kinetic-energy and XC energy level as well as a Hartree term describing the electrostatic interaction between subsystems,

(4) |

where , and in we have grouped the nuclear attraction term. A computationally attractive way to find subsystem densities [5] is the Frozen Density Embedding theory (FDE). In FDE the energy functional of all the subsystem densities is minimized in Eq. (4) w.r.t. only one subsystem density while keeping all the others frozen. This procedure can be repeated until self-consistency [46]. Freezing all densities but one during the minimization, is equivalent to using partial functional derivatives in the working KS equations and with that, mixed derivatives such as are always set to zero in FDE.

The aim of FDE is to represent the subsystems as a set of coupled Kohn–Sham systems. Hence, the subsystem densities must be non-negative, must integrate to a preset number of electrons, i.e. , and must be non-interacting -representable. In practical FDE calculations, the subsystem densities are constructed from subsystem molecular orbitals which are expanded in terms of localized atomic orbitals, often centered on atoms belonging to only one molecular fragment in a system (monomer basis set). Thus, the subsystems are identified as non-covalently bound molecules. Recently, however, embedding schemes based on optimized effective potential [47, 48, 49], projection operators [45, 50, 51, 52, 53], and fitted functionals [54, 55], are paving the way for the correct description of covalent bonds (high density overlap) by FDE.

Self-consistent solution of a set of coupled KS-like equations (also called KS equations with constrained electron density [4]) yields the set of subsystem KS orbitals

(5) |

with the effective subsystem potential given by

(6) |

In the FDE formulation of subsystem DFT [5, 4], the potential above, , is called embedding potential and is given by

(7) |

In the case of closed-shell subsystems, the density of the supersystem is thus found using Eq. (1) and Eq. (5) as

(8) |

The time-dependent extension of subsystem DFT [56, 57, 58, 59] is based on the Runge-Gross theorem [60] which was recently proved also for subsystem DFT by Wasserman et al. [61]. Recently, linear-response subsystem TD-DFT was formulated in terms of subsystem response functions by one of us [62]. The relevant Dyson-type equations for the subsystem response functions read as follows (contracting the integrals in a short-hand notation)

(9) |

where the correlated coupled subsystem response function, , is given by

(10) |

and the correlated uncoupled subsystem response function, , is found by solving the Dyson equation

(11) |

In the above equation, is the response of the KS subsystem. The kernel matrix, , in Eq. (9) and Eq. (11) is defined as

(12) |

where the kinetic energy kernels, expressed in the time domain, are defined as

(13) | ||||

(14) |

## 3 Fluctuation–Dissipation Theorem Subsystem DFT

The correlation energy within the Adiabatic-Connection Fluctuation–Dissipation Theorem DFT (ACFDT-DFT) can be written exactly as

(15) |

where is the -th spin-spatial coordinate, and is the correlated part of the electronic pair density which is defined as [25, 27]

(16) |

The functions and are the electron’s Kohn-Sham and interacting response functions, respectively. The interacting response is evaluated at coupling strength and the adiabatic connection formula is used [27, 63]. The central problem is finding analytical and/or computationally appealing ways to evaluate Eq. (15–16).

As the response functions explicitly depend on both occupied and virtual orbitals, for RPA and TD-DFT alike, a formal scaling with system size arises in their computation. For realistic systems this leads to a large computational expense and it is usually handled employing numerical techniques based on fitting of the response functions or the entire correlation energy expression [64, 35, 65]. An important aspect of this work is that the computation of Eq. (15–16) for the supersystem is completely avoided. We elaborate on this crucial aspect in the next section.

### 3.1 Short range, subsystem additive correlation energy

We define the short range, intra-subsystem correlation energy as the correlation energy derived with Eq. (15) with the pair density including only the response functions calculated in the uncoupled formalism of Eq. (11). This definition is inspired by the idea that when the subsystems are separated by an infinitely large distance, we have (where is the response function of the isolated subsystem). In this asymptotic case, the short range correlation energy depends only on the degrees of freedom of a single subsystem, and the correlated part of the interaction between subsystems vanishes.

According to our definition, the additive (or short range) part of the correlated pair density, , to be used in Eq. (16) depends on the difference of response functions

(17) |

where we have used the approximation , which was discussed previously [62]. To obtain the above uncoupled subsystem response function one needs to use Eq. (11) and evaluate the interaction kernel at coupling strength ,

(18) |

Even though the additive correlation energy can theoretically be calculated with this formulation, it is not computationally advantageous [33, 64] and poses problems related to a singular as due to the (semi) locality of commonly adopted XC functionals [25, 66]. Hereafter, we will make the assumption that the additive (or short-range) correlation is well captured by semilocal XC density functional approximants. Thus, the remainder of this work focuses on the inter-subsystem correlation energy, by finding an appropriate expression for the long range (or non-subsystem-additive) part of the correlated pair density, .

### 3.2 Long range non-additive correlation energy

The subsystem DFT density functional in Eq. (4) is comprised of subsystem-additive and non-additive terms. Specifically, the non-additive XC functional can be further split into exchange and correlation contributions, . The long range contribution to the correlation energy is defined as the energy term left subtracting the short-range contribution from the total correlation energy,

(19) |

The difference of response functions that yields the non-additive part of the correlated pair density is

(20) |

which by inspection of Eq. (9) is non-subsystem-additive. Thus, the long range correlation energy derived from Eq. (20) will naturally replace the non-additive correlation energy functional in Eq. (4). It is easily shown that

(21) |

In order to obtain the non-additive correlated pair density, integrals in both the coupling strength () and frequency () of the difference of subsystem response functions defined above need to be evaluated

(22) |

Thus, the long range component of the correlation energy can be calculated substituting Eq. (22) into Eq. (15). Eq. (22) is the key result of this work. Employing Eq. (9) in Eq. (20–22) results in

(23) |

## 4 Approximate treatments of the non-additive correlation energy

The non-additive correlation energy determined with Eq. (23) is still computationally prohibitive for realistic systems, as the coupled response functions must be obtained solving the Dyson-type equation in Eq. (9) which in theory requires operations, where is the size of the entire supersystem [62, 67]. It is, therefore, important to develop simplified expressions to compute the non-additive correlation energy. The first approximation is obtained by considering a perturbative solution to the Dyson-type equation Eq. (9) as

(24) |

and retaining only the first term of the expansion. Eq. (23) then simplifies to

(25) |

The non-additive correlation energy is now expressed in terms of the functions (i.e. the uncoupled subsystem response functions). In practice Eq. (25) has to be further approximated since is not known exactly. Motivated by the success of RPA for non-covalent interactions, we approximate by neglecting the frequency-dependent exchange-correlation kernel as well as the kinetic energy contributions to the kernel (, , and )

(26) |

Since the RPA kernel is frequency independent, it is computationally much more efficient to solve Eq. (25). We dub the approximate method GCP. One can reduce the computational cost further by applying subsequent approximations. We obtain GCP by evaluating the coupling strength integration for subsystems evaluated at . Finally, a further approximation can be carried out (although we never do in this work) by which the subsystems can be treated as isolated (GCP). We summarize the various approximations and their proposed acronyms in Table 1.

Description | Integrand of Eq. (25) | Acronym |
---|---|---|

Exact | ||

RPA & Perturbative | GCP | |

& RPA & Perturbative | GCP | |

Isolated-fragment response functions & RPA & Perturbative | GCP |

For only two interacting subsystems (), application of the GCP approximation and subsequent integration over (which gives a factor ) yields

(27) |

The use of a coupling-strength independent response function is analyzed in detail and justified in the appendix section.

An important test of the meaningfulness of the additive and non-additive correlation energy definitions given in this work can now be carried out by analyzing the behavior of the non-additive correlation in the limit of large subsystem separations (or equivalently negligible intersubsystem density overlap). In this limit, the correct espression of the van der Waals interaction is the generalized Casimir-Polder or Longuet-Higgins formula (first derived by McLachlan [7] for molecules) given in Eq. (3) of the introduction section. The GCP formula is analogous to Eq. (27) and is recovered by adopting the “GCP” approximation listed in Table 1. As mentioned before, for large inter-subsystem separations, . Thus, in the asymptotic limit, this theory reduces to the GCP formula – which is an important property, also termed Casimir-Polder consistency [68]. The GCP formula has a long history, both in wavefunction-based and DFT methods for correcting the mean-field approximation [69, 28, 70, 71, 72].

## 5 Practical Calculations

### 5.1 Integration of the fluctuation-dissipation formula

For closed–shell systems, the spectral representation of the response functions, assuming real solutions of the TD-DFT eigenvectors, is given by [31]

(28) | |||

(29) |

where are the KS orbitals of subsystem . Subscripts denote occupied orbitals and virtual orbitals. is the projection of the sum of the excitation () and de-excitation () TD-DFT eigenvector for the -th excited state of subsystem . The associated eigenvalue is given by . The Hartree–XC kernel used to determine is given by Eq. (12) with [56, 73]. Using Eq. (28) and employing the GCP approximation (see Table 1) we can express the difference of the coupled and uncoupled response functions

(30) |

Integrating over coupling strength and taking into account double counting arising from the summation over the subsystems, we obtain for Eq. (23)

(31) | |||

Finally, the frequency integration can be carried out analytically [32],

(32) | |||

### 5.2 Implementation of the GCP method

Eq. (5.1) was implemented in the Adf computer software [74] as a modification of the existing SUBEXCI code [75, 67, 57, 59]. The SUBEXCI code uses the so-called -vector (Hermitian) formulation of the TD-DFT equations. Eq. (5.1) is transformed to

(33) |

In the above,

(34) |

where , the vectors are the solutions of the Hermitian TD-DFT equations and, following the GCP approximation, in Eq. (34) is modified to only include the Coulomb kernel.

As Adf is a Slater-Type Orbital code, two-electron integrals are not available. For this reason, the matrix elements are calculated in density fitting. For each matrix element, the inducing density is fitted to yield a fitted induced potential, . The induced potential is then integrated over a monomer-based integration grid with . This is not a very efficient code for our purposes. In the recent works by Szalewicz et al. [65] regarding the implementation of the GCP formula, density fitting is applied to the entire response function before carrying out the frequency integration. Eventually, our Adf code will include a Szalewicz-type density fitting. However, as this work constitutes a proof-of-principle, we have left additional coding for a future work.

### 5.3 Binding energy combining the non-local correlation with the FDE energy: The FDE-vdW method

In FDE, the total binding energy of two subsystems can be calculated by subtracting from the total FDE energy the energy of the isolated subsystems [76]. In the case of only two subsystems, we have

(35) |

where we have introduced indicating the electron density of the isolated subsystem . Pernal et al. [77] have devised a strategy to first remove the dispersion interactions from the density functional, and then to augment the obtained binding energy by the non-local dispersion interaction obtained from the GCP formula. Rajchel et al. [45], instead, are able to compute dispersionless interactions by partitioning the HF or KS wavefunctions into fragment wavefunctions. Here we largely follow these ideas. However, instead of reparametrizing the density functional to remove the dispersion interactions from it, after the FDE procedure has completed and subsystem densities are recovered, we calculate the FDE energy omitting the correlation part of the non-additive XC functional. In this way, we completely remove the correlation part of the interaction energy. We then add back the correlation energy by adding to the “correlationless” binding energy the non-additive correlation obtained with the GCP. The resulting binding energy expression defines the FDE-vdW method and is given by

(36) |

It is important to remark that the subsystem DFT formalism greatly simplifies the theory by providing an algorithm that naturally includes the non-additive correlation energy in the DFT energy expression, without invoking approximations in the underlying theory. Another important observation is that the theory so far ignores the relaxation of the subsystem densities due to the non-additive, long range interaction. Since these interactions are typically small for non-covalently bound subsystems, this is a reasonable approximation.

### 5.4 Computational details

FDE and FDE-vdW binding energies are evaluated according to Eq. (5.3), and Eq. (36), respectively. FDE and FDE-vdW calculations are carried out employing the revPBE exchange density functional [78] and PBE correlation [79] unless otherwise stated. Eq. (11) is solved within the adiabatic approximation using the semi-local XC kernel corresponding to the density functional employed. The FDE non-additive kinetic energy functional PW91k [80] is used throughout. The calculations are carried out with a modified version of ADF [81], and employ a TZ2P Slater-type orbital basis set unless otherwise stated.

All excitations were included in the sum in Eq. (33) with exception of the large monomers (i.e. complexes containing a benzene monomer), for which excitations were sorted in descending order of contribution to the molecular polarizability, similarly to Ref. [82, 57]. Only the excitations accounting for 99.5% of the isotropic, static polarizability were used to compute Eq. (33), thereby reducing computational cost while maintaining a satisfactory level of accuracy.

## 6 Results and discussion

### 6.1 FDE vs. FDE-vdW binding energies

We calculated the binding energies of a selected set of weakly bound complexes from the S22 set [83] using GCP. The S22 set is widely used to benchmark the performance of methods. It consists of a set of dispersion bound systems, a set of hydrogen-bound systems, and a set which is both dispersion and hydrogen bound. Because the algorithm, as implemented, scales as , we are limited to apply it to relatively small test cases. We therefore only selected complexes from the S22 set with less than 20 atoms as well as the two benzene dimer structures, one with stacked monomers [(CH)(CH)] and one T-shape [(CH)(CH)]. This choice narrowed the set to 13 complexes out of 22. A summary of the resulting binding energies is found in Table 2.

Complex | FDE(1) | FDE-vdW(1) | FDE-vdW(2) | CCSD(T)/CBS |
---|---|---|---|---|

(NH) | 2.41 | 3.34 | 3.34 | 3.13 |

(HO) | 3.39 | 4.28 | 4.25 | 4.99 |

(HCOOH) | 12.32 | 16.70 | 16.54 | 18.75 |

(HCONH) | 12.73 | 16.77 | 16.74 | 16.06 |

(CH) | 0.01 | 0.58 | 0.58 | 0.53 |

(CH) | 0.07 | 1.28 | 1.27 | 1.47 |

(CH)(CH) | 1.01 | 1.66 | 1.67 | 1.50 |

(CH)CH | 0.24 | 1.74 | 1.73 | 1.45 |

(CH)HO | 1.04 | 3.39 | 2.58 | 3.28 |

(CH)NH | 0.35 | 2.92 | 2.07 | 2.31 |

(CH)HCN | 1.67 | 4.55 | 3.83 | 4.54 |

(CH)(CH) | 0.89 | 2.64 | 2.15 | 2.65 |

(CH)(CH) | 0.34 | 1.88 | 2.00 | 2.72 |

MUE | 2.31 | 0.46 | 0.57 |

from Ref. 94

The results presented in Table 2 indicate that the FDE-vdW method is consistently closer to the benchmark than semilocal FDE. The MUE is decreased from 2.31 kcal/mol when the non-additive correlation is evaluated with PBEc, to 0.46 and 0.57 when the GCP correlation is evaluated either with the revPBE or the SAOP model potential, respectively.

We now discuss the effect of the asymptotic behavior of the XC potential in the TD-DFT calculations. In the ideal case of employing a complete basis set, the GCP non-additive correlation energy overestimates the correlation energy computed with Eq. (23) for two reasons: (1) GGA functionals tend to underestimate excitation energies, which yields an overestimation of the non-additive correlation energy, since Eq. (5.1) contains excitation energies in the denominator, and (2) from arguments related to the approximate coupling-strength integration (see the appendix and Figure 2). Point (1) is known and currently dealt with in the literature for implementations of the GCP formula by employing hybrid XC potentials in combination with asymptotic corrections [65, 77, 84, 69, 28]. We also observe that in moving from revPBE to the asymptotically corrected SAOP model potential [85, 86, 87], the FDE-vdW binding energies generally decrease in magnitude.

Since the calculations were performed with a medium-sized basis set (the TZ2P set, see next section for a basis set analysis), which truncates the sum-over-states formula in Eq. (33), we are bound to profit from cancellation of error – although GCP should overestimate the correlation energy interaction, the finite basis set used will offset such overestimation and possibly result in an underestimation. The choice of basis set is analyzed in the following section.

### 6.2 Choice of basis set

It is well-known [33] that correlation energies calculated from RPA or TD-DFT response functions are significantly dependent on the choice of basis set. In the case of the GCP method, there is consensus that the larger the basis set used, the larger the calculated dispersion interaction [28]. This is what we also find for the GCP method as exemplified by the calculations presented in Table 3.

System | (NH) | (HO) | (HCOOH) | (HCONH) | (CH)(CH) |
---|---|---|---|---|---|

(TZ2P) | 1.57 | 1.64 | 7.13 | 5.97 | 0.99 |

(QZ3P) | 1.80 | 2.01 | 8.65 | 7.12 | 1.35 |

# functions (TZ2P) | 132 | 108 | 228 | 252 | 192 |

# functions (QZ3P) | 176 | 144 | 304 | 336 | 256 |

This behaviour can be rationalized by looking at Eq. (5.1). The larger the basis set, the more excitations one needs to include in the sum and thus a larger non-additive correlation is found, up to an asymptote. Assuming that the QZ3P calculations are close to the basis set limit (as we use Slater-Type Orbitals featuring correct exponential decay and electron–nucleus cusps, faster convergence of the calculations with respect to basis set size is expected than for Gaussian-Type functions), the values in Table 3 show that already the TZ2P set provides most of the non-additive correlation energy. Since only the long-range non-additive part of the correlation energy is calculated using GCP, we expect much less dependence on basis set size than for full RPA correlation energies. This is also observed in so-called range-separated RPA schemes [42, 43].

### 6.3 Choice of functional

Some approximate exchange functionals, such as the PW91 and PBE functionals, may artificially display minima for van der Waals complexes [88]. Others, such as the revPBE and the BP88, do not feature minima for van der Waals complexes and are often chosen to be the exchange counterpart of van der Waals density functionals [89, 90, 88].

In this work, we carry out a preliminary analysis of the GGA XC functional used in the FDE calculations, which then enters the FDE/FDE-vdW binding energy expression in Eq. (5.3–36). The binding energy of the van der Waals complexes featured in Table 2 was recalculated employing 7 XC density functional approximants. The results are summarized in Figure 1. In the figure, it is evident that PW91, PBE, BLYP, and RPBE introduce spurious attractive interaction between the subsystems when they are employed as non-additive exchange functionals. For this reason, they are not suitable for pairing with a van der Waals correction, as previously noted for KS-DFT calculation [89, 90, 77] and here is also noted in FDE-vdW calculations.

To conclude the XC functional analysis, the appropriate functionals for pairing with the GCP among the tested ones are: HTBS, revPBE, and BP86. In a future study, we will carry out a complete benchmark and a larger set of GGA functionals will be considered.

## 7 Conclusions

We have developed a new van der Waals inclusive DFT theory in the framework of subsystem DFT. Subsystem DFT is an ideal tool for identifying long and short range contributions to the binding energy, as they are associated respectively with the value of the additive and non-additive energy functionals. This is an important feature of the present theory which completely avoids the common parametrization of density functionals to classify short and long range interactions. Such parametrizations are often used to expedite the computational evaluation of van der Waals interactions [42].

As a result of this work, the needed equations to evaluate van der Waals interactions arising between overlapping electron densities are now available. This has been an outstanding issue in this field, as a DFT theory for van der Waals interactions between two distinct electron densities has been available only in the non-overlapping case [27].

A recent energy decomposition analysis of subsystem DFT calculations [44] has shown that the inclusion of empirical van der Waals terms in the interaction energy has the effect of improving the binding energies calculated with GGA functionals. In line with those findings, we show that the FDE binding energies are dramatically improved if the non-additive correlation energy functional is made non-local in character and is evaluated with the Fluctuation–Dissipation Theorem.

Preliminary results are promising and lead us to expect that FDE-vdW can deliver good binding energies across three binding regimes: dispersion, dipole, and mixed dispersion-dipole.

A very appealing property of FDE-vdW is its potential for systematic improvement. For example, by including a numerical integration over the coupling strength parameter, , which would fold in non-local interaction effects beyond RPA. Another interesting improvement is to use the second term in the perturbative expansion in Eq. (24). In this way, three-body terms would be included in the overall non-additive correlation.

There are two clear weak points in the current implementation of the theory. First, the computational scaling for evaluating the non-additive correlation is prohibitive and it is typical of a correlated wavefunction method. In addition, the exact diagonalization of the Casida equations needed for the evaluation of the subsystem response function is only possible for small systems. Larger systems, such as uracil, indol and pyrazine, incur into linear dependencies of the occupied-virtual product space that prevents a numerically stable diagonalization. The second weak point is that, in subsystem DFT, kinetic energy functionals are needed (at the non-additive level) for allowing a proper description of overlapping electron densities. The available kinetic energy functionals are mostly formulated at the GGA level [91, 80, 92, 1] which are known to describe well only certain types of weak interactions and not others [93, 76] (although GGA XC functionals had been used in the referred studies).

Thus, the present theory will need to be accompanied in the future with more accurate kinetic energy density functionals, and more computationally amenable procedures for evaluating the response functions and the associated non-additive correlation energy. Efforts to address both points are ongoing in our group and are based on the formulation of non-local kinetic energy functionals to address the former issue and on a density fitting approach to address the latter which will be inspired by the techniques developed by Szalewicz and coworkers [65].

## Acknowledgements

We are grateful to Alisa Krishtal for the feedback provided on an early version of the manuscript. Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund and the office of the Dean of FASN of Rutgers-Newark for partial support of this research.

## Appendix A Justification of the GCP approximation

In this section, we will show that the integrand of Eq. (25) after all integrals are carried out with the exception of the one over the coupling strength, , has positive and approximately monotonically decreasing derivative. Thus, if the response functions are chosen at a fixed value of , a good choice is .

### a.1 Derivative of the correlation energy with respect to the coupling strength

Let us define as the value of the integrand in the definition of in Eq. (25) when the imaginary frequency and the spatial integrations are carried out. Consider the first derivative of the integrand of Eq. (25), contracting the frequency and the spatial dependence,

(37) |

Realizing that

(38) |

and evaluating Eq. (37) in the RPA approximation, i.e.

(39) |

leads to

(40) |

which could have been achieved by simple evaluation of the derivative through finite differences of the response functions. Realizing that the terms involving the response function difference are one-order-of-magnitude smaller than the other term, the above equation leads to the following approximate derivatives

(41) |

In the above equation, by the short-hand notation “Tr” we mean integration over all variables with exception of .

### a.2 Linear extrapolation of

Recalling Eq. (39), and also that the evaluation of the integral in Eq. (25) with the KS response functions () yields much larger correlation energies than with correlated response functions [71], we find the following approximate inequality

(42) |

Because the integrand of Eq. (25) at zero coupling strength is exactly zero, the trend of can be inferred. Namely, always lies above the linear extrapolation line, see Figure 2. For comparison, see Figure 3 in Ref. 32 keeping in mind that the plot there is inverted on the -axis when comparing it to our Figure 2.

The above analysis, shows that a good approximation to the true is simply the linear extrapolation, which is realized by Eq. (27), using the RPA approximation.

From the graph, it is also inferred that the linear extrapolation of the integrand using as slope the derivative at a fixed will yield overestimated correlation energies, and that the largest overestimation is achieved when , and a much reduced overestimation is found if .

## References

- [1] C. R. Jacob and J. Neugebauer, WIREs: Comput. Mol. Sci. , Accepted for publication (2013).
- [2] M. H. Cohen and A. Wasserman, J. Phys. Chem. A 111, 2229 (2007).
- [3] O. V. Gritsenko, On the Principal Difference Between the Exact and Approximate Frozen-Density Embedding Theory, in Recent Advances in Orbital-Free Density Functional Theory, edited by T. A. Wesolowski and Y. A. Wang, chapter 12, pages 355–365, World Scientific, Singapore, 2013.
- [4] T. A. Wesolowski, One-Electron Equations for Embedded Electron Density: Challenge for Theory and Practical Payoffs in Multi-Level Modeling of Complex Polyatomic Systems, in Computational Chemistry: Reviews of Current Trends, edited by J. Leszczynski, volume 10, pages 1–82, World Scientific, Singapore, 2006.
- [5] T. A. Wesolowski and A. Warshel, J. Chem. Phys. 97, 8050 (1993).
- [6] P. Hohenberg and W. Kohn, Phys. Rev. 136, 864 (1964).
- [7] A. McLachlan, Proc. R. Soc. A 271, 387 (1963).
- [8] H. C. Longuet-Higgins, Discuss. Faraday Soc. 40, 7 (1965).
- [9] J. Klimeš and A. Michaelides, J. Chem. Phys. 137, 120901 (2012).
- [10] D. C. Langreth, M. Dion, H. Rydberg, E. SchrÃ¶der, P. Hyldgaard, and B. I. Lundqvist, Int. J. Quantum Chem. 101, 599 (2005).
- [11] O. A. Vydrov and T. Van Voorhis, J. Chem. Phys. 132, 164113 (2010).
- [12] Y. Y. Sun, Y.-H. Kim, K. Lee, and S. B. Zhang, J. Chem. Phys. 129, 154102 (2008).
- [13] S. Grimme, J. Comp. Chem. 25, 1463 (2004).
- [14] M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, and E. Kaxiras, J. Chem. Phys. 114, 5149 (2001).
- [15] P. Jurecka, J. Cerny, P. Hobza, and D. R. Salahub, J. Comp. Chem. 28, 555 (2007).
- [16] T. Sato, T. Tsuneda, and K. Hirao, Mol. Phys. 103, 1151 (2005).
- [17] O. A. von Lilienfeld, I. Tavernelli, U. Rothlisberger, and D. Sebastiani, Phys. Rev. Lett. 93, 153004 (2004).
- [18] T. Sato and H. Nakai, J. Chem. Phys. 133, 194101 (2010).
- [19] A. D. Becke and E. R. Johnson, J. Chem. Phys. 127, 154108 (2007).
- [20] A. Krishtal, D. Geldof, K. Vanommeslaeghe, C. V. Alsenoy, and P. Geerlings, J. Chem. Theory Comput. 8, 125 (2012).
- [21] S. N. Steinmann and C. Corminboeuf, J. Chem. Phys. 134, 044117 (2011).
- [22] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
- [23] A. Tkatchenko, R. A. DiStasio, R. Car, and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012).
- [24] F. Rob and K. Szalewicz, Chem. Phys. Lett. 572, 146 (2013).
- [25] F. Furche and T. Van Voorhis, J. Chem. Phys. 122, 164106 (2005).
- [26] W. Kohn, Y. Meir, and D. E. Makarov, Phys. Rev. Lett. 80, 4153 (1998).
- [27] J. F. Dobson and T. Gould, J. Phys.: Condens. Matter 24, 073201 (2012).
- [28] A. J. Misquitta and K. Szalewicz, Chem. Phys. Lett. 357, 301 (2002).
- [29] A. Heßelmann and G. Jansen, Chem. Phys. Lett. 362, 319 (2002).
- [30] M. E. Casida, Time-Dependent Density Functional Response Theory for Molecules, in Recent Advances in Density Functional Methods Part I, edited by D. P. Chong, pages 155–192, World Scientific, Singapore, 1995.
- [31] F. Furche, J. Chem. Phys. 114, 5982 (2001).
- [32] F. Furche, J. Chem. Phys. 129, 114105 (2008).
- [33] H. Eshuis and F. Furche, J. Chem. Phys. 136, 084105 (2012).
- [34] X. Ren, P. Rinke, C. Joas, and M. Scheffler, J. Mater. Sci. 47, 7447 (2012).
- [35] H.-V. Nguyen and S. de Gironcoli, Phys. Rev. B 79, 205114 (2009).
- [36] H.-V. Nguyen and G. Galli, J. Chem. Phys. 132, 044109 (2010).
- [37] J. E. Bates and F. Furche, J. Chem. Phys. 139, 171103 (2013).
- [38] J. Paier, X. Ren, P. Rinke, G. E. Scuseria, A. Grüneis, G. Kresse, and M. Scheffler, New J. Phys. 14, 043002 (2012).
- [39] H. Eshuis and F. Furche, J. Phys. Chem. Lett. 2, 983 (2011).
- [40] Y. Zhao, N. E. Schultz, and D. G. Truhlar, J. Chem. Theory Comput. 2, 364 (2006).
- [41] F. Bruneval, Phys. Rev. Lett. 108, 256403 (2012).
- [42] J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. Ángyán, Phys. Rev. Lett. 102, 096404 (2009).
- [43] J. Toulouse, F. Colonna, and A. Savin, Phys. Rev. A 70, 062505 (2004).
- [44] S. M. Beyhan, A. W. Götz, and L. Visscher, J. Chem. Phys. 138, 094113 (2013).
- [45] L. Rajchel, P. S. Żuchowski, M. M. Szczȩśniak, and G. Chałasiński, Phys. Rev. Lett. 104, 163001 (2010).
- [46] T. A. Wesolowski and F. Tran, J. Chem. Phys. 118, 2072 (2003).
- [47] J. D. Goodpaster, N. Ananth, , F. R. Manby, and T. F. Miller, III, J. Chem. Phys. 133, 084103 (2010).
- [48] J. D. Goodpaster, T. A. Barnes, and T. F. Miller, III, J. Chem. Phys. 134, 164108 (2011).
- [49] C. Huang, M. Pavone, and E. A. Carter, J. Chem. Phys. 134, 154110 (2011).
- [50] Y. G. Khait and M. R. Hoffmann, On the orthogonality of orbitals in subsystem kohn-sham density functional theory, in Ann. Rep. Comp. Chem., edited by R. A. Wheeler, volume 8 of Ann. Rep. Comp. Chem., pages 53 – 70, Elsevier, 2012.
- [51] T. A. Barnes, J. D. Goodpaster, F. R. Manby, and T. F. Miller, J. Chem. Phys. 139, 024103 (2013).
- [52] F. R. Manby, M. Stella, J. D. Goodpaster, and T. F. Miller, J. Chem. Theory Comput. 8, 2564 (2012).
- [53] J. D. Goodpaster, T. A. Barnes, F. R. Manby, and T. F. Miller, J. Chem. Phys. 137, 224113 (2012).
- [54] M. J. Gillan, D. Alfe, P. J. Bygrave, C. R. Taylor, and F. R. Manby, J. Chem. Phys. 139, 114101 (2013).
- [55] X. Hu, Y. Jin, X. Zeng, H. Hu, and W. Yang, Phys. Chem. Chem. Phys. 14, 7700 (2012).
- [56] M. E. Casida and T. A. Wesolowski, Int. J. Quantum Chem. 96, 577 (2004).
- [57] J. Neugebauer, Phys. Rep. 489, 1 (2010).
- [58] J. Neugebauer, J. Chem. Phys. 131, 084104 (2009).
- [59] J. Neugebauer, J. Chem. Phys. 126, 134116 (2007).
- [60] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
- [61] M. A. Mosquera, D. Jensen, and A. Wasserman, Phys. Rev. Lett. 111, 023001 (2013).
- [62] M. Pavanello, J. Chem. Phys. 138, 204118 (2013).
- [63] O. Gunnarsson, M. Jonson, and B. I. Lundqvist, Phys. Lett. 59A, 177 (1976).
- [64] H. Eshuis, J. Yarkony, and F. Furche, J. Chem. Phys. 132, 234114 (2010).
- [65] R. Podeszwa, W. Cencek, and K. Szalewicz, J. Chem Theory Comput. 8, 1963 (2012).
- [66] J. C. Kimball, Phys. Rev. B 14, 2371 (1976).
- [67] C. König, N. Schlüter, and J. Neugebauer, J. Chem. Phys. 138, 034104 (2013).
- [68] M. A. Marques, N. T. Maitra, F. M. Nogueira, E. Gross, and A. Rubio, editors, Fundamentals of Time-Dependent Density Functional Theory, Springer, 2012.
- [69] A. J. Misquitta, B. Jeziorski, and K. Szalewicz, Phys. Rev. Lett. 91, 033201 (2003).
- [70] B. Jeziorski, R. Moszynski, and K. Szalewicz, Chem. Rev. 94, 1887 (1994).
- [71] A. Heßelmann and G. Jansen, Chem. Phys. Lett. 367, 778 (2003).
- [72] A. Heßelmann and G. Jansen, Chem. Phys. Lett. 357, 464 (2002).
- [73] J. Neugebauer, M. J. Louwerse, E. J. Baerends, and T. A. Wesolowski, J. Chem. Phys. 122, 094115 (2005).
- [74] G. te Velde, F. M. Bickelhaupt, E. J. Baerends, S. J. A. van Gisbergen, C. Fonseca Guerra, J. G. Snijders, and T. Ziegler, J. Comput. Chem. 22, 931 (2001).
- [75] C. König and J. Neugebauer, Phys. Chem. Chem. Phys. 13, 10475 (2011).
- [76] A. Götz, S. Beyhan, and L. Visscher, J. Chem. Theory Comput. 5, 3161 (2009).
- [77] K. Pernal, R. Podeszwa, K. Patkowski, and K. Szalewicz, Phys. Rev. Lett. 103, 263201 (2009).
- [78] Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998).
- [79] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- [80] A. Lembarki and H. Chermette, Phys. Rev. A 50, 5328 (1994).
- [81] Amsterdam density functional program, Theoretical Chemistry, Vrije Universiteit, Amsterdam, URL: http://www.scm.com.
- [82] J. Neugebauer, C. Curutchet, A. MunÃÂoz-Losa, and B. Mennucci, J. Chem. Theory Comput. 6, 1843 (2010).
- [83] P. Jurecka, J. Sponer, J. Cerny, and P. Hobza, Phys. Chem. Chem. Phys. 8, 1985 (2006).
- [84] A. J. Misquitta and K. Szalewicz, J. Chem. Phys. 122, 214109 (2005).
- [85] O. V. Gritsenko, P. R. T. Schipper, and E. J. Baerends, Int. J. Quantum Chem. 76, 407 (2000).
- [86] P. R. T. Schipper, O. V. Gritsenko, S. J. A. van Gisbergen, and E. J. Baerends, J. Chem. Phys. 112, 1344 (2000).
- [87] O. V. Gritsenko, P. R. T. Schipper, and E. J. Baerends, Chem. Phys. Let. 302, 199 (1999).
- [88] J. Klimes, D. R. Bowler, and A. Michaelides, J. Phys.: Condens. Matter 22, 022201 (2010).
- [89] O. A. Vydrov and T. Van Voorhis, J. Chem. Theory Comput. 8, 1929 (2012).
- [90] M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
- [91] S. Laricchia, L. A. Constantin, E. Fabiano, and F. Della Sala, J. Chem. Theory Comput. 10, 164 (2014).
- [92] C. R. Jacob, S. M. Beyhan, and L. Visscher, J. Chem. Phys. 126, 234116 (2007).
- [93] R. Kevorkyants, M. Dulak, and T. A. Wesolowski, J. Chem. Phys. 124, 024104 (2006).
- [94] M. S. Marshall, L. A. Burns, and C. D. Sherrill, J. Chem. Phys. 135, 194102 (2011).
- [95] D. Chong, Mol. Phys. 103, 749 (2005).