Subsystem-DFT potential-energy curves for weakly interacting systems †

Kohn–Sham density-functional theory (DFT) within the local-density approximation (LDA) or the generalized-gradient approximation (GGA) is known to fail for the correct description of London dispersion interactions. Often, not even bound potential-energy surfaces are obtained for van der Waals complexes, unless special correction schemes are employed. In contrast to that, there has been some evidence for the fact that subsystem-based density functional theory produces interaction energies for weakly bound systems which are superior to Kohn–Sham DFT results without dispersion corrections. This is usually attributed to an error cancellation between the approximate exchange–correlation and non-additive kinetic-energy functionals employed in subsystem DFT. Here, we investigate the accuracy of subsystem DFT for weakly interacting systems in detail, paying special attention to the shape of the potential-energy surfaces (PESs). Our test sets include the extensive S22x5 and S66x8 data sets. Our results indicate that subsystem DFT PESs strongly vary depending on the functional. LDA results are usually quite good, but behave differently from their KS counterparts. GGA results from the popular Perdew–Wang (PW91) set of functionals produce PESs that are often, but not in general overbinding. Results from Becke–Perdew (BP86) GGAs, by contrast, show the typical problems known from the corresponding KS results. We provide some preliminary results for empirical corrections for both PW91 and BP86 in subsystem DFT.


Introduction
It has been well established during the past years that Kohn-Sham density functional theory (DFT) within the local-density approximation (LDA) and typical generalized gradient approximations (GGAs), but also DFT employing hybrid exchangecorrelation (XC) functionals is problematic for the description of London dispersion forces (for reviews, see e.g., ref. 1-3).Already in 1994, Kristya ´n and Pulay have noted that ''the total exchange-correlation energy of two non-overlapping charge distributions is the sum of the individual contributions for any (semi)local DFT'', 4 and consequently, such approximate XC functionals must fail for the long-range behavior of dispersion interactions.But also the so-called medium-range dispersion at distances with weak but non-zero density overlap is usually not well described.Consequently, several correction schemes have been developed within the DFT community.Among the commonly used ones are semilocal and hybrid XC functionals which are highly parametrized to recover medium-range dispersion, 5,6 (empirical) pair-wise potentials with the desired ÀC 6 /R 6 long-range behavior (e.g., the so-called DFT-D methods), [7][8][9][10] nonlocal van der Waals (vdW) functionals, 11,12 and one-electron potentials. 13,14ubsystem DFT 15,16 (abbreviated as sDFT hereafter; for a recent review, see ref. 17) is a density-functional theory variant in which interacting subsystems are described in terms of separate Kohn-Sham-like reference systems.It is related to the frozen-density embedding (FDE) method, 18,19 which, for the purposes of this study, can be regarded as a non-selfconsistent variant of sDFT (see ref. 17 for more details).In sDFT, the nonadditive part of the kinetic energy needs to be approximated in terms of a density-dependent functional.Since this is seemingly easier for weakly overlapping electron densities than for strong chemical bonds, it is commonly assumed (and confirmed by calculations [20][21][22][23] ) that sDFT works well for weakly interacting systems.This is easy to understand for polar molecules, where the intermolecular interactions are dominated by electrostatics: these interactions are treated essentially exactly in sDFT.For non-polar molecules, however, one should actually expect similar failures as in Kohn-Sham-(KS-)DFT.Especially the long-range dispersion should suffer from the same problems.In fact, already Senatore and Subbaswamy observed in 1986 that the proper -(C 6 /R 6 ) long-range distance dependence for London dispersion interactions cannot be obtained from approximate sDFT as employed in their work. 15 quite natural idea is thus to employ DFT-D corrections also in the context of sDFT, as has been used by one of the present authors. 24This can be further motivated by the fact that sDFT and KS-DFT become equivalent in the limit of exact nonadditive kinetic-energy functionals (given the same approximation for the XC functional).Also the concept of non-local van-der-Waals functionals has been successfully transferred to the FDE context. 25The situation for medium-range dispersion in sDFT employing standard GGA-type functionals [most often of Perdew-Wang91 (PW91) type] is more complicated due to the interplay of (attractive) exchange-correlation and (repulsive) non-additive kinetic-energy components.This calls for a more detailed analysis of the behavior at medium range, which directly affects the interaction energies and bond distances of van-der-Waals complexes.In order to use sDFT for geometry optimizations or molecular dynamics simulations, equilibrium distances as well as the shape of the potential energy surfaces and accurate interaction energies are of great importance. 20,26nterestingly, sDFT has been reported to provide better interaction energies for weakly bound systems.In fact, already the Gordon-Kim model, which is a predecessor method of sDFT, produced rather accurate interaction energies for raregas dimers. 27,28Extensive early studies of interaction energies from FDE were conducted by Wesołowski and co-workers, 23,[29][30][31][32][33][34] which have been reviewed in ref. 35.One particularly interesting set of results was obtained by Wesołowski and Tran in 2003. 23hey investigated 25 interacting closed-shell systems and found that sDFT resulted in much better interaction energies than KS-DFT for GGA-type functionals.Also for LDA-type functionals, very good results were obtained in such cases where the density overlap is very small (e.g., rare gas atom dimers).][38] Here it was found that sDFT using LDA works very well (in fact often better than KS-DFT) for several classes of weakly interacting systems, with the exception of p-stacking interactions and complexes with strong charge-transfer character.Later, this study was extended by (i) investigating a broader range of kinetic-energy functionals and (ii) including additional transition-metal complex benchmark sets. 39In that study, it could be confirmed that already the LDA often offers good interaction energies, while GGA-functionals usually lead to a significant improvement.
An apparent explanation for this often superior accuracy of sDFT over KS-DFT (using the same XC-approximation) for weakly interacting complexes is a fortunate error cancellation between the non-additive XC and kinetic-energy functional approximations.But even though several studies on (interaction) energies from sDFT have appeared in recent years, 24,25,40 there has been no systematic test of potential-energy surfaces for weakly bound systems with this method.For individual complexes like an H 2 Á Á ÁNCH, 29 an F-HÁ Á ÁNCH complex, 33 benzene dimer, 32 or cytosine dimer 21 some sDFT-PESs are available in the literature.
In addition, a study addressing equilibrium geometries of the complexes from the HB6/04, DI6/04, and WI9/04 test sets 36,37 has been conducted with sDFT in ref. 20.There, it was found that LDA leads to ''excellent intermolecular equilibrium distances for hydrogen-bonded complexes'', with a maximum error of 0.13 Å for the ammonium dimer.For rare gas atom dimers, somewhat larger errors of up to 0.32 Å were obtained.GGA-type functionals hardly offered any improvement in those cases.In fact, many structures were found to be worse than with LDA.Another aspect motivating our current work is that there are almost no applications of sDFT to the popular S22 41 and S66 42 test sets and their extended versions, which include several points along a potential-energy curve; only Pavanello and co-workers have included one specific combination of GGA-type functionals in their study for the S22 test set. 25Also, many previous studies relied on equilibrium structures for weakly interacting complexes obtained from accurate reference methods such as coupled cluster theory.
Here, we fill this gap and conduct a study focussing more on the characteristic features of (entire) potential-energy curves.In particular, bond lengths and interaction energies at consistently obtained equilibrium distances will be investigated.Furthermore, we will analyze the overall shape of the potential-energy curves and indicate possible empirical corrections to further improve them.One additional comment needs to be made here: many previous studies of sDFT use PW91type 43,44 approximations as a typical GGA, while this functional is somewhat unusual in its description of weakly interacting complexes (when compared to other typical GGAs) in the KS-DFT context.This needs to be kept in mind, since many dispersion corrections in DFT address the more common behavior found, e.g., in GGA functionals like BP86 45,46 or BLYP. 45,47In order to allow comparisons for these two different types of behavior, we will include both PW91 and BP86 in our tests, in addition to LDA.
This work is organized as follows: in Section 2, we present the methodology employed here, followed by an overview over the computational details in Section 3. Section 4 contains a (re-) investigation of potential-energy curves for rare-gas dimers.In Section 5, we present single-point interactionenergy calculations for a large set of test structures at geometries taken from reference electronic-structure methods.This includes the S22 and S66 test sets.Subsequently, we investigate potential-energy curves of closed-shell dimers in Section 6, where we employ the extended S22x5 48 and S66x8 42 test sets.Finally, we outline possible strategies for the improvement of sDFT potential-energy curves, and we conclude from our results in Section 7.

