Tuned range-separated hybrid functionals in the symmetry-adapted perturbation theory

Symmetry-adapted perturbation theory (SAPT) in its symmetrized Rayleigh-Schrödinger (SRS) formulation has become a successful method of calculating interaction energies in non-covalent molecular complexes.^1–5 In the most popular variant of SAPT, the description of monomers is provided by density functional approximations (DFAs). This approach is denoted as DFT-SAPT or SAPT(DFT).^6,7 For the latest review of the method, see the paper by Jansen.⁵ 

As has been recognized in the early days of the method, the proper asymptotic behavior of the exchange-correlation (xc) potential, v_xc, is crucial for DFT-SAPT.^8–11 In the case of the functionals that do not exhibit derivative discontinuities at integer numbers of electrons, so called continuum functionals, the behavior of the xc potential for large distances from the nuclei should be:^12,13

where Δ_∞ is the limit value of the potential at infinity and equals IP + ɛ_H [IP stands for the lowest ionization potential and ɛ_H is the Kohn-Sham (KS) eigenvalue of the highest occupied molecular orbital (HOMO)]. A solution that fulfills this condition is the asymptotic correction (AC) scheme.^13–15 Asymptotically corrected potentials provide an improved treatment of electron density and response properties essential for DFT-SAPT applications.^7,9–11 The AC method is based on combining a semilocal potential in the bulk region with a potential exhibiting proper long-range behavior in the asymptotic region of the density. The procedure referred to as “cutting and splicing” may employ different continuous switching functions. Misquitta and Szalewicz combined the long-range Fermi-Amaldi expression with the connection scheme of Tozer and Handy.^11,13 Hesselmann and Jansen merged the LB94 potential from van Leeuwen and Baerends¹⁶ with the PBE0^17,18 functional via the gradient-regulated seamless connection method of Grüning et al.^9,15 Finally, Misquitta proposed the use of the LB94 potential switched on by the function of Tozer and Handy.¹⁹ All three flavors of the AC have recently been validated regarding their performance for molecular properties, energies of excitations to Rydberg states, and interaction energies in the paper of Cencek and Szalewicz.^20,21 In this work, as an alternative to ensuring the positive value at infinity (Eq. (1)), we employ the AC that shifts down the bulk part of v_xc by Δ_∞, while keeping the zero limit at infinity.

Two problems in applications of asymptotically corrected potentials should be mentioned. First, AC potentials belong to the class of the so-called “stray potentials,” i.e., they are not functional derivatives of any density functional of the xc energy.^22,23 This prevents a unique definition of the exchange-correlation energy. Furthermore, energies from stray potentials lack the rotational and translational invariance.^23,24 Arguably, this is not a difficulty because monomer energies are not needed in SAPT. However, calculation of energy derivatives over the field or over fractional electron numbers is impossible and this, as we will argue, hinders checking the quality of the monomer description. Second, although asymptotically corrected potentials display the exact long-range behavior, the effects of many-electron self-interaction error (MSIE)^25,26 are still present if the AC scheme involves semilocal LDA and GGA functionals.²⁷ This manifests itself in erroneous dissociation of symmetric ion radicals²⁷ and improper response to externally applied electric field in extended molecular systems.²⁸ The latter has been traced to the absence of the field-counteracting term^29–31 in the exchange potential of the underlying functional.^32–35 

The long-range corrected (LRC) functionals offer an alternative way of improving the asymptotic distance dependence^36–41 of v_xc. Moreover, the tuning of the range-separation parameter assures the proper limiting behavior of the electron density.^42,43 LRC functionals provide the xc potential in a straightforward manner, i.e., as a functional derivative v_xc = δE_xc[ρ]/δρ, avoiding inconsistencies stray potentials are burdened with. Moreover, they allow for the minimization of MSIE,^44,45 and, as a consequence, reliable characterization of long-range-dependent properties including polarizabilities and hyperpolarizabilities of challenging π-conjugated systems.^36,44,46,47

