Finding chemical concepts in the Hilbert space: Coupled cluster analyses of noncovalent interactions

Noncovalent interactions (NCIs) play a major role in essentially all fields of chemical research. Energy decomposition analysis (EDA) schemes provide in‐depth insights into their nature by decomposing interaction energies into additive contributions, such as electrostatics, polarization, and London dispersion. Although modern local variants of the “gold standard” coupled‐cluster singles and doubles method plus perturbative triples (CCSD(T)) have made it possible to accurately quantify NCIs for relatively large systems, extracting chemically meaningful energy terms from such high level electronic structure calculations has been a long lasting challenge in computational chemistry. This review describes basic principles, interpretative aspects and applications of recently developed coupled cluster‐based EDAs for the analysis of NCIs. The focus is on computationally efficient methods for systems with a few hundred atoms, for example, the recently introduced local energy decomposition analysis. In order to draw connections between different interpretative frameworks, these schemes are compared with other popular approaches for the quantification and analysis of NCIs, such as Symmetry Adapted Perturbation Theory and supermolecular EDAs based on mean‐field as well as correlated approaches. Strengths and limitations of the various techniques are discussed.

Nowadays, a plethora of different quantum mechanical methods can be used to quantify NCIs at various degrees of accuracy. Although this is as a remarkable achievement of modern quantum chemistry, mere quantification does not improve per se our understanding of the abovementioned chemical processes. Typically, this requires the identification of families of NCIs with common properties (e.g., halogen bonds, hydrogen bonds, agostic interactions, etc.), the formulation of trends of interaction energies, and the development of simple models for their description. This kind of rationalizing activity is, as stated by Malrieu, 1 among the "most noble" tasks of quantum chemistry.
In order to facilitate the interpretation of the results obtained from quantum mechanical methods, a large variety of interpretative tools have been developed. In particular, energy decomposition analysis (EDA) schemes can be used to decompose NCI energies computed using a supermolecular approach into chemically meaningful additive components, such as electrostatics and London dispersion, thus providing a quantitative framework in which to describe the nature and properties of intermolecular interactions. As the energy contributions commonly found in EDA schemes are quantummechanically ill-defined quantities, different schemes typically provide different numerical values for the individual components of the interaction.
This review examines basic principles, interpretative aspects, applications and limitations of accurate correlated wavefunction-based methods for the quantification and analysis of NCI energies. Emphasis is placed on techniques applicable to systems with hundreds of atoms. A prominent example of this kind is the recently developed local energy decomposition (LED) scheme, [2][3][4] which relies on the Domain-based Local Pair Natural Orbitals (DLPNO) variant of the "gold standard" coupled cluster method with singles and doubles and perturbative included triples excitations, that is, DLPNO-CCSD(T). 5,6 This technique allows for the quantification and analysis of NCIs for virtually any system within the range of applicability of the CCSD(T) method, irrespective of the nature of the monomers, and of the strength of the interaction.
To put the discussion into a broader perspective, the fundamental theoretical aspects, strengths and limitations of the most popular perturbative and supermolecular schemes for the quantification and analysis of NCIs are also reviewed. Each scheme examined here is evaluated based on three criteria: 1. Accuracy. The observables of interest, for example, association or dissociation energies, should be computed with sufficient accuracy to describe the chemical problem of interest. There is not much to gain from the rationalization of a bad chemical trend. Ideally, all related quantities (one-and two-electron densities, geometries, etc.) should also be accurately described by the chosen level of theory. 2. Efficiency. The scheme should be efficient, so that it can be applied to small as well as to large molecular systems. In the context of the present review, a method is considered efficient if it can be routinely applied to systems of 100-200 atoms. 3. Range of applicability. Ideally, the scheme should remain valid for the broadest possible range of interaction types, irrespective of the nature of the interacting molecules. This is necessary in order to develop a universal understanding of NCIs.

| PERTURBATION THEORY OF NONCOVALENT INTERACTIONS
Typically, NCIs between atoms and molecules are weaker compared to covalent bonds. Hence, perturbative approaches appear to be a natural choice for their description. Besides their accuracy for the quantification of weak NCI energies, these schemes have the desirable feature of naturally providing an intuitive definition of their constituting electrostatic, induction and dispersion energy contributions, thus setting the framework for their description. These schemes have already been the subject of several reviews. [7][8][9] Hence, here we will only examine the basic principles of the various approaches and discuss their strengths and limitations.

| Basic principles
Consider two noninteracting monomers, X and Y. Their ground state wavefunctions are denoted as Ψ X 0 and Ψ Y 0 and the corresponding energies as E X 0 and E Y 0 : The exact wavefunction of the super system is given by the antisymmetrized product A Ψ X 0 Ψ Y 0 , whereÂ is the antisymmetry operator. The associated energy would thus be E XY 0 = E X 0 + E Y 0 . If the monomers are now brought into a weak interaction regime, we can use perturbation theory to estimate the variation in the total energy due to their interaction, that is, their interaction energy. The underlying idea here is to partition the Hamiltonian of the system as the sum of that of the monomers (Ĥ X andĤ Y ) plus the intermolecular interaction operatorV: A fundamental problem with this approach is that the eigenfunctions ofĤ 0 are simply the product of the eigenfunctions of H X andĤ Y , for example, Ψ X 0 Ψ Y 0 for the ground state. Hence, they do not satisfy the antisymmetry requirement. A perturbation series based on such zeroth-order solutions would not necessarily converge to the proper antisymmetric wavefunction of the molecular adduct. [10][11][12] Nevertheless, if Ψ X 0 and Ψ Y 0 do not overlap significantly, the interaction energy of X and Y can be expressed in a first approximation as a perturbation series within the Rayleigh-Schrodinger (RS) perturbation theory framework.
The corresponding zeroth-order energy in this case would be the energy of the noninteracting system, that is, E XY 0 : The first order reads as: where m and n label the excited states of X and Y, respectively. E 2 ð Þ ind denotes the so called "induction energy," which is a stabilizing component of the interaction. It represents the energy lowering due to the first-order changes induced in the wavefunctions of the monomers by the electric field generated by the other fragment. It can be understood as the electrostatic interaction between the permanent multipole moments on one fragment and the induced multipole moments on the other (induced by the electric field of the other fragment). Finally, E 2 ð Þ disp is the dispersion interaction between the monomers. This term is purely a quantum mechanical effect, always stabilizing, which is responsible for the stability and structure of a large variety of aggregates in gas and condensed phases, and whose importance in the various fields of chemistry is currently the subject of extensive research. 13,14 It can be understood as the interaction between instantaneous dipole moments forming on one fragment with the induced dipoles forming on the other, 15 and it is a pure dynamical electron correlation effect.