Methodology
The main idea of subsystem-DFT is a partitioning of the total electron density r tot into subsystem densities r I , where N is the number of subsystems.Each of these densities will be expressed in terms of an individual Kohn-Sham-like reference system of non-interacting particles, i.e., through a set of orbitals c I i , where n I is the number of electrons in subsystem I.This is in contrast to Kohn-Sham-DFT, where the total density is expressed in terms of one single system of non-interacting particles described by orbitals c tot i , where n tot ¼ P I n I is the total number of electrons in all subsystems.The total energy expression in subsystem DFT can be formulated such that the DFT nature of the approach is emphasized, where V nn is the nucleus-nucleus repulsion energy, and V en , J, E xc , and T s are the electron-nucleus interaction, the electronelectron Coulomb, the exchange-correlation, and the singleparticle kinetic energy, respectively, as defined in the context of Kohn-Sham DFT.The non-additive kinetic energy functional is defined as, Although the exact form of T nad s is only known in terms of orbitals rather than in terms of densities, a key aspect of most practical sDFT calculations is that T nad s is evaluated using explicit density-dependent approximations.Note that formally also the exchange-correlation energy can be partitioned into subsystem contributions E xc [r I ] and a non-additive part, Realizing that all Coulomb-like energy terms naturally consist of intra-and inter-subsystem terms, we can write the total energy in sDFT in a hybrid-energy fashion, where E KS [ r I ] is the Kohn-Sham energy of subsystem I (including the nucleus-nucleus interactions), and E ˜OF-DFT int is the interaction energy between the subsystems with densities as obtained in the complex, which is evaluated in an orbital-free-(OF-) DFT manner, Here, the sums over A I and B J run over all nuclei A and B in subsystems I and J (with nuclear charges Z I A , Z J B at positions R I A , R J B ), respectively, and v I nuc (r) is the nuclear potential of the nuclei in subsystem I.Note that the ''true'' interaction energy differs from this expression, as will be discussed below.Minimization of E sDFT [ r 1 ,. ..,rN ] with respect to the electron density of system I (keeping all other subsystem densities fixed) gives rise to the so-called Kohn-Sham equations with constrained electron density (KSCED), 19 where v I s (r) is the Kohn-Sham effective potential that would appear for the isolated subsystem I.The embedding potential for subsystem I is given as In sDFT, this energy minimization is carried out iteratively for all subsystem densities, until the total energy functional is minimized.
In the strict definition of FDE, only the electron density of one single (active) subsystem is optimized in a given, fixed background (environment) electron density of all other subsystems.This means that this density does not even need to be represented in terms of (one or several) sets of Kohn-Sham orbitals as in sDFT.Usually, however, also in FDE the densities of the environment systems are represented in terms of orbitals (e.g., of isolated systems), so that also an approximate total energy can be evaluated.This corresponds to a wider definition of FDE, which coincides with an approximate sDFT (see, e.g., the preface to ref. 49 or ref. 17 for a more detailed discussion of these formal aspects).The ''true'' interaction energy E int in a complex is defined as which in the case of sDFT means where E complex (or E sDFT , respectively) is the total energy of the complex and the total energy of the isolated monomer I is given by E KS [r iso I ], where r iso I is the density of the isolated subsystem I.This will be the definition used in the following sections.It can differ from the expression given in eqn (8) for two reasons: (i) the electron densities in the complex can be different from the isolated-molecule densities, and (ii) structural relaxation may take place.The latter point will be ignored in this study, i.e., we only consider fixed monomer structures.For the case of sDFT, a relation between the two definitions of the interaction energy can be established as follows: we rewrite E int in terms of E ˜OF-DFT int from eqn (8) and of an additional term taking the change in the monomer energies into account, Here, E KS I [r I ] is the total energy of the relaxed monomer I in the complex.If we sum up all promotion energies for all subsystems, we get the total promotion energy DE prom monomers .This definition of the promotion energy is in line with the definition as given in ref. 24 (also called ''internal strain energy'' in ref. 33).The energy decomposition used here 24  For further analysis of the interaction energy, it is useful to separate E ˜OF-DFT int into different contributions, where we combined all electrostatic terms [the first three terms in eqn ( 8)] into the electrostatic interaction energy For van-der-Waals (vdW) complexes of non-polar molecules, the electron densities may actually be well approximated by the densities of the isolated molecules.Consequently, an energy evaluation neglecting the mutual relaxation of the densities under the influence of the embedding potential may still give results close to the fully relaxed case, since promotion energies are usually close to zero in those cases.

