The nature of three-body interactions in DFT: Exchange and polarization effects

The decomposition of the interaction energy into many-body contributions provides insights into the effects beyond pairwise contacts. Although 3- and higher-body components are significantly weaker than 2-body ones, the number of the former increases rapidly with the cluster size.¹ They cannot be neglected when predicting the structure and properties of clusters of water,^2–8 clathrate hydrates,^9,10 molecular crystals,¹¹ and condensed-phase systems.^1,12 Finally, the quantitative description of nonadditive effects is indispensable when resolving the high-resolution spectra of molecular trimers and larger aggregates.^13–19 

For more than a decade, the applications of density functional theory (DFT) to noncovalent interactions have focused on the development of dispersion-correction schemes accounting for the dispersion effects missing from hybrid and semilocal density functional approximations (DFAs).²⁰ As a result, the dispersion-inclusive DFT methods achieved statistical errors within a few tenths of a kcal/mol with respect to ab initio benchmark databases for noncovalent dimers.^21–24 Nevertheless, the studies of larger aggregates drew attention to the importance of many-body errors in DFT and to their origin.^{4–6,9,25,26}

Recent work lists the following sources of DFA errors in binding energies:⁷ (a) errors in the intersystem exchange interaction, (b) the delocalization error, and (c) the monomer relaxation errors. The error in the exchange interaction stems from the inability of the semilocal DFAs to describe the electron exchange in the tails of overlapping densities, i.e., regions of significant change of reduced density gradients $s=|∇ρ|∕ ρ 4 ∕ 3$ upon formation of a noncovalent bond.^6,27,28 In the language of perturbation theory, this translates into the flawed treatment of the exchange terms arising from the Pauli exclusion principle. The delocalization error, also known as many-electron self-interaction error, refers to the fact that semilocal DFAs yield electron densities that are too delocalized.^23,29 This is pronounced in systems with strong induction and electrostatic interactions, although one should not forget that all interaction energy components are affected.³⁰ Finally, the monomer relaxation error, which contributes to the binding energies of the cluster, concerns the ability of the functional to account for the energetics of bond stretching due to the formation of the complex. The inclusion of pairwise dispersion corrections cancels out some of these errors.^7,31

Several studies have recognized the flaws in exchange functionals as a source of error in 3- and higher-body terms.^6,25,32 The most comprehensive test of DFAs encompasses a comparison of DFT and coupled-cluster total 3-body energies on the 3B-69 database of 23 molecular trimers in 69 configurations.³² A common suggestion in the aforementioned studies is that the errors arising from the incorrect many-body exchange may be more important than the missing many-body dispersion interactions. Clearly, this issue has to be addressed in order to improve the reliability of DFT predictions for large aggregates.

The conclusions in the existing literature (with only some exceptions—see Ref. 25) are based on inference from the supermolecular computations of many-body terms which bundle together effects of a different physical origin and geometry-dependence. Unlike in the wave function methods, such as Møller-Plesset perturbation theory, where the many-body interaction energy components residing in each order are known,¹⁸ such insights are missing in the case of DFAs.

The performance of DFT for 3-body interactions can be understood given a physically meaningful decomposition of the interaction energy. To this end, we define the decomposition of the total 3-body interaction energy into the exchange and polarization components. Our aim is to reveal which nonadditive contributions to 3-body interaction energy are captured by different DFAs and how accurate this description is in DFT. This will be done by comparing the DFT exchange and polarization components with rigorously defined many-body terms from the perturbation expansion of the interaction energy.

The assessment of the DFAs’ ability to account for the exchange nonadditivity involves the following strategy. To extract the exchange part from a supermolecular DFT 3-body term, we evaluate the interaction energy between the Löwdin-orthogonalized monomers, which is the DFT equivalent of the Heitler-London exchange nonadditivity.³³ Next, we substitute the Hartree-Fock formula for the intermolecular part of the exchange potential to investigate how switching from the DFT to HF description changes the exchange nonadditivity. This is carried out with the aid of the dispersion-free Pauli blockade (PBdf),³⁴ the method which we generalized here to the many-body case. Finally, nonadditive exchange energies extracted from both supermolecular DFT and PBdf are compared with the Heitler-London term, for which we provide benchmark data. This strategy serves as the first step toward revealing what 3-body contributions are captured in the supermolecular DFT interaction energies and how trustworthy their representation in specific DFT approximations is.