| Symmetry Adapted Perturbation Theory
Several strategies have been suggested to deal with the antisymmetry problem in perturbative approaches. In particular, the Symmetry Adapted Perturbation Theory (SAPT) 16 of Szalewicz and coworkers has proven particularly successful for the description of weak NCIs, as documented in many reviews. [7][8][9] In this approach, each perturbation energy term in the RS expansion (e.g., the electrostatic energy E 1 ð Þ el ) appears together with a repulsive "exchange" term (e.g., E 1 ð Þ exch ) representing the corresponding correction to the energy due to the antisymmetry requirement of the wavefunction. E 1 ð Þ exch is typically the dominant repulsive term in the SAPT energy expression. For instance, at first order we have: Several variants of SAPT exist. In its simplest formulation (SAPT0),Ĥ X andĤ Y are replaced by the corresponding Fock operators of the monomers and the expansion is truncated at the second order for the energy. The resulting interaction energy reads: exch − disp are the second-order exchange corrections. It is worth stressing here that, in the SAPT0 framework, E 1 ð Þ el represents the electrostatic interaction between the unperturbed electron densities of the monomers, which are kept at the Hartee-Fock (HF) level. This limits the accuracy of the approach when electrostatically bound complexes are considered. Analogously, E 2 ð Þ disp is the dispersion interaction between the monomers, which is essentially of second-order Møller-Plesset (MP2) quality. Finally, δHF is defined as: where ΔE HF int is the interaction energy computed at the HF level using a supermolecular approach. This term is often included in the SAPT0 energy expression as it typically improves the quality of E SAPT0 int with respect to benchmark CCSD(T) values, especially for polar molecules. It can be considered as a correction for the higher order terms. Although the accuracy of SAPT0 relies to a large extent on error cancellation, 17 recent implementations of this technique making extensive use of the density-fitting approximation have reduced the scaling of the method from N 7 to N 5 , thus making it applicable to systems of moderate size (up to 200 atoms and 2,800 basis functions). [18][19][20] More sophisticated SAPT approaches can be defined in which intramolecular correlation corrections to the various SAPT0 terms are included via the use of correlated wavefunction-based methods or by using density functional theory (DFT) for the description of the monomers, as discussed in the following section.

| Fifty shades of Symmetry Adapted Perturbation Theory
In the popular SAPT2 approach, a second-order correction with respect to the intramonomer electron correlation is added to the SAPT0 electrostatics, exchange and induction components, yielding interaction energies of essentially MP2 quality. 17 More accurate results can be achieved by including higher-order terms or by using more accurate wavefunction-based methods for the intramonomer electron correlation corrections. 17 For instance, Korona and coworkers have included CCSD intramonomer correlation corrections for all first and second-order terms, [21][22][23][24][25][26][27][28] leading to an accurate benchmark method with formal N 7 scaling.
In addition to the use of wavefunction-based methods, an effective approach for incorporating intramolecular electron correlation relies on DFT for the description of the monomers. Such DFT-based SAPT approaches are denoted as DFT-SAPT or SAPT(DFT). [29][30][31][32][33][34][35][36] The most recent implementations of these methods have achieved remarkable accuracy in the context of NCIs, 37 approaching CCSD(T) results for NCI energies for the popular S66 benchmark set of Hobza and coworkes. 38 The formal scaling of the technique is N 7 . However, similarly to SAPT0, its scaling can be reduced to N 5 using density-fitting approximations.
It is worth mentioning here that, although the SAPT equations for multiple fragments are known, 39,40 SAPT in its standard formulation is rarely applied to the case of more than two fragments. Recently, an approximated low scaling (N 3 ) SAPT variant, called XSAPT, has been introduced that allows for the definition of multiple monomers. [41][42][43] In this scheme, products of monomer wavefunctions computed at the XPol 44 level are used as zeroth-order wavefunctions in the SAPT treatment. Moreover, the SAPT dispersion and exchange-dispersion terms are replaced by an empirical force field term fitted from ab initio SAPT data for small molecules.
Finally, it is worth mentioning that open-shell SAPT variants have also been developed. 20,45,46 However, to the best of our knowledge, they have been applied only to the study of small model systems.

| Limitations of perturbative approaches
Wavefunction-and DFT-based SAPT approaches are among the most widely used computational methods for the study of weak NCIs. The success of these schemes relies on their remarkable accuracy and on the straightforward chemical interpretation associated with the various terms of the perturbation. In particular, the SAPT definition of the London dispersion energy is considered to be particularly reliable.
Unfortunately, being based on perturbation theory, the range of applicability of these schemes is limited to weak interactions. 47 Hence, they must be used with caution when dealing with moderately strong NCIs 7-9 and are not applicable to intermolecular interactions with a significant "covalent" contribution, for example, transition metal-ligand bonds in organometallic chemistry.
Another related, although less significant, limitation of perturbative approaches is that they do not naturally provide a quantification of the charge transfer (CT) contribution to the interaction energy, which is included in the induction energy. 48 It might be argued that CT is quantum mechanically an ill-defined quantity, which cannot be separated unambiguously from the induction energy. Although this argument appears as physically sound, one cannot deny the fact that CT is a fundamental staple of chemical reasoning, and it is often considered as one of the key interaction components in several strong and weak classes of NCIs. [49][50][51] Outside the realm of NCIs (and thus outside the range of applicability of SAPT), the role of CT is considered even more important. For instance, the concept of donor-acceptor interaction is used in essentially all inorganic chemistry textbooks to explain the properties of organometallic complexes. 52 Hence, given its ubiquitous nature, approaches have been suggested to define a CT contribution within the SAPT framework. 53,54 Finally, SAPT can only be used to compute interaction energies between molecules frozen in the geometry they have in the complex. Hence, it neglects the energy required for distorting the fragments from their equilibrium geometry to the ones they have in the adduct. This quantity is denoted here as "geometric preparation", ΔE geo − prep (also known as "strain energy" or "deformation energy"). In typical SAPT applications, ΔE geo − prep is computed using a standard quantum electronic structure method and then added to the SAPT interaction energy.