Computational details
We will provide results from three different choices of approximate functionals, where always conjoint approximations 54 for exchange and non-additive kinetic energy are employed: 1][22][23] The non-additive kinetic energy contributions are evaluated using the Thomas-Fermi (TF) model. 55,56These combinations will be denoted as KS-LDA and sDFT-LDA/TF for supermolecular KS-DFT calculations and sDFT, respectively.
(ii) the PW91-GGA for exchange, correlation, and, in reparametrized form, for the non-additive kinetic energy; 29,57 PW91k is one of the most widely used approximations for the non-additive kinetic energy in FDE and subsystem-DFT applications, and was very successfully applied for interaction energies in ref. 23.However, the corresponding XC functional is non-typical for standard KS-DFT calculations on weakly interacting complexes, as it often produces bound potential-energy curves (in contrast to many other GGAs).This leads to the abbreviations KS-PW91 and sDFT-PW91/PW91k, respectively.
(iii) the Becke88 GGA for exchange, 45 and the corresponding (reparametrized) functional for the non-additive kinetic energy by Lee, Lee, and Parr, 54 which is usually denoted as LLP91.
Here, we choose the Perdew86 46 functional for the correlation part, which leads to the common XC-GGA known as BP or BP86 with the typical failure for dispersion interactions.The corresponding calculations here will be labeled as KS-BP86 and sDFT-BP86/LLP91, respectively.
All sDFT calculations were carried out with the FDE/ subsystem DFT implementation 39,58 in the ADF 2014.01 59 program package.Usual KS-DFT calculations were also carried out with ADF 2014.01.All calculations were done without explicit use of point-group symmetry.The TZ2P 60 basis set was used in all calculations and we used the mono-molecular expansion in the case of sDFT.This expansion is free of the basis set superposition error (BSSE) and gives good results as long as no charge-transfer like interactions are important. 21,24Moreover, we want to assess the reliability of efficient and practically applicable sDFT setups suitable for production calculations, so that no super-molecular basis sets are considered.For a comparison of results from mono-molecular and super-molecular basis sets, we refer the reader to ref. 21 and 24.We have employed the pair fitting 61 (together with the standard auxiliary fit functions for the basis set) for the evaluation of Coulombpotentials and energy terms in all KS and sDFT calculations.KS-DFT interaction energies were corrected for basis set superposition errors using the counterpoise method. 62More detailed information about the calculations is provided in the ESI.† 4 Rare-gas dimers re-investigated We start our investigation with prototypical and simple test systems for van-der-Waals interactions, namely, rare-gas atom dimers.
It was shown in ref.21 that equilibrium distances for HeÁ Á ÁNe, HeÁ Á ÁAr, NeÁ Á ÁNe and NeÁ Á ÁAr derived with sDFT-LDA/TF are in good agreement with reference geometries.Also interaction energies for the above systems as well as for ArÁ Á ÁAr from sDFT-LDA/ TF were shown to be in good agreement with reference CCSD(T) ones. 20,23So far, there was no detailed comparison of entire rare-gas dimer PESs.A set of rare-gas atom dimer PESs is presented in Fig. 1.Here, we compare energies we derived from BSSEcorrected CCSD(T)/aug-cc-pVQZ (aug-cc-pVQZ-PP for Xenon; for more details see ESI †) and sDFT.For comparison, we also include the corresponding KS-DFT results.The choices for exchangecorrelation and non-additive kinetic-energy functionals are motivated in Section 3. Looking first at the KS-DFT results, we observe that KS-LDA curves are too attractive (see also, e.g., ref. 63 and 64 for analyses of KS-DFT for rare-gas dimers): they lead to shorter equilibrium distances and larger interaction energies when compared to the CCSD(T) reference.KS-PW91 leads to a bound curve with varying quality of equilibrium distances in all cases (reasonable for NeÁ Á ÁNe and KrÁ Á ÁKr, but too long for ArÁ Á ÁAr and XeÁ Á ÁXe), while KS-BP86 leads to the well-known 65 behavior of purely repulsive potential-energy curves.The sDFT interaction energies follow the overall trends of their KS-DFT counterparts, however, with some differences.The sDFT-LDA/TF potential-energy curves become less attractive and equilibrium distances are improved in most cases when compared to KS-LDA.sDFT-PW91/PW91k potential-energy curves are more attractive than their KS-PW91 counterparts, which also degrades the description of equilibrium distances.The sDFT-BP86/LLP91 curves are purely repulsive so that this combination of functionals cannot be used to describe rare-gas dimers adequately.
It becomes obvious that the quality of sDFT results for weak interactions strongly depends on the choice of the functionals.In fact, the potential-energy curves of sDFT-BP86/LLP91 are just as poor as KS-BP86.The bound character of the KS-DFT curves for LDA and PW91 is recovered in the sDFT calculations.sDFT-LDA/TF binding energies are lowered, while sDFT-PW91/ PW91k binding energies are increased compared to their KS counterparts.This also leads to larger equilibrium distances in sDFT-LDA/TF and smaller ones in sDFT-PW91/PW91k when compared to KS-DFT.At CCSD(T) equilibrium distances, the interaction energies from sDFT-LDA/TF are rather good.But this does not automatically imply the same accuracy for equilibrium distances.The sDFT-LDA/TF curves nicely follow the CCSD(T) ones at medium and larger distances.This observation is in line with the early results by Gordon and Kim, 28 who used an approach similar to sDFT-LDA/TF with unrelaxed densities.Although the densities in the present calculations are fully relaxed (and not derived from Hartree-Fock theory) the behavior of the sDFT-LDA/TF interaction energies is very similar to those presented in ref. 28.