The applicability of long-range corrected functionals to DFT-SAPT was studied in two recent papers. Cencek and Szalewicz studied the performance of LRC-SAPT against a set of four reference systems for which benchmark first- and second-order energy contributions were available.²⁰ In the Erratum²¹ to this paper, they show that the optimally tuned LRC functional affords accurate density for He. They further conclude that LRC functionals with optimally tuned range-separation parameter “perform not much worse than the best AC corrected functionals” and that further investigations are warranted.

Lao and Herbert conducted a more detailed study of the LRC-SAPT method.⁴⁸ Not only did they cover both S22⁴⁹ and S66⁵⁰ data sets for non-covalent interactions, but also constructed a new SS41 data set of benchmark energy components. In their approach, the second-order induction energy was treated at the uncoupled level of theory with only approximate treatment of orbital relaxation. Since the dispersion term calculated at the uncoupled level lead to unsatisfactory results, the authors replaced it with an empirical dispersion potential, the so-called SAPT(KS)+D approach.^51,52 Their conclusion was that the tuned LRC functionals provide quantitative predictions in DFT-SAPT, competitive with those obtained with the AC models. Furthermore, LRC-SAPT was in their view a step towards introduction of analytic gradients.⁴⁸ 

The aim of this work is to present a performance test of LRC-SAPT with second-order energy components calculated at the coupled level of theory, i.e., via frequency-dependent density susceptibilities (FDDSs) from time-dependent DFT, on a wider set of dimers than in Ref. 20. We also provide a comparison of SAPT based on different DFT approaches and relate their efficacy to the features of Kohn-Sham (KS) and generalized Kohn-Sham (GKS) theory.^53,54

The structure of this paper is as follows. In Secs. III A–III C, we reveal the weaknesses of the asymptotic correction scheme which are relevant to SAPT applications. In order to exemplify the effect of the missing field-counteracting term, we calculate static polarizabilities of polyacetylenes, similarly to Sekino et al.²⁸ In Sec. III C, on another model system, |$\rm {He_3^+}$|$He 3+$⁠, we demonstrate that when the electron distribution is examined, |$v_{\rm xc}^{\rm AC}$|$v xc AC$ may lead to the identically flawed delocalization of electron density as the underlying bulk semilocal xc functional.

Finally, in Sec. III D, a comparison with other DFT-SAPT schemes, namely, SAPT(PBE), SAPT(PBEAC), and SAPT(PBE0AC) is provided. In order to complement previous works on LRC-SAPT, we focus on total interaction energies for two sets of small van der Waals dimers, the A24 set of Hobza et al.⁵⁵ and DI6-04 set of Truhlar.⁵⁶ In particular, we address the behavior of the dispersion energy and recognize the possible limitations to the LRC-SAPT approach.

In DFT-SAPT, we calculate the interaction energy, E_int, as the sum of electrostatic (E_elst), first-order exchange |$\big(E^{(1)}_{\rm exch}\big)$|$E exch (1)$⁠, second-order induction |$\big(E^{(2)}_{\rm ind}\big)$|$E ind (2)$⁠, and dispersion |$\big(E^{(2)}_{\rm disp}\big)$|$E disp (2)$ energies as well as their exchange counterparts, exchange-induction |$\big(E^{(2)}_{\rm exch-ind}\big)$|$E exch − ind (2)$ and exchange-dispersion |$\big(E^{(2)}_{\rm exch-disp}\big)$|$E exch − disp (2)$⁠:

The last term, the so called δ Hartree-Fock term, approximately accounts for higher-order terms.⁵⁷ In our implementation, the second-order induction energy is computed using the coupled KS static response theory, and the dispersion energy via FDDSs from time-dependent DFT.

We chose a LRC functional composed of the PBE correlation and the short-range (SR) exchange ωPBE based on Henderson, Janesko, and Scuseria model of the exchange hole, LRC-ωPBE.⁴¹ We support such selection by noting that this is a range-separated variant of the PBE functional which is thoroughly tested in the context of DFT-SAPT. The following range-separation scheme into short-range and long-range of the form is applied:

where erf and erfc stand for the error function and its complement, and ω is the range separation parameter. A further decomposition of SR exchange into Hartree-Fock (HF) and DFT parts is employed,

In Sec. III D, we applied functionals with two different α values: 0.0 and 0.2, denoted as LRC-ωPBE and LRC-ωPBEh, respectively. In Sec. III E, α values from 0.0 to 0.6 were applied. The range-separation parameter, ω, was optimally tuned according to the IP-tuning procedure, i.e., by minimizing the difference between the negative of HOMO energy and the vertical ionization potential (IP):⁵⁸ 

where E denotes the self-consistent energy for N and N − 1 electron systems and ɛ_H stands for energy of the HOMO orbital. The tuning procedure was carried out separately for each monomer in its dimer geometry.

All SAPT calculations were performed in the developer's version of the MOLPRO package.⁵⁹ Our implementation utilized a hybrid xc kernel composed of the short-range LDA exchange kernel scaled by (1 − α), the LDA correlation kernel and the coupled Hartree-Fock kernel scaled by α. For both the A24 and DI6-04 data sets, the aug-cc-pVTZ basis set of Dunning et al. was applied.⁶⁰ The asymptotic correction scheme of Grüning et al. was applied to the PBE¹⁷ and hybrid PBE0 functionals, yielding PBEAC and PBE0AC functionals, respectively.¹⁵ 

Dipole polarizabilities were calculated analytically in the cc-pVDZ basis set,⁶⁰ except for coupled-cluster with singles and doubles (CCSD) data where the finite-field (FF) method was employed. We confirmed the agreement between the FF and time-independent CCSD results utilizing the CCSD(3) model implemented in MOLPRO (for polyacetylenes up to four –|$\rm {C{=}C}$|$C=C$– units).^61,62 The geometries and IP-optimized range-separation parameters were taken from Ref. 58 (−|$\rm {C{=}C}$|$C=C$− distances reflect bond alternation). The polarizabilities were obtained with the following methods and program suites: MP2 and CCSD with MOLPRO,⁵⁹ PBE0 and PBE0AC with DALTON 2013,^63,64 and LRC-ωPBE with Qchem.⁶⁵ 

For the model |${\rm He}_3^+$|$He 3+$ system, the aug-cc-pVTZ basis set was chosen. The MP2 and CCSD calculations were performed in GAUSSIAN 09,⁶⁶ all other results were obtained with MOLPRO.

The A24 set of noncovalent complexes comprises small model systems selected in order to cover a wide spectrum of noncovalent interactions, i.e., hydrogen bonds, mixed electrostatics/dispersion, and dispersion-dominated interactions (including π–π stacking).⁵⁵ The benchmark energies used to compare our SAPT results were based on the CCSD with perturbative triples extrapolated to complete basis set limit, CCSD(T)/CBS, with an additional correction for perturbative quadruples, CCSDT(Q).

To go beyond the interactions involving exclusively second-row atoms, an additional set of benchmark systems was chosen. The DI6-04 basis set features six dipole interacting complexes containing third-row elements, sulfur and chlorine.⁵⁶ Benchmark geometries and energy values for DI6-04 were taken from Ref. 51, except for (H₂S)₂, whose geometry was reoptimized at the CCSD(T)/avtz level and the benchmark CCSD(T)/CBS value was obtained according to the expression of Halkier et al. with X = 5 (Eq. (4) of Ref. 67).

Recently Tsai et al. conducted a comparison of the LRC and AC-based methods.²⁷ Their results demonstrate the superior performance of LRC functionals over the AC approach for a wide range of properties. In this work, we address the shortcomings of the AC scheme that should directly affect the SAPT terms, complementing the findings of Ref. 27. The shortcomings in question offer an indirect measure of the persistence of the self-interaction error. We emphasize that in the case of asymptotically corrected potentials, a direct test of the dependence of the energy with respect to fractional electron numbers²⁵ is impossible. This is due to the non-uniqueness of the AC energy, as stray potentials do not correspond to any density functional.