| GENERAL ASPECTS OF SUPERMOLECULAR APPROACHES
In addition to perturbative schemes, supermolecular approaches can be used to quantify and analyze intermolecular interactions. In this context, a quantum electronic structure method is used to compute the energy of a molecular adduct and that of its constituting monomers at their equilibrium geometry. Then, the binding energy is computed as that of the adduct minus that of the monomers.
In the case of a complex (XY) of two monomers (X and Y), the association energy can be written as: where E XY is the energy of the adduct and E X and E Y are that of the monomers. ΔE can be then divided into a geometric preparation ΔE geo − prep (see Section 2.4) and an interaction energy term ΔE int : The latter describes the interaction between the monomers in their in-adduct geometries, as in SAPT approaches. Essentially, any quantum mechanical method can be used for the purpose of computing intermolecular interactions within a supermolecular approach, going from fast and approximate semiempirical approaches 55 to accurate and typically computational expensive ab initio methods. The optimal choice depends on the interaction type, the accuracy target, and the system size.
In general, the main advantage that supermolecular approaches retain over perturbative schemes is their extended range of applicability. Moreover, several EDA schemes have been developed for decomposing ΔE int computed using a supermolecular approach, 47,56 thus providing a SAPT-like partitioning of the interaction energy. These schemes build upon the original works of Morokuma and coworker 57 and Ruedenberg. 58 In the following, the most commonly applied approaches for the quantification and analysis of NCIs within a supermolecular approach are examined.

| MEAN-FIELD APPROACHES
Mean-field approaches like HF or Kohn-Sham DFT (abbreviated as just "DFT" hereafter) rely on the solution of an independent particle problem in which the electron-electron interaction is replaced by the interaction with an effective potential. Due to a good balance between accuracy and computational cost, these schemes have found widespread applications for the study of intermolecular interactions.

| Accuracy of mean-field approaches for noncovalent interactions
DFT is perhaps the most popular method used in mainstream computational chemistry for the quantification and analysis of NCIs. Although it is well known that it does not describe long-range dispersion interactions properly, 59,60 several effective strategies have been developed to deal with this shortcoming. In the context of this review, these schemes are denoted as "dispersion-corrected" DFT approaches. For a thorough comparison of different correction schemes, the reader is referred to the recent review of Grimme et al. 61 For instance, a commonly applied strategy relies on the use of flexible exchange-correlation functionals, whose many parameters are determined by fitting on high quality benchmark databases. The popular Minnesota functionals developed by Truhlar and coworkers are a prominent example and have proven particularly successful for the quantification of NCI energies. [62][63][64][65][66][67] Alternatively, an ad hoc atom-pairwise force field term can be added to the DFT energy to incorporate the part of the London dispersion that is not described by the exchange correlation functional. 68 This strategy is adopted in the D2/3/4 69-71 methods of Grimme and in the Tkatchenko-Scheffler method, 72 as well as in many others. 61 In the simplest case, the atompairwise correction to the energy can be expressed as: where, for each pair of atoms A and B, R AB is the internuclear distance and C AB 6 is obtained from the dynamic dipole-dipole polarizabilities of the atoms, for example, using the Casimir-Polder formula. As the exchange-correlation functional already accounts for short-range correlation effects, distance-dependent damping functions f(R AB ) are included to avoid double counting of dispersion interactions. The damping functions have different parametrizations depending on the functional. Depending on the scheme, three-body corrections and dipole-quadruple corrections can also be included.
These simple yet powerful approaches allow for the calculation of NCI energies at the DFT level with an impressive accuracy for standard benchmark sets. 73 As an example, the performance of some popular DFT-D3 functionals on the S66 benchmark set is illustrated in Figure 1 (the data were taken from Reference 73). Although standard, "uncorrected" functionals show various performances and large errors, upon the inclusion of the D3 correction all functionals feature a similar Mean Absolute Deviation, that is, less than 0.5 kcal/mol. These results highlight the effectiveness as well as the functionaldependent nature of the D3 correction.
As a note of caution, one should be aware that the remarkable accuracy achieved by dispersion-corrected mean field approaches for standard benchmark sets 73,75 is the result of careful parametrization using available high quality benchmark values on relatively small systems, 76 typically obtained from accurate first principles calculations. Such methods include the "gold standard" CCSD(T) 77 and its local variants as well as Quantum Monte Carlo approaches. 78