Single-point interaction energies
As shown above, sDFT-LDA/TF interaction energies for rare-gas dimers are usually quite accurate, while those derived from sDFT-PW91/PW91k are generally too attractive, and sDFT-BP86/ LLP91 does not describe bound states.In order to further investigate the behavior for chemically more relevant systems, we re-investigate the test set employed by Wesołowski and Tran. 23ur study is based on a subsystem DFT energy implementation 39 in the Slater-basis code ADF, 59 while their previous study made use of a Gauss-basis implementation 33 in the DEMON code. 66,67hese results thus also serve the purpose of additional validation of the ADF implementation and a comparison of the two types of basis sets.We then continue with an application of sDFT to the popular S22- 41,68 and S66- 42,68 test sets.The structures investigated in the following are fixed structures from the reference articles, which have not been re-optimized in the present work unless noted otherwise.

The test set by Wesołowski and Tran
Our first test set is the set of 25 intermolecular complexes studied in ref. 23 with interaction energies between À0.3 and À12.0 kJ mol À1 .The structures were derived as follows: the monomer structures for complexes 1, 7, 8, 10, 14, 16, 19-25 in Table 1 were taken from references within ref. 23 and their respective supporting information.For all other complexes, the monomers were optimized with MP2/6-31G* in the original references, but coordinates were not reported.In order to get as close to the reference monomer structures as possible, we used the same method to derive optimized monomer structures using the GAUSSIAN 09 69 program package.The distances between the monomers were then taken from the respective references within ref. 23.Because the latter deviation seems quite large, we checked this interaction energy against results from ORCA 3.0.2 70ith RI-DFT/cc-pVTZ, 71 which supports our findings with an interaction energy of 3.26 kJ mol À1 (compared to our value of 3.08 kJ mol À1 ).The magnitude of the interaction is similar in ref. 23, but the sign is different (À3.05 kJ mol À1 ). Incontrast to ref. 23, we observe that the average absolute error of KS-LDA with respect to the wave-function reference results is smaller than that of sDFT-LDA/TF, which can most likely be attributed to the basis set.For this test set, KS-LDA is always in the regime of chemical accuracy (B1 kcal mol À1 ) and usually slightly too attractive.sDFT-LDA/TF, by contrast, yields interaction energies that are almost always less attractive than the reference wave-function data, and always less attractive than their KS-DFT counterparts (see Table 2).This situation is reversed for PW91: while KS-PW91 mostly underestimates interaction energies, sDFT-PW91/PW91k usually overestimates them.
For KS-BP86 and sDFT-BP86/LLP91, the interaction energies are very similar: the mean and mean absolute deviations (MD and MAD, respectively) between sDFT-BP86/LLP91 and KS-BP86 interaction energies are 1.06 and 1.14 kJ mol À1 , respectively.Hence, sDFT-BP86/LLP91 shows the best agreement with the corresponding KS-DFT results, which also means that it shows the typical problems of KS-BP86 in the description of weak interactions.But this, in turn, implies that empirical dispersion corrections developed for KS-BP86 could easily be transferred to its sDFT counterpart (see Section 6.3).