In the numerical section, we focus our attention on a range of challenging systems: from rare-gas trimers, to Ar₂–HX (X = F, Cl), and to selected hydrogen and dispersion-bonded molecular trimers. Based on the presented energy decomposition scheme and 3-body symmetry-adapted perturbation theory (SAPT) analysis, we draw conclusions on the performance of DFAs.

The total 3-body nonadditive interaction energy is defined in terms of trimer, dimer, and monomer total energies as

Since our focus is solely on nonadditive interactions in trimers, we omit the [3, 3] label usually applied in the literature. To analyze DFAs, we decompose E_int into the nonadditive exchange and deformation contributions,

The nonadditive exchange captures the energetic effect of applying the Pauli exclusion in the Kohn-Sham wave function upon formation of the trimer out of non-interacting monomers. The computation of $E nadd–ex$ for an arbitrary DFA consists of non-iterative trimer and dimers’ energy evaluation from the Löwdin-orthogonalized³⁵ orbitals of the isolated monomers.³³ E_int is computed directly from Eq. (1). The deformation energy corresponds to the energy lowering upon mutual self-consistent polarization of the monomers restrained by the Pauli exclusion principle and is computed as

An analogous energy decomposition for dimers was first proposed by Cybulski and Seversen.^33,36 E_nadd–ex is the DFT version of the nonadditive Heitler-London interaction energy, a well-established concept in the theory of intermolecular interactions,^33,36,37

where $Φ 0 = Φ 0 A Φ 0 B Φ 0 C$⁠, $Φ 0 μ$ is the wave function of the monomer $μ$⁠, $A$ is the anti-symmetrizer that exchanges electrons between monomers, and $E i n t H L$(XY) is the Heitler-London interaction energy of the XY dimer.

The energy partitioning introduced in Eq. (2) is physically appealing due to the fact that benchmark values for the E_nadd–ex term can be calculated directly from SAPT. The reference for E_nadd–ex, denoted as $E i n t H L ( K S )$⁠, is computed as the Heitler-London interaction energy with Kohn-Sham determinants.

In order to use the existing 2-body and 3-body codes for SAPT based on the Kohn-Sham description of monomers (DFT-SAPT), we compute this term as

where the first-order DFT-SAPT nonadditive exchange energy is

and $Δ M$ denotes the Murrell delta term,³⁸ also referred to as the zeroth-order exchange energy. (See the supplementary material for explicit formulas for the computation of $Δ M$⁠.)

The reference $E i n t H L ( K S )$ value depends only weakly on the underlying DFA, provided that an asymptotic correction³⁹ or a tuned range-separated functional^40,41 is applied (see Table S1 in the supplementary material for the sensitivity test of $E i n t H L ( K S )$⁠). The correct asymptotic behavior of the exchange-correlation potential yields reliable monomer densities, a well-known prerequisite for DFT-SAPT calculations.⁴² 

Based on the DFT-SAPT performance for the two-body first-order exchange energy,⁴³ we expect that the 3-body $E e x c h , S A P T ( 1 )$ term is most accurately reproduced when monomers are described with asymptotically corrected hybrid functionals. Consequently, we will use $E i n t H L ( K S )$ based on the PBE0^44,45 functional with the gradient-regulated connection (GRAC)⁴⁶ asymptotic correction as a reference for the nonadditive exchange.

Depending on the quality of the exchange-correlation functional, E_nadd–ex should approach the reference $E i n t H L ( K S )$ value. While both E_nadd–ex and $E i n t H L ( K S )$ are computed with approximate DFT methods, the latter replaces the exchange-correlation functional in the intermonomer region with the Heitler-London exchange. The description of the intramonomer correlation enters in both quantities through the use of Kohn-Sham orbitals.⁴⁷ 