| Basic principles of mean-field EDA schemes
An advantage of mean-field approaches for the study of NCIs is that the "wavefunction" in this case is a single Slater determinant, whose chemical interpretation in terms of its constituting molecular orbitals (MOs) is somewhat straightforward. Moreover, a large number of EDA schemes can be used to decompose interaction energies calculated with these approaches into various interaction types. Although the original scheme of Morokuma was based on the HF method, modern EDA schemes provide a decomposition of interaction energies computed using DFT.
A thorough comparison of EDA approaches is beyond the scope of the present review. However, some general features can be identified and discussed that are common among the various techniques. In many popular Morokuma-like EDA approaches, such as the approach of Bickelhaupt and Baerends, 79 the interaction energy can be written as: The first term, ΔE el represents the electrostatic interaction between the (typically) unperturbed densities of the monomers. It is conceptually and numerically 56 similar to the electrostatic term defined in SAPT approaches.
The ΔE Pauli term represents the so called "Pauli Repulsion," which essentially describes the repulsive component of the interaction due to the overlap of the densities of the interacting monomers. It is usually calculated as the difference between the energy of a "frozen state," in which polarization and CT effects between the fragments are absent, and the energy of the isolated monomers, minus the electrostatic interaction ΔE el . Although various definitions for the frozen state are possible, 80 in the most common scenario it is defined as the antisymmetrized product of the wavefunctions of the fragments in their in- The ΔE orb is denoted as "orbital interaction," and describes a stabilizing component of the interaction associated with CT and polarization effects. It is obtained as the energy difference between the SCF solution for the adduct, that is, the DFT energy, and the energy of the frozen state. Hence, it is intuitively similar to the DFT-SAPT induction term.
Finally, E disp is the dispersion energy, which is typically estimated from the magnitude of the D correction (if used) or with the aid of an auxiliary "dispersion-free" functional, as detailed in the next section.
Many different variants of the abovementioned approaches exist. For instance, ΔE orb can be further divided into orbital-pairwise contributions using the extended transition state (ETS) 81 method of Ziegler in conjunction with the natural orbital for chemical valence (NOCV) 82 approach of Mitoraj and Michalak, which is of particular relevance in the analysis of coordinate covalent bonds in transition metal chemistry. The approach is typically denoted as the ETS/NOCV method. 83 Other popular mean-field approaches include the block-localized wavefunction EDA, 84,85 and the closely related Absolutely localized molecular orbital (ALMO)-EDA scheme of Head-Gordon and coworkers. 86,87 Instead of an overall ΔE orb term, these schemes provide a separate quantification of polarization and CT contributions. Essentially, this is achieved by defining an additional intermediate state in which the fragment MOs are allowed to relax variationally in the presence of the other fragment but with the constraint of prohibiting MO delocalization between fragments. Correlated variants of the ALMO-EDA approach have been also defined, as detailed in the next sections.
Alternatively, electron density-based approaches have been developed such as the Natural Energy Decomposition Analysis (NEDA) 88,89 in the natural bond orbitals framework of Weinhold and coworkers, 90 or the Interacting Quantum Atom (IQA) method, 91,92 which relies on the Quantum Theory of Atoms In Molecules. 93 In principle, accurate electron densities obtained from correlated calculations can be used as input in these approaches. 91

| Strengths and limitations of mean-field approaches
The efficiency and the accuracy of mean-field approaches for the analysis of NCIs is determined by that of the underlying exchange-correlation functional. On the one hand, these schemes are typically much faster than perturbative approaches and are thus applicable to larger systems with hundreds of atoms. On the other hand, their main limitation is the heavy reliance on benchmark studies for parametrizing, assessing the accuracy and in some cases re-parametrizing the exchange-correlation functional. Considering that the large majority of reliable benchmark studies for NCIs focus on small systems with a few tens of atoms, one might question the reliability of these methods on computational predictions for large systems. 76 Fortunately, with the development of fast and accurate first principles methods like the modern local CCSD(T) variants, it is becoming possible to put these methods to a quantitative test for large systems. 73,94 As regards the interpretation of the various terms of the decomposition, the strongest point of mean-field EDAs is a chemically intuitive definition of CT and polarization effects (e.g., the ΔE orb term in the Bickelhaupt and Baerends scheme), as shown for instance by several studies of Frenking and coworkers 95,96 and many others. 47,79,97,98 However, an important limitation concerns the definition of the dispersion energy. In the large majority of cases, the magnitude of the force field term, for example, the D correction, is taken as a practical estimate of its value. Although this strategy has proven extremely useful for trend studies, 99,100 it has two practical shortcomings: 1. Being independent of the details of the electronic structure of the system, many force-field corrections are uniquely determined by its geometry. Although some recent approaches like D4 71 and the Tkatchenko-Scheffler method 72 incorporate density-dependent terms in the dispersion energy expression, it is clear that a force field term is not the ideal tool for describing charge, spin or electronic state-dependent effects on the dispersion energy. For instance, ions have different polarizabilities than their neutral counterparts and thus give rise to different dispersion interactions when interacting with another atom. Force field corrections like D3 are "blind" to such charge effects. The most rigorous way to address these problems is to rely on accurate first principles methods that naturally provide an electronic structure-dependent definition of London dispersion, as discussed in Section 6.4. As an example, recent DLPNO-CCSD(T)/LED calculations have elucidated recently the influence of the spin state on the strength and nature of various intermolecular interactions, 4,101 including hydrogen-and halogen-bonding interactions in Van der Waals complexes of singlet and triplet carbenes as well as Fe(II)-ligand interactions in heme derivatives. 2. As already mentioned, force field corrections do not account for the part of the dispersion energy that is already described by the exchange correlation functional, that is, short-to mid-range London dispersion effects. 68 Hence, to avoid double counting, damping functions are used in Equation (10) that are parametrized specifically for each functional. This makes the overall correction functional dependent, as shown in Figure 1, making the quantification of the relative importance of London dispersion on intermolecular interactions particularly difficult in those cases when short-to mid-range dispersion effects are important. As an illustrative example, consider the interaction between a bulky Lewis acid and a bulky Lewis base in the frustrated Lewis pair (FLP) shown in Figure 2. The inclusion of the D3 correction allows to obtain association energies with near-DLPNO-CCSD(T) accuracy for all functionals tested (the basis set superposition error (BSSE) corrected DLPNO-CCSD(T) value at the complete basis set limit (CBS) is −10.6 kcal/mol, as reported in Reference 102). However, the magnitude of the D3 correction ranges from −8.8 to −18.7 depending of the functional, that is, from the 83% to the 176% of the DLPNO-CCSD(T)/CBS energy.
F I G U R E 2 Binding energies for a frustrated Lewis pair computed using various DFT and DFT-D3 functionals. The Lewis acid is B(C 6 F 5 ) 3 , the Lewis base is the phosphine P(Mes) 3 . The contribution from the D3 correction is shown with a red arrow for the two extreme cases. The blue horizontal line represents the benchmark DLPNO-CCSD(T) value at the complete basis set limit (CBS). Data were taken from Reference 102 Finally, it is worth mentioning that another strategy for quantifying London dispersion within a mean field approach has been recently introduced by Head-Gordon and coworkers in the context of the ALMO-EDA scheme. For each exchange correlation functional, a corresponding ad hoc "dispersion-free" functional is initially defined (e.g., the HF). Then, the dispersion is approximately estimated as the part of the inter-fragment exchange-correlation energy which is not described by the dispersion-free functional. 87 The limitation of this approach is that the reliability of the dispersion estimate depends on both the accuracy of the exchange-correlation functional in describing the system and on the ability of the auxiliary dispersion-free functional to capture all other aspects of the interaction besides London dispersion.