A known failure of local and semilocal DFAs is the catastrophic overestimation of polarizabilities and hyperpolarizabilities for long-chained molecules.^30,68 Several studies proved that LRC functionals can rectify this problem, at least qualitatively.^{28,44,69–71} Figure 1 shows longitudinal static polarizabilities of polyacetylene chains up to |$\rm {C_{24}H_{26}}$|$C24H26$ predicted by different levels of theory (the remaining components of the polarizability remain coherent within all the methods as shown in Table S1 of supplementary material^72,73). The results show that PBE0 and PBE0AC fail identically, thus indicating the lack of the field-counteracting term. The same conclusion applies to the Casida-Salahub AC scheme, as observed by Sekino et al.²⁸ An overestimation in the linear response for long molecules should translate into inaccurate FDDSs and overestimated second-order induction and dispersion energy contributions in SAPT.

π-conjugated systems may pose a challenge for LRC-SAPT as well, as the IP-tuned ω values decrease rapidly with the chain length. This behavior was first observed by Körzdörfer et al. and attributed to the varying extent of the π-conjugated character of the system.⁵⁸ An alternative global-density-dependent (GDD) tuning procedure, termed ω_GDD, as proposed by Modrzejewski et al.⁷⁴ easily ameliorates this spurious behavior (Fig. 1). In this model, the optimal ω value is related to the average distance between the outer electron and an exchange hole. It is expected to lead to significantly improved response properties.

The lack of the field-counteracting term in the exchange potential should impact not only the static but also the dynamic response of the molecule. In particular, one expects incorrect behavior of dynamic polarizabilities at imaginary frequencies to yield systematically overestimated dispersion energy (see also Misquitta et al.).⁷⁵ In order to illustrate this effect we calculated isotropic C₆ coefficients for the interaction between two polyacetylene chains. The results presented in Table I follow the trends observed for longitudinal dipole polarizabilities: (i) the discrepancy between the DFA and CCSD reference values increases together with the chain length, (ii) PBE0AC constitutes no improvement over PBE0, and (iii) the optimally tuned CAM-B3LYP*⁷⁶ functional shows similar n-dependence as HF values. We conclude that for extended systems with fluxional electron distributions, SAPT based on PBE, PBE0, and PBE0AC will yield overestimated dispersion energies.

MSIE may adversely affect the distribution of the electron density. To demonstrate that AC does not cure the delocalization error of standard DFAs we applied it to the stretched |${\rm He}_3^+$|$He 3+$⁠, a model system studied by Körzdörfer and Brédas.⁴⁵ In this collinear system, the intermolecular distance R between the central and peripheral atoms is varied from 1.2 to 7 Å. Figure 2 shows the variations of the charge on the central He as a function of R obtained by different methods. As evidenced by the full configuration interaction (FCI),⁴⁵ with increasing R the positive charge should localize on the central helium atom. Both PBE0 and PBE0AC, as well as an untuned LRC-ωPBE(ω = 0.3), lead to uniform delocalization of electron density over three atoms (Fig. 2). Such erroneous charge density will result in inaccurate electrostatic interaction energy with any other subsystem. By contrast, the IP-tuned LRC-ωPBE correctly assigns whole positive charge to the central He atom ensuing correct electrostatic energy.

The distortion of electrostatics due to the delocalization error can be expected for many chemically interesting systems, e.g., in molecular electronics.^77–79 To better appreciate the consequences to SAPT, let us recall the case of the tetrathiafulvalene diquinone anion (Q-TTFQ)⁻ where Vydrov and Scuseria⁷⁷ demonstrated that LRC correctly localizes the negative charge on half of the molecule while both PBE and PBE0 erroneously delocalize the charge over the entire molecule. A SAPT(PBEAC/PBE0AC) calculation for this system as an interaction partner would no doubt be wrong.