It should be stressed that the major part of the exact nonadditive exchange energy is captured already at the Hartree-Fock level of theory.⁴⁷ In particular, the E_nadd–ex contribution which is included in the supermolecular 3-body Hartree-Fock interaction energy is equal to the $E i n t H L$ formula calculated with Hartree-Fock determinants. Obviously, Hartree-Fock misses the intramonomer correlation contribution to E_nadd–ex, which in the Møller-Plesset perturbation theory enters at the MP2 and higher levels.¹⁸ 

Several comments regarding the E_def component should be made. At the Hartree-Fock level of theory, the deformation contribution collects second- and higher-order induction energy components together with their exchange counterparts while the lowest-order intramonomer contributions to those terms enter at the MP2 level.¹⁸ We expect that the nonadditive many-body induction terms, including the intramonomer correlation contributions, should be recovered quantitatively by the existing DFAs and their description should systematically improve upon the minimization of the delocalization error.^7,32 On the other hand, the third-order induction-dispersion energy (along with its exchange counterpart) and second-order exchange-dispersion energy—both included at the MP2 level—are most probably not accounted for due to the semilocal nature of DFAs. We will investigate this problem in more detail in Sec. III B.

A well-known component of the many-body nonadditive interaction energy is the third-order dispersion energy, $E d i s p ( 3 )$⁠, which in the Møller-Plesset theory first appears at the MP3 level. Much attention has been devoted to the proper inclusion of this term in DFT.^20,48 We address the importance of $E d i s p ( 3 )$ relative to the remaining second- and third-order terms in Secs. III B and III C.

To expose the DFA errors in the description of intermonomer regions, we introduce a scheme in which the Kohn-Sham exchange in the intermonomer region is replaced by the exact Hartree-Fock expression (Fig. 1). This approach eliminates the dispersion component of the interaction energy and is referred to as dispersion-free Pauli blockade.³⁴ In this section, we briefly introduce the idea behind many-body PBdf and comment on the contents of the nonadditive 3-body PBdf interaction energy.

In the PBdf formalism, the Kohn-Sham equation for the orbitals of monomer $μ$ in a system composed of L weakly interacting monomers is

where $f μ ( r )$ is the Kohn-Sham operator of the monomer $μ$⁠, $v ν ( r )$ describes the nuclei-electron attraction, $j ν ( r )$ is the Coulomb electron-electron repulsion, and $k ν ( r )$ is an exact HF-like exchange operator. Note that $f μ ( r )$ corresponds to the unaltered Kohn-Sham operator for a given DFT approximation.

Equation (7) is solved in an iterative fashion, keeping the monomer orbitals orthogonal at all times. This may be achieved with the use of either the penalty function⁴⁹ or the exponential ansatz of orbital rotation.⁵⁰ The latter scheme is more efficient and numerically stable and therefore has been used in the current implementation of the PBdf method. The iterative process starts from the Löwdin-orthogonalized³⁵ orbitals of the isolated monomers.

The converged orbitals are used to compute the dispersion-free energy of the complex,

where $E μ [ ρ ̃ μ ]$ is the monomer energy functional, $ρ ̃ μ$ denotes the density of the monomer $μ$ obtained from the converged orthogonalized orbitals, the electrostatic interaction term reads

with $W μ ν$ denoting the nuclear-nuclear repulsion energy, and E_exch is the HF exchange interaction between monomers.

The PBdf energy formula, Eq. (8), is used to compute the dispersion-free total nonadditive 3-body interaction energy, $E i n t d f r e e$⁠, according to Eq. (1), which requires converging the orbitals in Pauli blockade equations for the trimer and for the dimers. The energy partitioning of Eq. (2) yields dispersion-free counterparts of nonadditive exchange and deformation energy contributions,

where the $E n a d d − e x d f r e e$ term is obtained in the same manner as E_nadd–ex for the full DFA functional, i.e., employing orbitals of the isolated monomers which are Löwdin-orthogonalized but otherwise kept unperturbed. As will be shown, because of the enforcement of the exact exchange, $E n a d d − e x d f r e e$ varies much less than E_nadd–ex among different DFAs.

While the prescription for the PBdf energy is based on a simple model, it provides a useful diagnostic for the insight into the source of qualitative errors in the DFT exchange. A detailed derivation of the Pauli blockade and dispersion-free Pauli blockade schemes may be found in the supplementary material and in Refs. 34 and 49.