The S22 and S66 test sets
We continue our investigation with the S22 41,68 and S66 42,68 test sets.][74] In both sets, the complexes are grouped into three subsets in ref. 41 and 42, based on findings from symmetry adapted perturbation theory calculations.We use the same classification scheme here.The S22 and S66 complexes are explicitly listed for the three different categories in Tables 3 and 4. Both sets contain only the equilibrium structures (obtained with CCSD(T) and MP2) of the respective complexes.Because the interaction energies for these reference structures are similar for both sets, they are discussed together here.The results are summarized in Fig. 3 and 4.
From these figures, we see that KS-LDA interaction energies are usually too attractive, while this is very rarely the case for sDFT-LDA/TF compared to the reference.Compared to their KS-DFT counterparts, sDFT-LDA/TF interaction energies are commonly less attractive, which is in line with the findings from Section 5.1.The MDs of sDFT-LDA/TF from KS-LDA are quite large (cf.Table 5), but the MADs of both methods compared to the reference are similar (8.65 vs. 9.69 kJ mol À1 with KS-LDA and sDFT-LDA/TF for S22, respectively and 7.07 vs. 7.27 kJ mol À1 respectively for S66).sDFT-LDA/TF interaction energies are quite poor when p-interactions play a role, which is also in agreement with previous results from the literature. 20,21,23owever, p-p-stacked interactions are quite well described by KS-LDA (also in line with findings from ref. 75).KS-LDA quite nicely fits to the reference interaction energies of the vdW and mixed bonded complexes.
Interaction energies from KS-BP86 and sDFT-BP86/LLP91 are most often too repulsive compared to the reference for these test sets.The MADs between sDFT and KS-DFT are similar for BP86 and PW91 (cf.Table 5).For both approximations, the sDFT interaction energies are mostly more attractive than the KS-DFT ones for vdW and mixed interactions.For hydrogenbonded systems, sDFT-BP86/LLP91 is usually less attractive than its KS-BP86 counterpart.None of the sDFT approximations show a general overestimation of dispersion-dominated interaction energies, while simple KS-LDA does in almost all  cases.Interaction energies from sDFT-LDA/TF and sDFT-BP86/ LLP91 show similar MADs from the reference compared to their KS-DFT counterparts.Interaction energies from KS-PW91 and KS-BP86 are usually not attractive enough when compared to the reference.Since interaction energies from sDFT-PW91/ PW91k are quite systematically more attractive than KS-PW91 interaction energies, the MAD of sDFT-PW91/PW91k interaction energies is the best for these test sets (4.72 and 2.90 kJ mol À1 for S22 and S66, respectively).
6 Potential-energy curves from FDE

S22x5 test set
We now discuss systematics in PESs derived from sDFT and start with the S22x5 48,68 test set, which consists of the complexes listed in Table 3.In addition to the equilibrium structures, this set contains four extra structures displaced from equilibrium for each complex in the set.In the displaced structures, the intermolecular distance is varied, so that the S22x5 set concerns representative parts of the intermolecular potential-energy curves.The MADs for each set of sDFT calculations for the S22x5 test set are presented in Fig. 5 along with their KS-DFT counterparts.Not just for equilibrium structures, but also for PESs, the average MAD of sDFT-PW91/PW91k with a magnitude of 2.92 kJ mol À1 is the smallest among the DFT  The weaker the interactions between the complexes are (i.e., at large intermolecular distances), the less pronounced should be the effect of the non-additive contributions to the embedding potential and to the interaction energy.Hence, the results of KS-DFT and sDFT converge for larger distances.This can clearly be seen for all curves in this test set.In order to clarify the origin of binding in sDFT for weakly interacting systems, a decomposition of the sDFT interaction energy contributions according to eqn ( 12) and ( 13) for the three different sDFT approximations used here is presented for the ammonia and methane dimer complexes in Fig. 7 and 8, respectively.
Around a center of mass distance of 3.7 Å and shorter distances for the ammonia dimer, the non-additive exchange-correlation energy E nad xc has values comparable to the total electrostatic interaction in the sDFT-LDA/TF case.The magnitude of the non-additive kinetic energy T nad s is similar to E nad xc for all but the shortest distances.But T nad s is of opposite sign.Therefore, all non-additive contributions play a non-negligible role in the description of the interaction of the two subsystems.The nonadditive terms more or less cancel for center-of-mass distances of 3.5 Å and larger.Hence, the total interaction energy appears to be almost identical to the total electrostatic interaction energy.This journal is © the Owner Societies 2015 A similar picture is observed for sDFT-PW91/PW91k, while sDFT-BP86/LLP91 shows a smaller (in magnitude) XC contribution, thus leading to a more repulsive total interaction energy curve.
For the methane dimer (Fig. 8) it is clearly demonstrated that binding occurs mainly because of the non-additive exchange-correlation parts in sDFT-LDA/TF and sDFT-PW91/ PW91k.The non-additive exchange-correlation energy is even positive for sDFT-BP86/LLP91 for distances larger than 3.75 Å.The magnitude of the non-additive contributions is clearly larger than that of the electrostatic interaction energy, as one could expect for these non-polar molecules.In both examples presented here, the polarization of the subsystems plays a minor role for the interaction energy.Depending on the polarity and polarizability of the monomers in the complex, this is not generally the case for all intermolecular complexes.For the equilibrium structures of the S22x5 test set (which is, in fact, the S22 test set) a maximum monomer polarization energy of E pol monomer E 50 kJ mol À1 can be observed for the formic acid dimer.‡ Thus, monomer polarization energies can be of significant magnitude.For the S22 and S66 test sets this is typically the case for hydrogen-bonded complexes.All promotion energies for the S22 and S66 test sets can be found in the ESI † (Tables 3 and 4).The overall average MAD of sDFT-PW91/PW91k is the best for this benchmark set (2.92 kJ mol À1 ), but taking only the second and third block into account (complexes 8-22), KS-LDA shows, on average, even smaller deviations with respect to the reference.The average MADs of KS-LDA and sDFT-LDA/TF are almost equivalent (5.95 and 5.98 kJ mol À1 , respectively) although their descriptions of the PESs differ largely from each other.As explained before, for large distances the KS-DFT and sDFT interaction energies coincide.Hence, the MADs of both methods are very similar, because most of the data points belong to this region.KS-LDA is mostly too attractive for the complexes in this set, so that equilibrium distances are generally too short.For complexes 11-15, which are all p-p stacked complexes, KS-LDA interaction energies across the whole PES are usually very good.This is in line with the findings in ref. 75.For the interaction energies at equilibrium distances (see Section 5) we have seen that sDFT-LDA/TF is usually less attractive than KS-LDA.This picture is reversed for PW91 and BP86, where PESs from sDFT are usually more attractive than their KS-DFT counterparts.This usually remains true for the whole PES (see Fig. 6), except for cases where sDFT and KS-DFT are very similar anyway [cf.Fig. 6(e)].For hydrogen-bonded complexes, sDFT-BP86/ LLP91 is most often less attractive than KS-BP86 for the whole PES.While KS-PW91 usually does not describe vdW complexes sufficiently well, interaction energies of sDFT-PW91/PW91k are better at equilibrium distances of the reference method.But the repulsive part of the PESs at shorter distances is usually poorly described.For the first group of molecules in the S22x5 test set (hydrogen-bonded complexes), sDFT-PW91/PW91k is more attractive than KS-PW91 when the X-H bond (where X can be N or O) is pointing directly towards the donating lone pair of the other monomer.In the second group it can be seen that sDFT-PW91/PW91k interaction energies are not attractive enough to fit the reference, but usually the functional form is qualitatively acceptable.At short distances, however, interaction energies are usually more attractive than the reference.Regarding the last group, consisting mainly of X-HÁ Á Áp-interactions, sDFT-PW91/PW91k interaction energies are still more attractive than KS-PW91 and the error with respect to the reference is typically in the range of chemical accuracy for the whole PES.