Having discussed challenges for DFT-SAPT applications, we now move to the performance test of different approaches to the method. In Table II, we present mean absolute errors (MAEs) and root-mean-square deviations (RMSDs) for DFT-SAPT calculations of the A24 and DI6-04 data sets, as compared to benchmark energies. The entries from the left to the right gradually move from a KS to a GKS description of monomers.

For the A24 data set, there is a clear improvement going from PBE to PBEAC to PBE0AC to LRC functionals, i.e., from KS to GKS description of monomers. We observe that optimally tuned LRC functionals stay in close agreement and outperform asymptotically corrected PBE0AC, considered as the functional of choice for the DFT-SAPT method.^5,80

In the case of the dipole-bound systems in the DI6-04 data set, the situation is different, as the increased fraction of the HF exchange correlates with a considerable overbinding. There is a clear deterioration from PBEAC to PBE0AC and analogously from LRC-ωPBE to LRC-ωPBEh, with PBE0AC and LRC-ωPBE performing almost identically. Similar overbinding in DFT-SAPT was reported for dispersion dominated systems containing phosphorus, (P₂)₂ and (PCCP)₂, in Ref. 80.

An analysis of the SAPT energy components attributes this overbinding to the imbalanced reduction between |$E^{(1)}_{\rm exch}$|$E exch (1)$ and |$E^{(2)}_{\rm disp}$|$E disp (2)$ in GKS with respect to SAPT(PBEAC) (see Fig. S1 in the supplementary material⁷²).

In Figs. 3 and 4, we provide percentage errors of the DFT-SAPT total interaction energies for systems of A24 and DI6-04 data sets, respectively. The underbinding observed for the several dispersion-bound systems of the A24 data set, i.e., CH₄…CH₄, CH₄…C₂H₆, Ar…CH₄ and Ar…C₂H₄ results from the aug-cc-pVTZ basis set being too small to saturate the dispersion interaction. Extending the basis set to aug-cc-pVQZ reduces the errors (not shown). One should bear in mind that the basis set error introduced by the application of aug-cc-pVTZ influences the overall comparison of the DFT-SAPT approaches. Nevertheless, a similar performance of LRC-SAPT and AC-SAPT in larger basis sets should hold, as shown in the work of Cencek and Szalewicz.²¹ 

Let us more closely examine the behavior of the dispersion energy in different realizations of DFT-SAPT. Williams and Chabalowski were first to notice a correlation between induction and dispersion terms and the gaps between occupied and virtual orbitals in the uncoupled second-order SAPT.^6,8 In the KS scheme, these differences are underestimated, which manifests itself in HOMO-LUMO (HL) gaps being significantly narrower than fundamental gaps. (We note that there is a difference in physical interpretation of HL gaps obtained from Kohn-Sham and Generalized Kohn-Sham approaches, see Refs. 81 and 82). In the GKS scheme, realized by hybrid and LRC functionals, HL gaps widen as a fraction of Hartree-Fock, the self-interaction-free exact-exchange term, is introduced.⁸³ In hybrid functionals, application of the AC ensures proper asymptotic behavior of the v_xc potential and shifts the HOMO level to match the vertical ionization potentials. This, however, is insufficient to raise LUMO to match HOMO of the negative ion because of the derivative discontinuity (DD). Consequently, even in asymptotically corrected hybrid functionals HL gaps remain too narrow. On the other hand, absorbing the DD by the optimal tuning procedure in LRC functionals leads to almost perfect correspondence between the HL orbital gaps and the (accurate) fundamental gaps from Δ self-consistent field (SCF) computations.^58,74,84,85