| MP2-BASED APPROACHES
MP2 is the simplest useful post-HF method for incorporating dynamic electron correlation. Unlike mean field approaches, it naturally describes London dispersion forces, even though it typically overestimates NCI energies for dispersively bound systems. Hence, MP2 generally outperforms HF and "dispersion-uncorrected" exchange-correlation functionals in the calculation of interaction energies for noncovalently bound systems. Importantly, MP2-based EDA schemes have set the stage for the development of more sophisticated EDA approaches based on accurate wavefunction-based methods, for example, CCSD(T). Basic principles, strengths, and limitations of MP2-based EDAs are reviewed in the following sections.

| Basic principles
The energy obtained from any correlated wavefunction-based method can be written as a sum of a reference (E ref ) and correlation (E C ) contributions: where E ref is typically the HF energy E HF in the closed shell case. Hence, the total interaction energy in Equation (9) can also be written as a sum of reference and correlation contributions: Mean field approaches of the type described in Section 4.2 can be used to decompose ΔE ref int into chemical meaningful terms, for example, as in Equation (11). However, the development of EDA schemes for the correlation part of the interaction energy (ΔE C int ) has proven much more challenging. As London dispersion is a pure correlation effect, it is customary to use ΔE C int as a practical estimate of its magnitude. This strategy is adopted in various EDA schemes, such as the popular LMO-EDA method of Su and Li. 103 However, this choice is not satisfactory because ΔE C int incorporates a number of different physical effects besides London dispersion, 3 as it will be discussed in detail in Section 6.4.
A useful starting point for any correlated EDA scheme is the observation that the correlation energy of a closed-shell system can always be written as: where i, j denote the occupied orbitals, a, b the virtual orbitals andτ ij ab the contravariant amplitudesτ ij ab = 4τ ij ab − 2τ ji ab 1 + δ ij . Note that, in the coupled cluster theory, τ ij ab = t i a t j b + t ij ab , where the doubles t ij ab and singles t i a amplitudes are obtained by solving a set of nonlinear equations. Hence, Equation (14) shows that the correlation energy can always be expressed as a sum of double excitation contributions from pairs of occupied orbitals to the virtual space. Moreover, it follows from Equation (14) that E c can be partitioned into pair correlation energies ε ij : Equation (15) has several practical applications. In local correlation schemes, this partitioning is used to develop efficient local approximations to the correlation energy. [104][105][106][107][108][109] For instance, upon localization of the occupied space, for example, by using the Foster-Boys scheme, 110 it can be demonstrated that ε ij decays very fast with the distance between the centroids of i and j. Hence, many pair correlation energies can be neglected or treated at a low level of theory without significant loss of accuracy. Importantly, Equations (14) and (15) have proven extremely useful for interpretative purposes, as discussed in the next section.

| MP2-based EDAs
By assigning each orbital to the fragment in which it is dominantly localized, the correlation energy can be divided into intraand inter-fragment contributions using Equation (15). As an example, this simple approach was used by Grimme et al. 111 to study the interaction between aromatic systems at the Local MP2 (LMP2) level. The resulting inter-fragment correlation energy was also divided into contributions from pairs of different orbital symmetry, that is, σ-σ, σ-п, and п-п. An analogous approach is used in the framework of the fragment molecular orbital approach of Kitaura et al. [112][113][114] Unfortunately, this strategy only provides limited insights as it does not allow for discriminating between dispersive and CT excitations, leading to an unsatisfactory definition of London dispersion. 115 However, if a local representation of the virtual space is used, virtual orbitals can also be assigned to the fragment in which they are dominantly localized (i.e., X or Y). This allows regrouping the double excitations in Equation (14) into various families, corresponding to the different physical components of the interaction. This strategy was initially suggested by Schütz et al. using the LMP2 method. 115 In this case, virtual orbitals were obtained from the atomic orbitals by projecting out the occupied space. A graphical representation of the various excitations that can be defined in the presence of two fragments is shown in Figure 3. It is worth mentioning here that CT excitations of the type (e) and (f) were not included in the decomposition and were regarded as BSSE excitations. Modern local EDAs, 2-4 which rely on a different definition of the virtual space, incorporate explicitly these excitations in the decomposition, as it will be detailed in Section 6.2.
The approach of Schütz et al. has been applied to the study of metal-metal interactions in small model systems of coinage metals. 116 In the same study, a CCSD variant of the scheme was also employed. Moreover, the scheme has found applications in the study of aurophilic interactions. 117,118 The LMP2 dispersion energy can be also visualized by means of the so called dispersion interaction density plot introduced by Mata and coworkers. 119 Using an analogous partitioning of the correlation energy, Head-Gordon and coworkers proposed an energy decomposition scheme for the decomposition of MP2 quality energies within the ALMO-EDA framework. 120,121 Importantly, their method can be used to identify a series of chemically reasonable components in ΔE C int , including attractive as well as repulsive contributions. So far, this scheme has been applied for the study of intermolecular interactions in small model systems, 122 including halogen-bonded complexes. 123

F I G U R E 3 Classes of double excitation
contributions in which the correlation energy can be decomposed in post-HF EDA schemes. X and Y denote the fragments

| Limitations of MP2-EDA schemes
The strongest point of MP2-EDA schemes is that they naturally describe long-range correlation, that is, London dispersion. However, their main shortcoming is their limited accuracy and range of applicability. In fact, MP2 is known to overestimate dispersion interactions, it describes poorly open-shell systems and fails dramatically for strong interactions, especially when transition metals are present. To overcome these limitations, coupled cluster-based EDA schemes have to be used. These schemes are reviewed in the next section.