S66x8 test set
We extend our investigation on PESs from sDFT with the larger S66x8 42,68 test set, which consists of equilibrium structures of the complexes described in Table 4 plus 7 displaced structures for each complex in the set.The MADs are visualized in Fig. 9.The complete set of figures for all PESs together with the values corresponding to Fig. 9 are provided in the ESI.† Here, we choose to provide examples for six complexes for each of the three categories (hydrogen bonding, vdW, mixed interactions).The first set is presented in Fig. 10 for hydrogen-bonded complexes.
As we have seen before, KS-LDA is too attractive (with respect to the reference) for all hydrogen-bonded complexes.Moreover, interaction energies derived from its sDFT counterpart are in general more repulsive than the reference.Therefore equilibrium distances from sDFT-LDA/TF are usually larger than those from KS-LDA, but also fit better to the reference.KS-PW91 interaction energies are very similar to those from sDFT-PW91/PW91k for systems involving one hydrogen bond between two hydroxy-groups.If an N-H-group is the acceptor of a hydrogen bond from an electron donating group (e.g.AcNH 2 dimer), the differences become larger and sDFT-PW91/PW91k is usually more attractive than KS-PW91 (and mostly also the reference) around distances smaller than the reference equilibrium.The same behavior can be recognized for sDFT-BP86/ LLP91 and its KS-DFT counterpart.If N-H-groups take part as the accepting partner in a hydrogen bond, sDFT-BP86/LLP91 interaction energies are close to KS-BP86, but at smaller distances they tend to be more attractive.However, BP86 in any formulation is not able to recover reference interaction energies, but equilibrium distances are modeled quite well for hydrogenbonded complexes in this set.
The second set, containing six examples from vdW-bonded complexes, is shown in Fig. 11.Again KS-LDA describes potential energy surfaces of p-p-stacked complexes quite accurately, while equilibrium distances are often too small compared to the reference.For larger distances KS-LDA is always quite accurate for this group of vdW-dominated interactions (also due to the smaller magnitude of interaction energies).With sDFT-LDA/TF the bonding character of KS-LDA is kept, but interaction energies are more repulsive than the reference.Equilibrium distances, however, are often quite accurate.KS-PW91 gives bounding curves, but they are usually too repulsive around reference equilibrium geometries.Hence, equilibrium distances derived from KS-PW91 are usually too large.Interaction energies derived from sDFT-PW91/PW91k for p-p-stacked interactions are often close to the reference, but the repulsive parts of the PESs are again too attractive.Thus, equilibrium distances are usually too small.KS-BP86 results in repulsive curves without bonding character, except for some rare cases where equilibrium distances are very poor.sDFT-BP86/LLP91 follows the trends of its KS-DFT counterpart, but tends to be more attractive around reference equilibrium geometries.
The six examples of the last set consisting of mixed-type interactions are presented in Fig. 12.The trends discussed above for LDA and PW91 also hold for this set of PESs, so that sDFT-LDA/TF often describes reference equilibrium distances better than its KS-DFT counterpart.Interaction energies from sDFT-LDA/TF are again mostly not attractive enough, while KS-LDA is much too attractive at short distances (always compared to the reference).KS-PW91 is in general too repulsive.However, interaction energies as well as equilibrium distances for the ethyne-involving complexes in this group are very often quite accurate.Also for this set, sDFT-PW91/PW91k is systematically more attractive than its KS-DFT counterpart, which leads to equilibrium distances that are usually shorter than the reference.As we have seen in Section 5, interaction energies from sDFT-PW91/PW91k can be quite accurate at reference equilibrium distances.Nonetheless, the qualitative picture is usually not correct, because the repulsive part is not described correctly.KS-BP86 shows bounding curves for complexes where electrostatic interactions are more pronounced.Apart from that, interaction energies from KS-BP86 are usually too repulsive and equilibrium distances are often too large.In those cases where equilibrium distances from KS-BP86 nicely fit the reference, PESs derived from sDFT-BP86/LLP91 are too attractive at shorter distances.This trend for mixed-bonded interactions is similar to the one observed for PW91 when comparing KS-DFT and sDFT.Here this is mostly the case if some kind of hydrogen-bonded interaction without proper orientation of the donating lone-pair towards the accepting X-H-bond is present.