In Fig. 5, we present correlation between the dispersion energy, |$E^{(2)}_{\rm disp}$|$E disp (2)$⁠, and HL gaps for homodimers from the A24 data set. As expected, a reduction of the dispersion interaction at the uncoupled level of theory follows a broadening of the HL gap from KS to AC KS and AC GKS to LRC GKS. A similar dependence is observed for the coupled treatment of dispersion (Fig. 5). It is expected that both KS and GKS formalisms would yield “exact” response properties, provided that the exact underlying (local or nonlocal) xc potentials and kernels in both theories were known and could be calculated. Results presented in Fig. 5 reveal that the LRC description of orbital energies corresponds the observed underestimation of the dispersion energy in LRC-SAPT with respect to AC-SAPT and KS-SAPT methods (see Sec. III E. and Table S1 in the supplementary material⁷²). Still, the LRC functionals were shown to provide a reliable description of response properties,^27,86 such as static^{28,44,69–71} and dynamic polarizabilities^20,46 which benefits the quality of dispersion energy at the coupled level of SAPT.^87,88

Nevertheless, systems characterized by large ω values pose a particular challenge for LRC-SAPT, as exemplified by dimers containing rare-gas atoms. Increasing the total fraction of HF exchange too much may distort the balance in error compensation between correlation and exchange functionals, resulting in the incorrect behavior of the energy contributions.^89,90 We show this by plotting the first-order, E⁽¹⁾ = E_elst + E_exch, and dispersion, |$E^{(2)}_{\rm disp,tot} = E^{(2)}_{\rm exch-disp}\break + E^{(2)}_{\rm disp}$|$E disp , tot (2)=E exch − disp (2)+E disp (2)$ energies, for the Ar…CH₄ and CH₄…CH₄ complexes (see Fig. 6). The first-order contributions remain in excellent agreement between SAPT(PBE0AC) and LRC-SAPT, indicating that the tuned ω values yield good electron densities. The problem lies in the dispersion energy, which in the case of LRC-SAPT acquires a more HF-like character, leading to the observed underbinding.

Previous studies of the LRC-SAPT method focused on He₂ and Ne₂ dimers as the simplest representatives of dispersion-bound systems.^20,48 Though excellent agreement with benchmark values was found for He₂, the results for the neon dimer were less satisfactory. Lao and Herbert interpret this as the problem of unsaturated dispersion interaction for the Ne₂ system. However, this could not be the case in the work of Cencek and Szalewicz who calculated the second-order energy contributions at the coupled level of theory in a large basis set. Our results provide the following explanation: due to the large ω value, the dispersion energy becomes significantly shifted towards the Hartree-Fock limit, as seen in the Ar…CH₄ system (Fig. 6). For He₂, this shift toward Hartree-Fock fortuitously improves the result with respect to the full CI benchmark, whereas the opposite is true for Ne₂ (see Fig. S2 in the supplementary material⁷²).

In this section, we examine the role of the short-range HF exchange in LRC-SAPT in hydrogen-bonded systems as exemplified by the water dimer. The data in Table III have been recomputed in the geometry of Ref. 91 to match that used by Cencek and Szalewicz.²⁰ We intend to demonstrate that the accurate performance of the LRC method extends to a much wider class of IP-tuned functionals with variable short-range HF exchange. To this end we present a family of LRC-ωPBEα functionals [Eq. (3)] that increase the total percentage of the HF exchange (%HF) so as to gradually morph into the SAPT based on HF monomers. We recall that α denotes the constant fraction of HF in the functional. To avoid basis set issues, the SAPT(CCSD)^92–94 reference values have been recomputed in the same basis set as the other SAPT data in the table.

The total fraction of the HF for a given geometry, i.e., originating from both long- and short-range exchange, was calculated in three steps. First, the full functional was applied saving the vectors. In the second and third step, we performed no-SCF calculations starting from the saved vectors, with the DFT and HF exchange switched off, respectively.

As shown in Table III, the LRC-ωPBEα functionals afford total interaction energies in very good agreement with the reference data for a wide range of α values from 0 to 0.4. This range corresponds to the total percentage of HF from 27% to 50%. Above 50% E_int begins to drop and quickly becomes HF-like. We have discussed similar deterioration in rare gas dimers above. The first-order exchange energy decreases gradually as the %HF rises. By contrast, the dispersion energy experiences a drop in binding between PBE0AC and LRC with ω = 0.498 despite nearly equal %HF (25% vs. 27%). Superficially, this drop can be rationalized by the increase in the HL gaps. More fundamentally, it can be attributed to the transition into monomers being described fully within GKS.