For He we chose the doubly augmented quintuple-zeta basis of Ref. ⁵¹. The argon trimer and Ar₂–HX (X = F, Cl) systems were studied with the aug-cc-pVQZ basis of Dunning.⁵² The molecular trimers from the 3B-69 basis³² and Ref. 7 were studied with the aug-cc-pVTZ basis. All 3-body energies were counterpoise-corrected.⁵³ The calculations with MVS⁵⁴ and SCAN⁵⁵ functionals were done on a large grid of 300 radial and 1202 angular points of the Lebedev grid. The IP-optimized values of the range-separation parameter $ω$ for the $ω$PBE⁵⁶ functional (also known as LRC-$ω$PBE), $ω$ = 0.55 for Ar and $ω$ = 1.0247 for He, were taken from Ref. 57. For molecular trimers studied in Sec. III C, $ω$PBE was used with $ω=0.5$⁠, as recommended in Ref. 7. The LC-PBETPSS functional⁵⁸ represents the class of meta-generalized-gradient-approximations (meta-GGA) range-separated hybrids. All LC-PBETPSS functional calculations were performed with the default $ω=0.35$ value. The DFT-SAPT calculations employed the PBE0 functional^44,45 with the GRAC⁴⁶ asymptotic correction. The first-order nonadditive exchange and second-order nonadditive exchange-dispersion SAPT contributions were obtained without overlap expansion, i.e., they are correct to all orders of the intermonomer overlap.^59–63 Additionally, for molecular trimers analyzed in Sec. III C, we calculated second-order exchange-dispersion energies at the uncoupled Hartree-Fock level [$E e x c h − d i s p ( 2 ; 0 ) ( S 2 + S 3 )$] with the SAPT2012 code.^47,64–66 All 3-body PBdf and DFT-SAPT calculations were done in the developer’s version of Molpro.⁶⁷ Calculations with the $ω$B97XD3⁶⁸ functional were done in the Orca program.⁶⁹ 

The failure of supermolecular DFT in predicting 3-body interactions in rare gas trimers was studied in the past.^6,70 The wrong description of regions with high reduced density gradients translates into the flawed nonadditive exchange interaction energies in the cluster.^6,7,25 This assessment of the DFT performance is evidenced in Table I, where the total 3-body interaction energies seriously deviate from both Hartree-Fock and coupled-cluster values. Because the tested functionals do not include any dispersion correction, it is reasonable to expect the DFT energies to be close to the Hartree-Fock ones.

For the comparison of DFAs, we chose two functionals of the GGA rung which feature a different asymptotic behavior of the exchange enhancement factor F_x(s): BLYP⁷¹ and PBE,^44,72 as well as their hybrid counterparts: B3LYP⁷³ and PBE0. In addition, the optimally tuned $ω$PBE functional represents the class of range-separated functionals.

The functionals characterized by the rapid increase of F_x(s) (BLYP and its hybrids) predict much too attractive nonadditive exchange energy, while the slow increase of F_x(s) in PBE-based functionals results in much too repulsive E_nadd–ex (Table I). It is worthwhile to note that a similar behavior of the 3-body exchange terms in the B88- and PBE-based functionals was recognized before in the analysis of water hexamers.²⁵ 

Two factors determine the accuracy of the DFT nonadditive exchange energy: the quality of density in the region of large reduced density gradients and the approximations in the exchange-correlation energy formula.^6,7 The progression of results from semilocal PBE, to a global hybrid PBE0, and to a range-separated $ω$PBE functional shows that including a fraction of Hartree-Fock exchange is mandatory and the long-range corrected form of the functional brings the biggest improvement (Table I). For the PBE-based functionals, the long-range corrected form of the functional is imperative to account even for the correct sign of E_nadd–ex. Next, the effect of the HOMO, which determines the tail of the density, is also important. For example, for He₃ the $ω$PBE functional with the default $ω$ parameter corresponds to E_nadd–ex = −3.9 $μ$hartree, whereas the variant of this functional, which enforces Koopmans’ theorem, corresponds to E_nadd–ex = −2.3 $μ$hartree. Note that the sign of E_nadd–ex is correct in B88-based functionals, but the inclusion of the exact exchange results in a further substantial improvement.