Applying DFT-D3(BJ) to sDFT
Already in previous work, the idea of improving sDFT interaction energies by empirical DFT-D corrections was explored by one of the present authors. 24In that work, the DFT-D2 scheme 76 and the sDFT-BLYP/TW02 45,47,77 approximations were used, but only reference equilibrium structures were considered.For PESs, we have seen in the previous sections that sDFT-BP86/LLP91 results are often very similar to their KS-DFT counterparts, and thus show the typical weaknesses of non-dispersion corrected KS-DFT approximations.Hence, it can be anticipated that the DFT-D scheme should lead to similar improvements in KS-BP86 and in sDFT-BP86/LLP91.To test this assumption, we apply the latest version of the DFT-D scheme, DFT-D3(BJ), 10,78 in its original parametrization for KS-BP86 to sDFT-BP86/LLP91 interaction energies.The DFT-D3(BJ) correction was calculated using the dftd3 program (V3.1 Rev 0). 79As most of the data points in the S22x5 and S66x8 test sets describe large distances, the average MAD with sDFT-BP86/LLP91-D3(BJ) gets systematically lower.The average MAD drops down from 9.68 to 3.84 kJ mol À1 and from 10.32 to 3.10 kJ mol À1 for the S22x5 and S66x8 test sets, respectively.Two examples are shown in Fig. 13 to illustrate the effect of DFT-D3(BJ) for sDFT-BP86/LLP91 for a hydrogen-bonded and a vdW-bonded complex.The complete set of sDFT-BP86/ LLP91-D3(BJ) PESs is provided in the ESI.† Clearly, interaction energies at larger distances are improved significantly.But in cases where the description of the PES with sDFT-BP86/LLP91 at short distances deviates from KS-BP86, the DFT-D3(BJ) correction usually further worsens the result since it adds additional attractive contributions.Therefore, the description of the PESs is qualitatively only improved for large distances.An adaption of the damping function for short distances is thus needed for further improvement of sDFT-BP86/LLP91-D3(BJ) interaction energies.This journal is © the Owner Societies 2015