| COUPLED CLUSTER ANALYSES OF NONCOVALENT INTERACTIONS
In this section, the basic principles and selected applications of a recently introduced EDA scheme for the decomposition of NCI energies computed with local variants of CCSD(T) are reviewed. In the application part, emphasis is placed in discussing the role of the all-important London dispersion interaction, whose quantification has proven challenging for supermolecular approaches.

| Accuracy and efficiency of CCSD(T)-based schemes
High level correlated wavefunction-based methods can be used to compute NCI energies with great accuracy, provided that dynamic electron correlation is properly taken into account. 9,124,125 In this context, CCSD(T) is known to provide the most accurate energies and properties among the practical single reference ab initio methods applicable to systems of moderate size. 77 Unfortunately, like in standard SAPT, the formal scaling of the CCSD(T) method in standard implementations is N 7 . This has stimulated the development of approximations aimed at reducing the computational cost of this technique without degrading its accuracy. In particular, by exploiting the short-range nature of the dynamic electron correlation, 104-109 local variants of the CCSD(T) method have been developed 5,6,126,127 that allow for the calculation of CCSD(T) quality energies for systems with hundreds of atoms. In fact, benchmark studies by Neese and coworkers 128 as well as by Werner and coworkers 129 have shown that modern local CCSD(T) methods, 5,6,126,127 in which the virtual space is spanned by pair natural orbitals (PNOs), 130 can be used to compute NCIs with essentially canonical-coupled cluster accuracy.
In particular, the DLPNO-CCSD(T) method 5,6,126 developed by Neese et al. allows one to compute NCI energies with an accuracy of~0.2 kcal/mol on standard benchmark sets with respect to its canonical counterpart, provided that TightPNO settings are used. 130 Recently, DLPNO-CCSD(T)/TightPNO calculations in conjunction with large basis sets have been applied to the study of large systems and/or with complicated electronic structure, 4,102,131,132 including for instance quantum mechanical/molecular mechanical (QM/MM) calculations with a QM region of more than 300 atoms and 6,800 basis functions. 132 Moreover, the recently introduced LED scheme 2-4 can be used to decompose DLPNO-CCSD(T) interaction energies into chemically meaningful contributions for virtually any system within the range of applicability of the coupled cluster method. As an example, DLPNO-CCSD(T)/LED/TightPNO calculations at the CBS limit were used to investigate the intermolecular interaction between the Lewis acid and the Lewis base in 14 FLPs and the results were used to propose designing principles for FLPs with tailored bonding features and reactivity (4,500 basis functions were used for the largest system in the set). 102 Hence, the DLPNO-CCSD(T) method in conjunction with the LED scheme is a powerful tool for the accurate quantification and analysis of NCIs in large systems. Moreover, as it will be discussed in the following, the scheme is not limited to the study of weak interactions and is applicable to the simultaneous analysis of multiple fragments and for closed-shell 2 as well as for open-shell 4 systems. The theory underlying the LED scheme is reviewed in the next section and selected applications are examined.

| Local energy decomposition of Total DLPNO-CCSD(T) energies
In contrast to the large majority of the schemes discussed in the present review, the LED analysis 2-4 provides a decomposition of the total energy of a system, computed at the DLPNO-CCSD(T) level, into an arbitrary number of fragment-pairwise contributions. Using this approach, DLPNO-CCSD(T) energies can be decomposed exactly into additive components, irrespective of the nature of the system and of the strength of the interaction between the fragments. For instance, this partitioning can be used to discuss intramolecular interactions 133 or, more generally, any interaction for which single point energies on the isolated fragments are not feasible. Above all, the LED scheme has found widespread applications in the analysis of strong and weak intermolecular interactions within a supermolecular approach. 99,101,102,[133][134][135][136][137][138] Consider our prototype complex (XY) consisting of two interacting fragments X and Y. The DLPNO-CCSD(T) energy for the adduct can be written as a sum of a correlation E C and a reference contribution E ref . Typically, E ref is the HF energy in the closed shell case 2 while it is the energy of the quasi-restricted orbital (QRO) 139 determinant in the open shell case. 4 By exploiting the fact that occupied orbitals are localized in the DLPNO-CCSD(T) framework, one can assign each orbital to the fragment in which it is dominantly localized. This allows us to exactly decompose E ref as the sum of four terms: where E The dominant DLPNO-CCSD(T) part of the correlation energy, E C − CCSD , can be decomposed into fragment-pairwise terms using a partitioning similar to that of Figure 3, 2-4 while the contribution from the perturbative triples correction is decomposed into intra-and inter-fragment terms based on the localization of the occupied orbitals. This leads to the following expression for the correlation energy: represents the contribution from the intra-fragment correlation, E C− CCSD disp is the DLPNO-CCSD (T) dispersion interaction between the fragments, E C− CCSD CT is the contribution from CT excitations between the fragments (including all CT excitation types shown in Figure 3), and E inter C − T ð Þ is the inter-fragment correlation contribution from the triples.
In the case of more than two fragments, a delocalized correlation contribution, which is typically negligible, is also included. It collects the contribution from the delocalized double excitations shown in