The Heitler-London exchange energy can be reproduced accurately provided that the interaction between the monomers is described fully at the Hartree-Fock level. This can be achieved by means of PBdf, where the monomers interact with the exact Hartree-Fock potentials, but the intramonomer contributions to the interaction energy are described at the DFT level. Applying the $ω$PBE functional in the PBdf framework reduces the error in the Heitler-London energy by an order of magnitude with respect to $ω$PBE for He₃ and Ar₃ (Table II). Somewhat surprisingly, the remaining discrepancy between PBdf and $ω$PBE suggests that the full Hartree-Fock exchange at long range in the range-separated hybrid is not sufficient to capture the effect of the nonadditive exchange. The reference values which PBdf should approach, provided that the underlying functional accurately describes intramonomer correlation, are the 3-body interaction energies at the MP2 level minus the second-order exchange-dispersion term. We reiterate that while MP2 does not include the 3-body dispersion energy, it contains the uncoupled second-order exchange-dispersion.^18,74 The latter was shown to be important in rare gases.^47,75–77 (The induction correlation effects are also present but for rare gas trimers are negligible.)

The paradigm trimers for studies of nonadditive forces, Ar₂–HF and Ar₂–HCl, were the first systems for which it was possible to extract the 3-body potential based on the microwave^78,79 and far-infrared^80,81 spectroscopic data, as done in the studies of Hutson et al.¹³ and Cooper and Hutson.¹⁴ The nature of nonadditive interactions in Ar₂–HX was later thoroughly studied in a series of papers combining supermolecular MPn calculations and direct calculations of nonadditive energy components based on perturbation theory.^77,82–85 3-body interactions manifest themselves most clearly in the regions of in- and out-of plane rotations of HF and HCl.

As shown in Refs. 82 and 83, the dominant contribution to the anisotropy in Ar₂–HX comes from the exchange nonadditivity, yet the induction and dispersion components cannot be neglected. The leading contribution to the third-order induction comes from the moments induced on two argon atoms by the field of HX. The exchange and induction nonadditivities are more pronounced in Ar₂–HF as the distance between HF and the center of mass of Ar₂ is considerably smaller than in Ar₂–HCl, and the HF dipole is larger than that of HCl.

In the present work, we focus only on the in-plane rotations of HX, where $Θ$ is varied from $0 °$ (the global minimum) to $18 0 °$ (Fig. 2). $Θ=$ $0 °$ corresponds to H pointing toward the center of the triangle. The bond lengths for Ar₂–HCl are^82,83 r_HCl = 1.275 Å, r_Ar–Ar = 3.861 Å, and $R=3.4736Å$⁠. The values for Ar₂–HF are r_HF = 0.917 Å, r_Ar–Ar = 3.826 Å, and $R=2.9798Å$⁠.

The performance of DFT methods for the total 3-body interaction energies of Ar₂–HX is qualitatively the same as in the case of noble gas trimers. The BLYP and B3LYP functionals are in better agreement with both Hartree-Fock and MP2 results than the PBE-based functionals (Fig. 3, the results for PBE and PBE0 are given in Fig. S1 of the supplementary material). The behavior of the PBE-based approximations improves upon introducing range separation of the exchange functional. Finally, PBdf and Hartree-Fock, which do not account for the second-order exchange-dispersion and induction-dispersion effects, closely follow each other for Ar₂–HX.

The semilocal (PBE) and global hybrid (PBE0) functionals based on the PBE model heavily overestimate the nonadditive exchange, E_nadd–ex, with respect to the reference $E i n t H L ( K S )$ results based on DFT-SAPT (Fig. 4). As in noble gas trimers, the errors are ameliorated in $ω$PBE. The PBE-based range-separated functional on the meta-GGA rung, LC-PBETPSS, performs similarly to $ω$PBE, which confirms that our conclusion on the role of the range-separated exchange is general. In addition, the B88-based functionals stay in better agreement with the reference Heitler-London energy than the PBE-based approximations.

It is tempting to check if going beyond PBE- and B88-type enhancement factors leads to any improvement in the description of the nonadditive exchange. To this end, we examined two recently developed meta-GGA functionals: MVS⁵⁴ and SCAN.⁵⁵ The enhancement factors built into these formulas have a radically different asymptotic behavior from those in PBE and B88.⁵⁴ The MVS/SCAN enhancement factor satisfies