Repulsive empirical corrections for sDFT
Contrary to the behavior of sDFT-BP86/LLP91, we have seen that sDFT-PW91/PW91k interaction energies are often very good at reference equilibrium distances.But they are usually more attractive than the reference at shorter distances for the examples in the S22x5 and S66x8 test sets.In those cases, equilibrium distances tend to be too short.Here, we provide a first test for a possible empirical correction of interaction energies (and consequently also equilibrium bond distances) of the widely used sDFT-PW91/ PW91k approximation.Inspired by the corresponding expressions for dispersion corrections in the DFT-D methods, 2 we model the difference and sDFT-PW91/PW91k (E sDFT-PW91/PW91k int ) interaction energies [defined according to eqn (12)] with an empirical repulsive energy term.We decide to use a Buckingham-like exponential, which is motivated by the fact that the repulsive correction should be related to the overlap of the subsystem densities.Preliminary studies also revealed that the repulsive 1=R 12 component of a Lennard-Jones potential is not suited here.Finally, the atom-pair dependent energy correction G is given by Here, A and B refer to atoms of subsystems 1 and 2, respectively and a IJ and b IJ are empirical atom-pair (IJ) specific parameters.
The function has a value of one if the internuclear distance equals zero, and it decays smoothly to zero for larger distances.The parameter R AB 0 is a pre-defined cut-off radius for atom pair AB taken from ref. 10.The parameter s r,n is a radius scaling factor (in principle functional dependent) and g is a parameter that defines the steepness of the function.For fitting our test parameters here, we only used those complexes from the S22x5 test set for which sDFT-PW91/PW91k was clearly overbinding, i.e., where interaction energies at equilibrium distances were lower than the reference (complexes 1, 2, 4, 6-9 and 16, cf.Table 3 and Fig. 3).The model function G [eqn ( 14 The final results from the fit are summarized in Table 6.Interaction energies from sDFT-PW91/PW91k-R are by construction more repulsive than those from sDFT-PW91/ PW91k.For the chosen training set of 8 complexes (vide supra) the MAD of the interaction energies for the respective PES was significantly improved (see Fig. 14), showing that the form of the correction is adequate for these cases.Also for the other complexes in the S22x5 test set, equilibrium distances as well as the qualitative picture of the PESs are often improved (cf.Fig. 15).Since we only add repulsive corrections, however, interaction energies often get worse.In conclusion, this fit is not a general improvement, although it can improve the description of equilibrium distances and the general shapes of the PESs.

Discussion and conclusions
In this work, we have investigated sDFT interaction energies for three different conjoint approximations 54 (sDFT-LDA/TF, sDFT-PW91/PW91k, sDFT-BP86/LLP91) for several test sets including the S22 41 and S66 42 sets.We have confirmed that sDFT interaction energies from sDFT-LDA/TF and sDFT-PW91/ PW91k are often better than their KS-DFT counterparts at reference (i.e., high-level wave-function) equilibrium geometries.But they are not generally better in reproducing equilibrium distances or shapes of PESs.We then investigated entire potential-energy curves from sDFT for 5 rare-gas dimer systems as well as for the S22x5 48 and S66x8 42 sets.The shapes of the PESs strongly depend on the choice of approximations used for E xc /E nad xc and T nad s .The bonding character of KS-DFT PESs (if existent for the respective complex) is usually reproduced by the sDFT counterpart.We have seen that sDFT-BP86/ LLP91 shows the best agreement with the corresponding KS-DFT results.This similarity implies that the sDFT-BP86/ LLP91 interaction energies share the weaknesses (especially for vdW-bonded complexes) with their KS-DFT counterpart, but this can easily be remedied by empirical DFT-D corrections.In the recent literature [21][22][23]31,33 sDFT-PW91/PW91k is often mentioned for its good interaction energies at reference equilibrium geometries. This as been confirmed here.But equilibrium geometries from sDFT-PW91/PW91k are not necessarily of the same quality, since potential-energy curves are often not repulsive enough at short distances.Hence, we briefly explored the idea of using an empirical repulsive correction in the spirit of DFT-D.Our simple parametrization here demonstrates that this can improve bond distances and shapes of PESs from sDFT-PW91/PW91k.But since the error in sDFT-PW91/PW91k interaction energies does not show a uniform behaviour for different complexes, the correction in its present form often leads to worse interaction energies.Further work along these lines is clearly needed for high-quality PESs from sDFT.Finally, it should be pointed out that the results obtained here indicate only a minor dependence of the interaction energy on the embedding potential for vdW-bonded (and mixedbonded) complexes.This can be seen from the fact that the promotion energy components of the interaction energies are often small (see the examples in Fig. 7 and 8 and the complete list of DE prom monomer in the ESI †).The differences in the embedding potentials for different sDFT approximations presented here are usually very small even for polar complexes, when judging from the differences caused in DE prom monomer : maximum differences in DE prom monomer occur between sDFT-LDA/TF and sDFT-BP86/LLP91 with 4.46 and 5.51 kJ mol À1 for the S22 (formic acid dimer) and S66 (acetic acid dimer) test sets, respectively.
Overall, we conclude that sDFT certainly has the potential to produce accurate PESs, if similar effort as in KS-DFT is spent for optimizing the corresponding energy contributions.
arises naturally from the different energy contributions in sDFT.It shows some similarities (e.g. in the form of the electrostatic terms) to the energy decomposition analyses (EDAs) developed by Morokuma 50,51 and by Ziegler and Rauk 52 (for a review, see ref. 53).The promotion energy arises because of mutual relaxation in the monomer densities due to electrostatic and short-range quantum-mechanical effects upon forming the complex.It is, by definition, always positive.
We define the error DE of the DFT interaction energy E DFT int with respect to the reference interaction energyE ref int as DE = E ref int À E DFT int .Hence, positive values represent a too negative interaction energy, i.e., an overestimation of the interaction.In view of the different types of basis functions used, the results obtained with ADF are in good agreement with data presented in ref. 23.The errors DE are graphically presented in Fig. 2. The sDFT values have a maximum absolute deviation of 2.21 (C 6 H 6 Á Á ÁC 6 H 6 ) and 1.18 kJ mol À1 (C 2 H 4 Á Á ÁC 2 H 4 ) for sDFT-LDA/TF and sDFT-PW91/PW91k, respectively, when compared to the results from ref. 23.The mean absolute deviations (MADs) relative to ref. 23 for sDFT-LDA/TF and sDFT-PW91/PW91k are 0.49 and 0.41 kJ mol À1 , respectively and 0.70 and 0.61 kJ mol À1 for KS-LDA and KS-PW91.Maximum absolute deviations are 3.05 and 6.13 kJ mol À1 (both C 6 H 6 Á Á ÁC 6 H 6 ) for KS-LDA and KS-PW91, respectively.

Fig. 2
Fig. 2 Errors between sDFT-and reference interaction energies from ref. 23 for the test set by Wesołowski and Tran.Errors are defined as D E = E ref int À E DFT int .Errors from KS-DFT have been included for comparison.All values are given in kJ mol À1 .

Fig. 3
Fig. 3 Errors between sDFT-and reference CCSD(T)/CBS interaction energies from ref. 68 for the S22 test set.Errors are defined as DE = E ref int À E DFT int .Errors from KS-DFT have been included for comparison.All values are given in kJ mol À1 .

Fig. 4
Fig. 4 Errors between sDFT-and reference CCSD(T)/CBS interaction energies from ref. 68 for the S66 test set.The different pictures describe (a) hydrogen-bonded complexes, (b) vdW complexes and (c) mixed complexes.Errors are defined as DE = E ref int À E DFT int .Errors from KS-DFT have been included for comparison.All values are given in kJ mol À1 .

Fig. 5
Fig. 5 MADs per molecule of sDFT interaction energies with respect to reference CCSD(T)/CBS energies taken from ref. 68 for the S22x5 test set.MADs from KS-DFT have been included for comparison.
)] without f damp was fitted to DE int using the locally modified Levenberg-Marquardt algorithm from ref. 80. Afterwards, parameters for the damping function f damp were obtained by fitting the complete model function G to DE int while fixing the optimized parameters a IJ and b IJ , resulting in s r,n = 0.51 and g = 6.00.Optimized parameters with negative sign were discarded in order not to add artificial attractive terms.The corrected version of sDFT-PW91/PW91k interaction energies, will be termed sDFT-PW91/PW91k-R in the following.It is calculated as E sDFT-PW91/PW91k-R int = E sDFT-PW91/PW91k int + G.

Fig. 14 Fig. 15
Fig. 14 MADs of interaction energies per molecule for the empirically repulsion-corrected version sDFT-PW91/PW91k-R in comparison to noncorrected sDFT-PW91/PW91k results for the S22x5 test set.All values given in kJ mol À1 .

Table 1
Errors between sDFT-and reference [CCSD(T) or MP2] interaction energies taken from ref. 23.Errors are defined as DE = E ref int À E DFT int , where E DFT int refers to the interaction energy derived from KS-DFT or sDFT.Errors from KS-DFT have been included for comparison.All values are given in kJ mol À1 .MD and MAD stand for mean deviation and mean absolute deviation, respectively This journal is © the Owner Societies 2015 Phys.Chem.Chem.Phys., 2015, 17, 14323--14341 | 14329

Table 2
Mean deviations (MDs) and mean absolute deviations (MADs) of sDFT interaction-energy errors with respect to their KS-DFT counterparts for the Wesołowski-Tran test set.The error DE = E KS-DFT This journal is © the Owner Societies 2015

Table 4
Complexes in the S66 test set.The ''peptide'' is N-methylacetamide.(TS) stands for T-shaped

Table 5
MDs and MADs of sDFT interaction-energy errors with respect to their KS-DFT counterparts for the S22 and S66 test sets.The error DE DFT =