| Local energy decomposition of DLPNO-CCSD(T) interaction energies
The LED scheme can be used to decompose interaction energies ΔE int between two monomers using a supermolecular approach. By subtracting from the decomposed DLPNO-CCSD(T) energy discussed in the previous section the ground state reference and correlation energies of the monomers in the geometry they have in the dimer, the following expression for ΔE int is obtained: involving more than two fragments where ΔE ref el − prep is called electronic preparation and corresponds to the energy required to distort the electron density of the monomers from their ground state to the one that is optimal for the interaction, that is, the difference between E ref and the ground state energy of the monomers at the reference level. E C− CCSD disp is the DLPNO-CCSD(T) dispersion contribution. ΔE C− CCSD no − disp is the non-dispersive part of the DLPNO-CCSD(T) correlation contribution to the interaction energy and incorporates all CT and intra-fragment excitations. Its physical meaning is discussed in detail in the next section. ΔE C − T ð Þ int is the triples correction contribution to the interaction energy. Note that a strategy has been recently suggested to further decompose ΔE C− T ð Þ int into dispersive and non-dispersive terms. 136 It is worth mentioning that all reference terms, that is, ΔE ref el − prep , E ref elstat and E ref exch , are obtained from the fully relaxed localized MOs of the adduct. Hence, they incorporate permanent as well as orbital relaxation components to the interaction energy. A strategy to disentangle these contributions was recently proposed in Reference 3. Accordingly, each of these terms is further decomposed into a "frozen" and "orbital relaxation" contribution. For instance, the frozen component of E ref elstat describes the permanent electrostatic interaction between the monomers, while the associated orbital relaxation contribution corresponds to the induction energy. Similarly, the frozen component of ΔE ref el − prep is conceptually similar to the "Pauli repulsion" term in mean-field EDAs, while the corresponding orbital relaxation term describes the energy investment associated to the polarization of the fragments.

| LED study of prototype interactions
As London dispersion is a pure electron correlation effect, it might be argued that ΔE C int already provides a reasonable estimate for its magnitude. In fact, this strategy is adopted in many popular energy decomposition schemes. 103 Unfortunately, this definition of London dispersion is unsatisfactory because electron correlation contains a number of physically different contributions besides London dispersion. In fact, it also provides a correction for all the other interaction components that are approximately described at the HF level. As an example, Head-Gordon and coworkers observed that ΔE C int assumes positive values in the long-range region for the interaction of two water molecules. 140 This effect was attributed to the characteristic overestimation of dipoles at the HF level, which leads to an overestimation of the electrostatic attraction at the reference level.
One of the main features of the LED scheme is a clear-cut separation between dispersive and non-dispersive contributions to the interaction energy at the coupled cluster level. In order to discuss this aspect, we have examined in Reference 3 the nature of the correlation energy in different interaction types. Our study encompassed a series of systems held together by intermolecular interactions of different nature, such as the Ar dimer, Ar─Li + , the beryllium dimer and the water dimer. The HF and DLPNO-CCSD(T) binding energy profiles for these systems are shown in Figure 5, together with a decomposition of the correlation contribution at the minimum geometry into dispersive and non-dispersive terms. The available experimental data [141][142][143] for the binding energy are also shown with a red cross.
DLPNO-CCSD(T) and experimental binding energies were found to be in striking agreement in all cases. Moreover, the binding energy profiles demonstrate that electron correlation describes entirely different effects in different systems. For the dimers of Ar and water, ΔE C int is dominated by the dispersion component (in the water dimer case the situation changes in the long range, when electrostatics dominate the interaction; for the Ar dimer dispersion dominates the correlation for all interatomic distances, as discussed in Reference 3). Interestingly, a numerical comparison between the LED and DFT-SAPT estimates of the London dispersion energy for the water dimer case demonstrated that the two methods converge to the same values in the long range, while they slightly differ in the short range. 134 In contrast to the abovementioned cases, London dispersion constitutes only the 8% of the overall ΔE C int in ArLi + . The large non-dispersive contribution in this case originates from the characteristic underestimation of polarizabilities at the HF level of theory, which leads to an underestimation of the induction component at the HF level. Finally, in the Be dimer case, the correlation contribution to the binding energy, which is responsible for the stability of the adduct, incorporates dispersive and nondispersive effects in equal proportions. In this case, the non-dispersive correlation provides a correction to the "electron sharing" phenomena, which is approximately accounted for at the HF level. This effect was originally described by Ruedenberg 58 and leads to a substantial accumulation of charge in the middle of the bond between the two beryllium atoms.
Finally, it is worth mentioning that all the above considerations were found to be consistent with the detailed analysis of the charge displacements occurring in the DLPNO-CCSD(T) electron density upon bond formation. Moreover, the dispersion energies obtained from the LED scheme are in qualitative agreement with those obtained using the simple but physically meaningful London formula. 144

| London dispersion effects in organometallic chemistry
London dispersion is known to play a fundamental role in organometallic chemistry as demonstrated by a large number of experimental and theoretical works. 13,99,100 However, the quantification of the London dispersion component of the interaction between a transition metal (TM) and a coordinating ligand (L) has proven to be a challenge for theoreticians. In fact, SAPT was designed specifically for the study of weak NCIs and it is thus not applicable for the study of moderate or strong intermolecular interactions involving transition metals.
Grimme and coworkers have recently compared the D3 values obtained with various functionals for a series of intermolecular interactions in various organometallic complexes. They observed that, if functionals like BLYP or B3LYP are used, 99 the magnitude of the D3 correction is close to the LED estimate of the dispersion energy in a large variety of cases. These findings are consistent with previous results from our group. 102,135 However, the agreement between D3 and LED dispersion is of course not universally valid and there are notable exceptions. 136 Hence, LED appears as a natural choice for analyzing TM-L dispersion forces. As an example of its usefulness in this context, we examine here the case of the interaction between a TM and a CH group in a series of agostic complexes and alkane σ-complexes. These systems are considered as model intermediated of CH activation reactions 145,146 and their importance is demonstrated by the impressive number of experimental and theoretical studies devoted to their characterization. 147 Herein, we discuss the main findings of these studies on two representative systems, shown in Figure 6.
The intermolecular ReÁ Á ÁCH 4 binding energy profile for the experimentally characterized σ-complex [CpRe(CO) 2 (CH 4 )] ( Figure 6, left panel) shows that London dispersion dominates the interaction between the alkane and the transition metal. Remarkably enough, these results hold true for a large number of alkane σ-complexes, demonstrating that this fundamental dynamical correlation component is to a large extent responsible for the thermodynamic stability of these systems. Based on these and on other results concerning the bonding features of these systems, a series of designing principles were suggested for stable σ-complexes. In particular, the magnitude of the London dispersion energy can be modulated by increasing the size of the alkane and of the ancillary ligands on the TM, as discussed in detail in Reference 136. A similar result concerning the fundamental importance of London dispersion for the interaction between the TM and a CH group was found for a series of agostic complexes. 133 As an example, we consider here the [EtTiCl 3 (dmpe)] system, that was the first β-agostic complex to be experimentally characterized by X-ray spectroscopy. 149 In the right panel of Figure 6, the rotational barrier for the rotation of the agostic methyl group around the C α ─C β bond of the ethyl ligand is shown at the DLPNO-CCSD(T) and HF levels of theory and the correlation contribution to the barrier is decomposed into dispersive and non-dispersive terms. These results demonstrate that, without London dispersion, this textbook case of agostic complex would not be stable.