where $α$ is a dimensionless variable that depends on the kinetic energy density.⁵⁴ In contrast, the BLYP enhancement factor is approximately linear for large s, whereas for PBE it approaches a constant.

The new F_x(s), however, does not lead to an improvement for 3-body energies in Ar₂–HX. For $Θ∈ ( 0 ° , 6 0 ° )$⁠, both MVS and SCAN predict wrong angular dependence of the total interaction energies (Fig. 5). The error can be attributed to the behavior of the nonadditive exchange part (for Ar₂–HCl results, see the supplementary material).

The spread of deformation energies computed with different DFT methods is much smaller than in the case of the exchange nonadditivity (Fig. 4). Long-range corrections no longer move the results towards Hartree-Fock ones. The curves corresponding to $ω$PBE and LC-PBETPSS deviate from both Hartree-Fock and global hybrids in the range of $Θ= 0 °$ to $6 0 °$⁠.

To understand whether the difference between range-separated DFT and Hartree-Fock is related to any physical effect, we compare

with two quantities derived from the MP2 interaction energy,

and

$Δ E ( 2 )$ is the correlation contribution to MP2 interaction energy, which contains intramonomer correlation contributions to nonadditive exchange and induction terms, second-order exchange-dispersion, and third-order induction-dispersion and exchange-induction-dispersion terms.^18,74 As shown before,^77,85 the induction-dispersion terms are partly canceled by their exchange counterparts.

Importantly, the correlation contribution to the nonadditive exchange is small in Ar₂–HX. This can be inferred from Fig. 4 where the difference between the Heitler-London exchange energy at the DFT-SAPT level of theory $E i n t H L ( K S )$ and Hartree-Fock nonadditive exchange is on the order of 1 $μ$hartree (see also Fig. S3 in the supplementary material). Among the physical contributions to $E d e f c o r r , D F T$⁠, we expect that the DFAs reliably account for the intramonomer correlation contributions to the assorted induction terms.

The range-separated functionals qualitatively agree with $Δ E ̃ ( 2 )$⁠, while semilocal functionals and global hybrids are close to $Δ E ( 2 )$ (Fig. 6). This indicates that the former excludes the second-order exchange-dispersion term, while the latter includes it. The $E exch–disp ( 2 )$ contribution, which is the exchange counterpart of (the additive) second-order dispersion, may appear in semilocal functionals only as an artifact, similar to the fictitious exchange binding of noble gas dimers. Note that the addition of the Hartree-Fock exchange in global hybrids partially removes the spurious angular dependence observed for the semilocal exchange at $Θ<6 0 °$⁠.

It should be stressed that although the second-order exchange-dispersion is a sizeable effect^77,85 (see Table S5 of the supplementary material), it should not contribute to the damping of the dispersion nonadditivity, the third-order Axilrod-Teller-Muto (ATM) term. This role should be reserved for the exchange counterpart of the third-order ATM dispersion (e.g., Refs. 47, 86, and 87). If $E e x c h − d i s p ( 2 )$ proves necessary in the context of DFT, it should be employed in a controllable manner. Therefore, range-separated functionals provide an appropriate starting point for the development of new dispersion-corrected approaches.

We apply our scheme for energy decomposition on a subset of the 3B-69 data set of Řezáč et al.³² (water, formaldehyde, methanol-ethyne, and acetonitrile trimers), three dispersion-bonded (CO, CH₄, N₂) and two H-bonded (HF, NH₃) trimers studied in Ref. 7. Only a subset of 3B-69 data set is included due to our computational limits. Our set includes seventeen configurations in total, see the supplementary material for details.

A different description of the exchange nonadditivity is the major source of variance among DFT approximations (see Tables III and V, Tables S7 and S8 in the supplementary material). The deformation energy, on the other hand, is much more consistent even between functionals based on different models of the enhancement factor, e.g., PBE and BLYP yield almost identical E_def for seventeen tested trimers.