In the water dimer, as well as in the other studied dimers, the dispersion energy obtained in LRC functionals systematically underbinds compared to PBE0AC. The underbinding is partially offset by the less repulsive exchange-dispersion term. The LRC dispersion also departs from the SAPT(CCSD) benchmark by −8%. Cencek and Szalewicz²¹ report a similar departure of −7.6% in aug-cc-pVQZ basis set. Although such a dispersion is less accurate, one may argue that this is the price one pays for the removal of the delocalization error within the monomers.

LRC functionals have led to significant improvements over regular semilocal functionals in the description of a number of ground-state^47,85,95 and excited-state properties^96–99 that are relevant to the description of intermolecular forces. One should stress that these improvements were achieved, ultimately, as a result of the removal of MSIE. The SAPT formulation of the interaction energy relies on the quality of both types of properties. It is therefore pertinent to ask: what is the preferable DFT treatment of monomers to build SAPT theory on?

In this work, we examined the performance of SAPT built on GKS with optimally tuned LRC functionals. This treatment of monomers leads to the removal or the minimization of the effect of MSIE. This in turn corrects the erroneous delocalization of the ground-state density. As far as the excited-state properties are concerned, the optimized LRC functionals feature HL gaps being equal to fundamental gaps and a restored field-counteracting term in the v_xc potential. This in turn leads to an improved description of static and dynamic response properties. Our results show that the AC scheme does not overcome deficiencies of the underlying functional, such as spurious delocalization of density in |$\rm {He_3^+}$|$He 3+$ and overestimated response, both static and dynamic, in π-conjugated chains. Thus, it is justifiable to expect that SAPT based on LRC functionals will be more predictive and of wider applicability than that based on the present AC schemes.

A comparison of the two approaches confirms that tuned LRC functionals yield results comparable or better than the best AC functional based on PBE0. Nevertheless, a detailed study of the water dimer reveals that PBE0AC energy components compare better with the SAPT(CCSD) values. To conclude, while the AC-SAPT interaction energy components taken separately are closer to the benchmarks, the total LRC-SAPT interaction energies are better than the AC-SAPT ones. One may infer from this observation that while AC-SAPT may provide better energy components, the cancellation of errors seems more evenly balanced within LRC-SAPT. Clearly, a more thorough comparison with SAPT(CCSD) energies is necessary to explain the origin of the observed error cancellation in LRC-SAPT for its future applications.

Two main requirements of LRC-SAPT should be recognized. The first one is related to the strong dependence of SAPT terms on the value of range-separation parameter which makes optimization of ω with good precision a necessity.²⁰ For large molecular systems, combining LRC-SAPT with more effective functional tuning schemes, such as ω_GDD,⁷⁴ is recommended. The second prerequisite is related to the applicability of the optimally tuned LRC functionals themselves.^{36,77,100,101} Exceeding ∼50% of the HF exchange may distort the balance in error compensation between exchange and correlation within the functional, even in carefully tuned functionals. This happens in a limited number of cases, such as rare gases and other molecules with tightly bound electrons, where the optimal ω is exceptionally large, that is, the switch to HF exchange occurs at relatively short distances. This makes the interaction too HF-like.

M.H. was supported by the project “Towards Advanced Functional Materials and Novel Devices—Joint UW and WUT International PhD Programme” operated within the Foundation for Polish Science MPD Programme, implemented as a part of the Innovative Economy Operational Programme (EU European Regional Development Fund). Ł.R. performed his work under the auspices of a scholarship provided by the Alexander von Humbold Foundation. G.Ch. wishes to acknowledge the support by the Polish Ministry of Science and Higher Education, Grant No. N204 248440, and by the National Science Foundation (US), Grant No. CHE-1152474.