| Analysis of open-shell systems
As already mentioned, only a few EDA [150][151][152][153] and SAPT 20,45,46 schemes are available for the study of NCIs in open-shell species and have found so far limited applications. Hence, the open-shell variant of the LED scheme is particularly promising in this context and has been already applied to the study of a number of different intermolecular interactions. 4,101 In particular, it F I G U R E 6 Left: Binding energy profiles for the ReÁ Á ÁCH 4 interaction in [CpRe(CO) 2 (CH 4 )]computed at the DLPNO-CCSD(T)/def2-TZVPP and HF/def2-TZVPP level of theory as a function of the inter-fragment distance. Data were taken from Reference 136. Right: DLPNO-CCSD(T)/ def2-QZVP HF and HF/def2-QZVP energy profiles associated to the rotation of the agostic methyl group around the C α ─C β bond in the agostic [EtTiCl 3 (dmpe)] (dmpe = 1,2-bis(dimethylphosphino)ethane). The reference energy corresponds to the equilibrium geometry (θ = 0). Data were taken from Reference 133 was used to elucidate the mechanism that governs the change in the singlet-triplet energy gap of carbenes 4,101 and of heme 4 upon formation of a complex.
It is experimentally known that heme, that is, a ferrous iron-coordinated porphyrin (denoted hereafter as PFe), has a triplet ground state. However, its CO-bound PFeÁÁÁCO analogues have a singlet ground state. [154][155][156][157] Thus, CO binding reverses the spin-state ordering of heme derivatives. However, the chemical mechanism behind this experimental observation was not clear. This aspect was investigated in Reference 4 using the open-shell variant of the LED scheme. As a case study, the simplest PFe model without additional axial ligands or side chains was considered.
In Figure 7, the binding energy profiles for the interaction between CO and PFe computed at the DLPNO-CCSD(T) and QRO levels are reported for three different spin states of PFe, namely the 1 A 1g / 3 E g / 3 A 2g states. The DLPNO-CCSD(T)/LED analysis clearly shows that PFe( 3 E g )ÁÁÁCO and PFe( 3 A 2g )ÁÁÁCO adducts can be described as van der Waals complexes stabilized by weak dispersion forces. In fact, their interaction is repulsive at the QRO level and correlation is dominated by London dispersion in both cases. In contrast, the PFe( 1 A 1g )ÁÁÁCO bond can be described as a strong interaction in which dispersion is present but plays an auxiliary role, being only 41% of the correlation binding energy. In fact, PFe( 1 A 1g ) features an empty d z 2 orbital on the central iron atom which is capable of accepting electrons from the CO lone pairs (Fe(II) CO σ-donation). At the same time, the filled d xz and d yz orbitals can donate charge into the π* antibonding orbitals of the CO molecule (Fe(II) ! CO π-backdonation). Hence, the PFe( 1 A 1g )ÁÁÁCO bond can be described as a strong donor-acceptor interaction. The different nature of the interaction between PFe( 1 A 1g )ÁÁÁCO and PFe( 3 E g / 3 A 2g )ÁÁÁCO is thus responsible for the selective spin state stabilization of the singlet state of heme upon CO binding.

| Drawbacks
It might be argued that a possible shortcoming of local approaches for the decomposition of correlation energies is the heavy reliance on localized MOs. In fact, unlike electron densities or total energies, MOs are not quantum mechanical observables. A simple rotation of the occupied space could potentially change the energy components obtained from such decomposition schemes, even in those cases for which the total binding energy is not significantly affected. However, this aspect has been investigated in detail and it has been demonstrated that, in the vast majority of cases, the components of the decomposition are largely independent of the localization scheme employed for the occupied and virtual space, especially for NCIs. 2 Of course, in trend studies, one has to be consistent with the choice of the localization scheme.
However, a significant limitation of the LED scheme is the computational cost associated with the underlying DLPNO-CCSD(T) method. In fact, although DLPNO-CCSD(T) calculations with several hundreds of atoms have been reported, the accurate quantification of NCI energies typically requires the use of relatively large basis sets and TightPNO settings. In certain critical cases, 131 even more conservative PNO thresholds are required in order to achieve quantitative accuracy. Although the controllable accuracy of the method is definitely a strong point of the scheme, these drawbacks limit its applicability to systems with 200-300 atoms if DLPNO-CCSD(T)/TightPNO/CBS accuracy is the target. This aspect is currently under investigation in our group and new efficient LED-based strategies for the quantification and analysis of NCI energies are under development.

| CONCLUSIONS
We have discussed fundamental theoretical aspects, strengths, and limitations of computational methods for the analysis of NCIs in large systems. After introducing the most popular perturbative and supermolecular approaches, we have reviewed recent developments and applications of an EDA scheme based on a local variant of the CCSD(T) method, that is, the LED analysis. For systems of hundreds of atoms, LED can be used to quantify accurately and analyze intra-and intermolecular interactions between an arbitrary number of fragments for both closed-and open-shell systems. The scheme remains valid irrespective of the nature of the monomers and of the strength of the interaction. An important feature of LED is a clear-cut definition of London dispersion, which is numerically consistent with SAPT calculations on weakly interacting systems as well as with simple physical models like the London dispersion formula. While LED is rapidly becoming popular in computational chemistry, there is still room for further progress in the methodology. The main challenge is the development of strategies to further increase the efficiency of the method.