Remarkably, in contrast to noble gas trimers and Ar₂–HX discussed above, in H-bonded systems, the exchange nonadditivity is well reproduced by range-separated hybrids. For the water trimer in cyclic configuration, $ω$PBE gives E_nadd–ex = −0.267 kcal/mol, which is in reasonable agreement with the reference $E i n t H L ( K S ) =−0.302$ kcal/mol (Table III). Both semilocal and global hybrid approximations predict exchange nonadditivities less accurately than range-separated hybrids. The MVS functional yields better E_nadd–ex than any other tested semilocal methods. Moreover, out of seven systems of predominantly H-bonded character, in four cases (HF, $NH 3$⁠, methanol-ethyne 3a, and 3c), MVS is more accurate than the global hybrids and approaches the accuracy of $ω$PBE (Tables S7 and S8 in the supplementary material).

The remaining error of $ω$PBE stems mainly from the lack of the exchange-dispersion and dispersion contributions as well as inaccurate deformation due to the residual delocalization error. For example, for the water trimer, adding the coupled Kohn-Sham $E e x c h − d i s p ( 2 ) =0.196$ kcal/mol moves $ω$PBE close to the MP2 result (Table III). Note that exchange-dispersion and other correlated contributions at the MP2 level largely cancel each other in this system, as indicated by the excellent agreement between Hartree-Fock and MP2.

To test how the addition of the exchange-dispersion term affects DFT performance, we added this contribution (at the uncoupled SAPT Hartree-Fock level) to a set of pure, global hybrid, and range-separated functionals (Table S9 in the supplementary material). As expected from our discussion of Ar₂–HX, range-separated functionals benefit the most from the inclusion of $E e x c h − d i s p ( 2 )$⁠; the average errors of $ω$PBE and $ω$B97XD3 on the set of seventeen trimers are reduced by one-third. In contrast, for PBE, PBE0, and BHLYP, the average errors become larger. This finding shows that the 3-body DFT dispersion correction should depend on the type of the base functional.

As seen in Table III, due to the delocalization error, pure GGA functionals yield inaccurate deformation in H-bonded clusters. The error diminishes upon admixture of a large portion of exact exchange (recommended 50%^7,32) or application of the long-range Hartree-Fock exchange in range-separated hybrids. Accordingly, $ω$PBE and BHLYP yield similar deformation energies for water trimers: E_def = −2.427 kcal/mol and −2.392 kcal/mol, respectively.

The good performance of PBE and PBE0 for H-bonded clusters relies on the cancellation of errors between the nonadditive exchange and deformation components (Table III). This is in accordance with the observation of Ref. 7. One should note, however, that the error cancellation in supermolecular calculations with PBE-based semilocal and global hybrid approximations does not work equally well for every configuration: both PBE and PBE0 typically overestimate the repulsive 3-body energy contributions.³² On the whole test set, dispersion-corrected PBE stands out with its large average error of MAE = 0.093 kcal/mol, which is only slightly improved by the introduction of the long-range Hartree-Fock exchange in $ω$PBE (Table IV).

The B88-based semilocal and global hybrid approximations improve the nonadditive exchange with respect to the PBE-based methods. Still, without the long-range Hartree-Fock exchange it is impossible to realistically describe E_nadd–ex. BHLYP supplied with a dispersion correction gives the most accurate total 3-body interaction energies on the whole set of seventeen trimers (Table IV). It was also identified as the best-performing functional in the paper of Řezáč et al.³² with the root-mean-square error of 0.045 kcal/mol compared to 0.059 kcal/mol for MP2. However, our energy decomposition has shown that BHLYP relies on error cancellation between the inaccurate exchange nonadditivity and other contributions. For example, if we substituted $E i n t H L ( K S )$ for the BHLYP exchange nonadditivity of the water trimer (Table III), the total 3-body interaction energy would be the same as that of $ω$PBE.

The 3-body effects in dispersion-bound trimers (⁠$N 2$⁠, CO, and $CH 4$⁠) are beyond the reach of DFAs even at the qualitative level. None of the DFT functionals correctly accounts for 3-body nonadditive exchange, similarly to noble gas trimers (Table V). In contrast, the dispersion-free PBdf approach gives a good description of exchange interactions for both dispersion- and H-bonded clusters. A plausible explanation for the distinction between the DFT treatment of hydrogen- and dispersion-bonded systems is that in the case of the latter, the energetic contributions to the interaction energy come from regions of relatively small densities.⁸⁸ 

This work has examined the components of nonadditive interactions obtained in the supermolecular DFT calculations. We proposed a physically meaningful decomposition of the total 3-body nonadditive DFT interaction energy into the exchange interaction term, which captures the effect of the Pauli exclusion in the wave function of the trimer, and the deformation energy, which originates from the mutual polarization of the monomers. This scheme proved to be a useful tool for understanding the performance and future improvements of DFT models for noncovalent many-body systems.

The major source of variance among DFT methods is the nonadditive exchange term. The behavior of this contribution can be discussed for two limiting cases: dispersion and hydrogen-bonded clusters. For the former, no DFT approximation reproduces nonadditive exchange even at a qualitative level, as shown, e.g., for rare gas trimers, Ar₂–HX, and methane trimer. For H-bonded systems, range-separated functionals are the only type of DFT approximations capable of describing the nonadditive exchange (to within 20% of the reference value). For the paradigm Ar₂–HX (X = F, Cl) systems, range-separated hybrids yield correct anisotropy of E_nadd–ex, but the interaction curves are significantly shifted from the reference values. This pertains to both GGA and meta-GGA range-separated hybrids, as nonempirical $ω$PBE and LC-PBETPSS, respectively, are tested. The MVS and SCAN functionals employing the latest ideas in the enhancement factor design fail in these trimers.

The DFT errors in the deformation energy are smaller than those of the nonadditive exchange, and the differences between functionals are less pronounced. The correlation effects on induction terms could be indirectly probed in specific instances such as Ar₂–HX. We showed that range-separated functionals can recover these effects semiquantitatively. Other DFA types do not permit a more detailed interpretation because their deformation effects are obscured by artifacts of the exchange functionals.

The second-order exchange-dispersion nonadditivity was shown to be very important in the wide range of systems.⁴⁷ This nonadditivity represents the exchange counterpart of the additive second-order dispersion term and is completely ignored in the existing ATM force-field-like dispersion corrections due to a completely different functional dependence. In the majority of studied situations, the second-order exchange-dispersion either exceeds or nearly equals the third-order ATM nonadditivity. For example, in the cyclic water trimer (1c), the coupled Kohn-Sham $E e x c h − d i s p ( 2 )$ equals 0.196 kcal/mol, whereas an approximate damped ATM term reported in Ref. 32 amounts to 0.028 kcal/mol. Our results for Ar₂–HX suggest that range-separated hybrids do not account for $E e x c h − d i s p ( 2 )$⁠, while semilocal and global hybrid functionals may imitate this contribution in the short range. Consequently, on a larger set of molecules, supplementing the existing dispersion corrections with the $E e x c h − d i s p ( 2 )$ term considerably improves the error statistics for range-separated functionals, but not for pure and global hybrid approximations.

The dispersion-free Pauli blockade scheme, in which the mutual polarization of the monomers occurs via HF operators, fully corrects the description of the nonadditive exchange energy for every functional and every system, in particular, for both hydrogen- and dispersion-bonded trimers. By comparison with DFT-SAPT, we have shown that PBdf can be used as a reference for the assessment of nonadditive exchange. It should be noted that in PBdf, the nonadditive exchange-dispersion term is not accounted for.

To conclude, in designing DFT approximations for the description of nonadditive 3-body interactions in weakly bound systems, the major challenges are related to the reliable modeling of the nonadditivities of the first-order as well as the exchange-dispersion effects. The former can be ameliorated in the particular case of hydrogen-bonded systems by the use of exact exchange at long range. The latter presents a challenge for all types of systems and warrants going beyond the expressions currently used for 3-body dispersion corrections.

See supplementary material for the derivation of the $Δ M$ term and PB(df) methods, sensitivity test of the $E i n t H L ( K S )$ term, and detailed numerical data for systems presented in Sec. III.

M.H. and M.M. were supported by the Polish Ministry of Science and Higher Education (Grant Nos. 2014/15/N/ST4/02179 and 2014/15/N/ST4/02170). This work was partly supported by the National Science Foundation (Grant No. CHE-1152474).