Perspective: Found in translation: Quantum chemical tools for grasping non-covalent interactions

Because of the indisputable importance of non-covalent interactions (NCIs) in both the micro- and macroscopic worlds, it is little wonder that chemists are interested in harvesting the power of these amazing forces. However, thorough exploitation requires understanding not only of the physical origin of the phenomenon but also ways to accurately quantify their magnitude. Gaining such insight is nontrivial. Of course, NCIs do not arise from odd, unpredictable, or mysterious phenomena but from electrical forces acting on charge distributions. A logical entry point to gain insight is examining the charge distribution. This distribution, in particular the electron density, can be accessed through quantum mechanics by approximately solving the Schrödinger equation. Using classical electrodynamics,^1,2 the exact non-covalent forces can be computed from the electron density (Fig. 1) and the exact interaction energy can be extracted directly from the eigenvalue equations. In fact, without relying upon any computational analytical tools (vide infra), an accurate quantum chemical method provides the wavefunction, density, and interaction energy, which are the most comprehensive pieces of information available for describing any interaction.

The question remains, however, does knowing the total interaction energy provide sufficient “insight”? Predictably, this question might be met with widely varying answers, although most chemists would likely respond with a firm “no.” But why? In contrast to other theoretical domains for which a mathematical proof is sufficient, chemists tend to rely upon intuitive concepts, such as dipole-dipole interactions, charge transfer, and polarization to compare and contrast molecular behavior (see Coulson’s quote, Box 1). In this respect, the “energy” represented as a single number along with the mathematically intricate “wavefunction” and the insensitive “electron density” are quite unsatisfactory since they provide little chemical insight. Ultimately, a type of “rosetta stone” is needed that capably translates the results of quantum computations into intuitively useful concepts, which appeal to chemists. Ideally, this tool would treat all possible interactions and molecular systems while being fully self-consistent, thereby providing an answer to any conceivable chemical question. While a single all-encompassing tool would be highly desirable, in reality, the practicing chemist must utilize a number of different tools each having a limited purpose. Some frequently encountered questions might include the following:

1. How do these two molecules interact?

2. Which region of the system is affected by this interaction?

3. Can this interaction be partitioned into easy to understand concepts such as electrostatics/polarization?

4. What is the contribution of each of these ingredients?

These four questions outline the scope of different classes of tools available while also naturally dividing them into two categories: quantitative [including energy decomposition analysis (EDA) schemes] and qualitative ones that permit visualization and identification [such as the quantum theory of atoms in molecules³ (QTAIM), the density overlap region indicator⁴ (DORI), and the non-covalent interaction (NCI) index⁵] (see Fig. 2). Of course, the specific answers to these questions will vary based upon the tool employed. As an example, the concept of “charge transfer” is fairly intuitive (although ill-defined in the complete basis set limit) in the orbital-based picture used by the most variational EDAs. It is defined as the result of excitations of electrons localized on a monomer to a virtual orbital localized on its partner. Since the energy operator does not distinguish between the monomers, that picture does not hold for most perturbation-based EDAs (e.g., SAPT⁶). Likewise, definitions of “dispersion” in tools like SAPT and ALMO-EDA⁷ are non-equivalent. As could be imagined, these are but two examples of the types of problems commonly encountered. As such, it is imperative to understand and define the scope of the computational tools because the underlying object (i.e., density, wavefunction, or Hamiltonian) is often as important as the result itself for gaining knowledge into the nature of a specific non-covalent interaction.

Taken separately, each tool for analyzing interactions can be seen as a type of pocket dictionary which contains a rather limited vocabulary that is only consistent with the object (i.e., the wavefunction, the Hamiltonian, and the density) upon which the tool was constructed. The situation is even more complicated for non-covalent interactions that, as a result of their subtle nature, require additional consideration ensuring that the chosen computational methods are sufficiently accurate. However, it is precisely these two factors, the need for both accuracy and computational efficiency, which makes computational analysis tools so valuable.

Aside from the inherent educational value gained in converting quantum mechanical objects into a chemical language, the tools used to analyze NCIs also play a vital role in constructing new, efficient, and accurate computational methods. It is quite logical to imagine that decomposing a single interaction into separate terms and building a model upon this information is far simpler than to construct the wavefunction of a system from the bottom-up. Indeed, the former idea of modeling an interaction from individual energy terms was the foundation upon which force fields such as UFF,⁸ X-Pol,⁹ AMBER,¹⁰ or  CHARMM¹¹ were built. Likewise, many ab initio methods, including density functional theory (DFT), also benefit from similar treatments, with a noteworthy example being dispersionless density functionals corrected for dispersion (dlDF + D).¹² Computational methods constructed using this philosophy include specialized^13,14 or general⁸ force fields as well as electronic structure approaches targeting a class of interactions (e.g., dlDF + D for dispersion-driven complexes).¹² 

The delineation of those different classes itself relies upon the utilization of these tools. For example, the use of energy decomposition schemes clearly identifies that water-water interactions as dominated by electrostatics, alkane chains interact with one another through London dispersion interactions, and certain molecular pairs are likely to form charge-transfer complexes. Thus, the energy decomposition schemes allow chemists to classify systems based on a clearer set of criteria than just geometrical patterns and functional groups, which further facilitates the development of tailored computational methods discussed earlier (e.g., force fields for water or charge-transfer targeted density functionals). Finally, creating links between properties and different types of interactions can assist in the design of new materials including catalysts, nanowires, superconductors, semiconductors, and pharmaceuticals to name a few. Maximizing desired properties such as conductivity, regioselectivity, and yield becomes a much simpler task when working with a well-defined descriptor like dispersion or charge transfer energies. This illustrates the third goal of those tools (particularly those used to analyze non-covalent interactions) that is to move beyond the simple goal of gaining additional chemical insight. Overall, these fundamental tools can now be used for very practical purposes such as developing some novel computational methodologies and rationally designing new materials. The objective of this perspective article is to provide a map that facilitates the navigation of the methodological toolbox with particular emphasis placed on identifying the best methodology for a specific purpose.

As discussed in the Introduction, the wavefunction provides the maximum possible amount of information about a system. To this end, the wavefunction is the central concept in most variational Energy Decomposition Analysis schemes. Because systems interact to lower their overall energy, variational EDA schemes associate each specific interaction with a lowering of energy due to a certain aspect of wavefunction relaxation. The oldest scheme of this type, the Kitaura-Morokuma (KM) EDA,^16,17 employs a Hartree-Fock wavefunction constructed from doubly occupied (assuming a closed-shell system) orbitals. Beginning from the dimer, a computation begins with four sets of orbitals, an occupied and a virtual set on each monomer. Relaxing the dimer wavefunction using progressively larger orbital sets introduces interactions between the monomers. Specifically, mixing only the occupied or only the virtual orbitals gives rise to electrostatic and exchange interactions, whereas mixing virtual and occupied orbitals together introduces induction and charge transfer (CT) terms (Fig. 3). To separate exchange and electrostatic interactions, KM EDA (as well as some later developed EDAs including CAFI¹⁸ and PIEDA¹⁹) uses the wavefunctions that are not fully antisymmetric, which can lead to the overestimation of the induction energy²⁰ that arises from this physically unjustified choice. Related as well as competing approaches including the Extended Transition State (ETS²¹ and NOCV-ETS^22,23) approach, Bickelhaupt-Baerends EDA,^24,25 RVS,²⁶ CSOV,²⁷ BLW-EDA,^28–30 singles-CI EDA,³¹ de Silva-Korchowiec EDA,³² and ALMO-EDA⁷ have addressed this issue by choosing not to separate the exchange from the electrostatic term. Unfortunately, by not distinguishing the exchange from the electrostatic terms, these methods are hindered in describing problems involving steric clashes. However, for non-covalent interactions, the exchange interaction problem becomes less important as it primarily describes interactions that are short-range in nature.

An alternative approach to separate the exchange and Coulomb contributions is using a DFT wavefunction in which the definition of exchange follows naturally from the exchange component of the density functional, as is done in LMO-EDA,³³ GKS-EDA,³⁴ and steric-based EDA³⁵ and a procedure involving projection operators introduced by Head-Gordon and co-workers.³⁶ 

Kitaura and Morokuma assert that the wavefunction contains the maximum possible amount of information, which should be extracted using orbitals. Most variational EDAs follow this or a similar approach, yet treating orbitals as a link between quantum mechanics and chemical insight produces unintended consequences. Initially, each orbital must be assigned to one fragment, which can be achieved by using localized orbitals [e.g., natural bonding (NBO³⁷) in NEDA^38,39 or absolutely localized⁴⁰ in ALMO-EDA^7,41 (which is equivalent to the block-localized wavefunction²⁸)] or fragment molecular orbitals (FMOs⁴² in PIEDA¹⁹). Localized MOs have the additional advantage of permitting the study of through-bond phenomena, such as conjugation and hyperconjugation.⁴³ On the other hand, these localized orbitals also result in an increased sensitivity to the choice of basis set.^44,45

From a philosophical perspective, it could be argued that orbitals are not observables. Rather, they are an auxiliary concept that is potentially helpful for constructing a wavefunction. As a result, some have claimed that extracting information from these mathematical constructs, as opposed to taking information directly from the actual physical quantities, may lead to false conclusions.^46,47 In particular, it could be imagined that the same total binding energies might be obtained in ALMO-EDA (using a DFT wavefunction) and from a coupled cluster wavefunction in which the corresponding orbitals are entirely different. The resulting decomposition, $Δ E = E F R Z + E P O L + E C T$⁠, would contain very different contributions from the frozen density term (i.e., the sum of the electrostatic and the exchange interaction, E_FRZ), the polarization energy, E_POL, and the charge transfer energy, E_CT, based on the chosen wavefunction. In practice, however, the divergence has only been pointed out when diffuse basis functions were employed. Relying on orbitals also causes results to be particularly susceptible to the quality of the wavefunction, which is further complicated by the fact that assigning orbitals to fragments precludes the use of extrapolation techniques to obtain the complete basis set (CBS) limit.

For chemists, interpretation based on orbitals is a very intuitive process. Orbital-based definition of, for example, the charge-transfer term is quite natural and aligns well with experimentalists’ ideas. As shown in Fig. 4, the value of the charge-transfer term depends on both the EDA schemes used (with NBO-EDA results deviating significantly from the other approaches) and the wavefunction. Specifically, the amount of exact exchange in the density functional dramatically impacts the magnitude of the CT term. Despite this, each EDA scheme (including those that are unable to capture the binding nature of the interaction) identifies CT as significantly stabilizing the system.

An alternative to the orbital-based approaches is the density-based EDA (DEDA) proposed by Wu et al.,⁴⁹ which relies on a constrained relaxation of the electron density. By replacing the orbitals with the electron density (an observable that seems particularly natural within the framework of density functional theory), the sensitivity of basis set choice on the CT term is diminished and the E_FRZ term is introduced in a variational way.⁵⁰ However, separating the CT term in DEDA employs constrained DFT⁵¹ that places an arbitrary constraint on the charge distribution that potentially influences the outcome.⁵² Thus, employing the density, as opposed to orbitals, moves the arbitrariness from one location to another, rather than eliminating it entirely.

The intuitive definition of charge transfer (CT) is amongst the strongest points of the variational EDAs. In contrast, London dispersion is ill defined within this framework. Variational EDAs either completely neglect this interaction (as in KM EDA) or it is added on top of the decomposition in an ad hoc manner.

The neglect of dispersion is, in part, a relic from the time in quantum chemistry when the only available wavefunctions were uncorrelated, i.e., they employed a mean field approximation. Moreover, the original intent of EDAs was to describe covalent bonding. Because London dispersion interactions are relatively weak, long-range interactions were believed to only minimally contribute to the total energy of a covalent bond. However, for neutral complexes possessing symmetric charge distributions, dispersion is the only attractive interaction present. As shown by one of us, those weak interactions even play a dominant role in ground state charge transfer complexes.⁵⁰ Since the contribution from London dispersion forces is strictly attractive in the ground state, the magnitude of the interaction grows with the system size and can even increase to the level of covalent bonds in magnitude.⁵³ 

In wavefunction based EDAs, the terms “dispersion” and “correlation” are often used interchangeably, while in a majority of DFT based EDAs “dispersion” is defined directly from the chosen correction that is applied to the chosen density functional. Proceeding in this manner is problematic primarily because most atom pairwise dispersion corrections are fitted to a specific functional. Thus, the “correction term” may not only comprise the actual dispersion contribution but also include a portion of the correlation energy that was not captured by the original functional.⁵⁴ In fact, the energetic contribution to the “dispersion correction” is better used to gauge the inadequacy of the original functional than to measure any physical interaction (see Box 2). This is because such corrections were originally developed to cure the inability of DFT to describe long-range interactions that, in regions of overlapping electron density, are generally treated directly by the exchange-correlation functional. For functionals that explicitly account for medium-range correlation effects (e.g., the Minnesota density functionals⁵⁵), the situation becomes more opaque. The dispersion contributions to the benzene and water dimers presented in Figs. 5 and 6 are illustrative: for the benzene dimer, the dispersion component is correctly transferred to the correlation term (CORR) using GKS-EDA, but confusingly for the water dimer the dispersion contribution is redistributed amongst each of the other components. More recently, Head-Gordon and co-workers^36,56 investigated this problem using ALMO-EDA. They defined the dispersion term as the difference between the interaction energies from the exchange-correlation component of the functional that comprises dispersion and from exchange-correlation component of a “dispersion-free” (i.e., exact-exchange of rev-PBE⁵⁷) functional. While this certainly assists in providing a clearer definition of the van der Waals term, it does not completely resolve the problem.

Variational EDA schemes are generally cheap tools that provide sound chemical pictures of different interactions. Because of their efficiency (in particular, the FMO approaches^19,59 and the linear-scaling EDA of Phipps et al.⁶⁰) and their ability to treat many-fragment systems,⁶¹ they are instrumental in the continuing effort to better understand liquid water and ice,^62–64 cooperative effects,^25,65–68 and other interactions in solution.^69,70 Moreover, they treat both inter- and intramolecular^18,19,43,71 (including through-bond) interactions on equal footing (including through bond) and, thus, capably describe phenomena that are of a mixed or non-covalent nature, e.g., dative bonds,^27,72,73 as well as chemisorption and physisorption processes ^221,222 (through a “periodic EDA” by Tonner et al.²²³). Owing to the numerous levels of accuracy possible, including multireference variants,^41,74 it is also possible to treat strongly correlated systems⁷⁵ and molecules having elongated bonds. Owing to an unambiguous definition of CT, EDAs are often used to identify and quantify charge transfer interactions,^27,22,48,76 frequently in conjunction with ALMO-CTA,^77,78 which is complementary to ALMO-EDA.

While variational EDAs ably describe many interesting facets of inter- and intramolecular interactions, they are not without their own shortcomings. For example, the principle drawback is their poor performance for describing dispersion interactions. Addition of a “–D” style correction or identification using correlation makes EDAs suitable for describing van der Waals complexes by masking, rather than solving, the problem. EDAs are able to describe non-covalent interactions, but not sufficiently reliable to describe dispersion driven interactions or provide models for further method development. Moreover, this class of approaches is highly sensitive to the wavefunction quality. Thus, the specific values of different energy components can only be interpreted and compared within a single EDA scheme. Even within a single scheme, the values obtained from different density functionals are not comparable with one another. Finally, the supermolecular nature of these computations introduces a significant dependence on the basis set. The basis set superposition error⁷⁹ (BSSE) also becomes problematic. Conversely, that same supermolecularity prevents the theory from breaking down for strong interactions, such as doubly hydrogen bonded systems⁸⁰ and, recently, permitted the formulation of “adiabatic EDA,” which directly investigates the influence of different interaction types on structural and other molecular properties.⁸¹ 

An alternative to studying a system’s wavefunction is to examine the other ingredient of the Schrödinger equation: the energy operator. The basic concept of perturbational EDAs is to “correct” the Hamiltonian of the non-interacting system to include the desired interaction. Each component of the correction is then given a physical interpretation. Thus, the energy operator becomes a central figure for this class of approaches. Because perturbation theory assumes that the magnitude of the interaction is small, these approaches are restricted to non-covalent interactions. By far the most successful perturbation-based EDA is symmetry-adapted perturbation theory (SAPT),⁶ which is inspired by the Rayleigh-Schrödinger method where the Hamiltonian of dimer AB is given as

where $H ^ A$ and $H ^ B$ are the unperturbed Hamiltonians of systems A and B and V is the interaction operator. In its simplest formulation, SAPT0⁸² [or SAPT(HF)], the Hamiltonians of the operators are simply Fock matrices and the zeroth-order wavefunction is a product of the Slater determinants representing the monomers. In order to recover the exchange interaction, dimer wavefunctions of an order higher than zero are antisymmetrized. The sophistication level can be increased by including intramonomer correlation through, for example, Møller-Plesset type fluctuation operators in the monomer Hamiltonians or by using Kohn-Sham⁸³ or coupled-cluster wavefunctions.⁸⁴ Thus, SAPT can be employed both as an energy decomposition analysis scheme and as a method providing accurate interaction energies.

The energy components found in SAPT include the usual terms: electrostatic, exchange, and dispersion. Rather than charge-transfer and polarization, SAPT uses a component that comprises both of those terms, namely, the induction energy, which arises from the singly excited monomer wavefunctions. Because of symmetry, additional terms appear in the higher-order expansions of the interaction energy, such as induction-exchange and dispersion-exchange.

The treatment of the exchange interaction distinguishes SAPT from polarization theory,⁸⁵ where exchange is neglected but also from modern perturbation-based approaches that either avoid this problem by decomposing only the correlation energy for a fully antisymmetric wavefunction of the dimer (see Refs. 86 and 87) or taking advantage of localized molecular orbitals (mainly biorthogonal approaches by Mayer^88–90).

The form of the Hamiltonian, specifically the perturbation series of the interaction operator, dictates energy decomposition in SAPT. Exchange-type components are extracted in each perturbation order as the difference between the interaction energies obtained with antisymmetrized wavefunctions and those obtained from the Rayleigh-Schrödinger theory. The remainder of the first-order energy comes from polarization (or simply electrostatics) that corresponds to the Coulomb interaction of the non-deformed monomers. For the second-order SAPT, the induction (accompanied by the repulsive induction-exchange component, identified with a correction to Pauli repulsion, and interpreted as a polarization of one monomer in the static electric field produced by the other monomer) and dispersion (accompanied by a repulsive dispersion-exchange and interpreted as the correlated instantaneous multipole moments of monomers) components can be distinguished from one another,

The three main components, E_ELST, E_IND, and E_DISP, are also defined within the Rayleigh–Schrödinger perturbation theory and each of them can be expressed through monomer properties: the electrostatic term relates to electron densities, the dispersion components to dynamic multipole polarizabilities, and induction (approximately) to multipole moments and polarizabilities.⁶ However, higher-order terms, such as induction-dispersion and third-order dispersion, lose any conventional physical meaning. In contrast to the variational EDAs, the SAPT components arise naturally rather than being extracted from the energy. The components also sum to the total interaction energy, which means they can be systematically improved through the expansion of the interaction operator and the intramonomer correlation operator.⁸⁰ 

The greatest asset of SAPT is certainly its ability to quantify dispersion at all ranges (see Box 2), including regions of overlapping electron density. Generally, SAPT weights dispersion more heavily than other approaches (for example, see the water and benzene dimer examples presented in Figs. 5 and 6 and Refs. 92 and 93). In contrast to dispersion, induction, which replaces the polarization and charge-transfer terms, is not as chemically intuitive as charge-transfer. Some argue that this is actually advantageous, since polarization and CT cannot be distinguished in the CBS limit.^2,94 For practical purposes, one can simply consider the entire second-order perturbation energy bar the dispersion (E_CTPOL)—like in the ClF-NF3 complex (see Fig. 4), for which the E_CTPOL values are quantitatively different from those obtained with the variational EDAs, but they do follow the same trends. Alternatively, the CT can be extracted from a SAPT computation using an approach akin to ALMO-CTA,⁷⁷ where the CT is the difference between the induction energies computed in the dimer-centered basis set and monomer basis set.⁹⁵ This scheme is strongly basis set dependent to the point that the CT term can even become positive (i.e., the interaction is viewed as destabilizing).⁹⁶ This dependence, however, can be overcome through the regularized SAPT approach⁹⁷ or constrained DFT.⁹⁶ Thus, values of charge-transfer can be computed within SAPT, but this quantity is not consistent with the theory.

For standard formulations of SAPT, dimers and trimers^106–108 make perfect model systems which can be easily described. Of course, the full realm of non-covalent interactions is far vaster. As an example, today the importance of both cooperativity^109–111 and intramolecular dispersion^112–115 in various aspects of chemistry is becoming increasingly clear. To date, no SAPT variant which is able to treat systems consisting of more than three fragments has been proposed and other perturbation based many-fragment approaches are also scarce (the efficient, but approximate XSAPT,¹¹⁶ the Chemical Hamiltonian approach of Halász et al.¹¹⁷ and expensive, highly basis-set dependent supermolecular Møller-Plesset PT methods¹¹⁸ are noteworthy, but scarcely applied, examples). For standard SAPT, the treatment of many-fragment systems is both expensive and complicated, as it necessitates imposing the antisymmetry condition multiple times. Of course, the entire cluster of molecules could be treated as interacting dimers or trimers, at which point SAPT would again become applicable.^119,120

The analysis of intramolecular interactions is another field of applicability that remains challenging for perturbation theory. In these situations, obtaining the interaction operator that first requires extracting the non-interacting Hamiltonian and wavefunction is nontrivial considering that the electrons are indistinguishable and that it is necessary to preserve the covalent bonds. One option is utilizing localized orbitals and decomposing the correlation energy consisting of the interaction between two electron domains with the local MP2 framework.⁸⁶ Interestingly, a tool permitting visualization of the spatial contributions to the dispersion energy was recently established.¹²¹ Similar to the LMP2 approach but applied to the entire interaction energy is intra-SAPT,^122,123 for which the zeroth-order Hamiltonian is inspired by the chemical Hamiltonian of Mayer.⁸⁸ Unfortunately, this method lacks a CBS limit and involves somewhat complicated mathematics. An alternative to intra-SAPT is “ISAPT”¹²⁴ which derives its origins from two methods (A-SAPT¹²⁵ which associates contributions to the SAPT interaction energy with atoms and orbitals, and F-SAPT¹²⁶ which connects those contributions with functional groups within molecules) for visualizing interaction contributions. ISAPT uses a similar partitioning to A-SAPT and F-SAPT, along with HF-in-HF density matrix embedding¹²⁷ to construct the zeroth-order wavefunction. ISAPT is computationally efficient and easy to use but may suffer from the consequences of using a single-exchange approximation and spurious dipole-dipole interactions that originate from the fragment partitioning. Importantly, both ISAPT and intra-SAPT are in their infancy and their continued development would likely benefit from novel fragmentation approaches¹²⁸ along with new ways of partitioning the Hamiltonian.¹²⁹ 

While other approaches are still being developed, SAPT clearly remains the most prominent perturbation based EDA and popular method for achieving highly accurate interaction energies.^130–133 This method is a self-contained and consistent theory of molecular interactions built upon rigorously defined energy components and, unlike variational EDAs, has a non-empirically defined short-range dispersion energy. Thus, SAPT is an essential tool for obtaining benchmark values useful for designing and building new computational methods.^{13,12,58,134,135} Moreover, the approach can be systematically improved owing to the double-perturbation nature of its standard formulation. Currently, the coupled-cluster variant of SAPT has been shown to be very accurate for many systems⁸⁴ (and the much more computationally efficient SAPT(DFT) is not far behind),^83,136 although the treatment of multi-reference systems remains problematic.^137,138 Fortunately, in recent years an open-shell variant of SAPT has been developed,^139,140 which is an important prerequisite for a version capable of treating systems with a multi-reference character, including molecules reacting with one another and metallic systems.

An earlier bottleneck of SAPT was its computational cost, but efficient implementations^92,141–143 have seen SAPT applied to increasingly larger systems like fullerenes¹⁴⁴ and nanotubes.¹⁴⁵ Still, SAPT generally remains more costly than variational EDAs and visualization tools (vide infra). The work to improve efficiency has produced the aforementioned XSAPT^116,146 as well as linear-scaling SAPT,¹⁴⁷ but these approaches use modified or semiempirical dispersion expressions. Nevertheless, there are reasons to believe that the computational cost is becoming of diminishing importance.

Finally, as in all perturbation theories, SAPT may break down when dealing with particularly strong interactions. A frequent origin of the failure of SAPT is the single-exchange (or S²) approximation. In order to reduce the cost of solving the SAPT equations, the consequence of having many electrons exchanging places is ignored during the antisymmetrization. Generally, this approximation is very good, although it tends to break down for charged systems at close distances.^148–150 This failure can be overcome by scaling the exchange contribution⁸⁰ or by completely abandoning the S² approximation.¹⁵¹ 

Overall, SAPT’s lack of versatility remains perhaps its only true weakness. The “comfort zone” of the approach involves moderately sized neutral dimers and trimers. This lack of versatility contrasts virtually with all other aspects of SAPT, which are very strong. Although no unambiguous way to partition the interaction energy into components exists, SAPT fulfills all the criteria that make a decomposition scheme appealing: all the terms are physically interpretable, they reduce to simple models in the asymptotic limit and they add up to full interaction energy, which is exact in the infinite order of SAPT.

What is more, even when using a wavefunction of relatively poor quality, reasonable interaction energies and components can generally be obtained. For these reasons, SAPT is a much more “black-box” approach than the variational EDAs. Interestingly, it is also a unique tool that provides benchmark dispersion data for areas of overlapping densities, which is critically important for developing new computational methods.

Finally, we focus on a final class of approaches that reveals the presence and nature of non-covalent interactions in real space. This class of visualization tools, which are complementary to those that provide quantitative feedback, is particularly attractive for identifying specific molecular regions that are most affected by an interaction. For visualization purposes, the electronic density is a logical and well-suited variable. It is a quantum observable that yields the ground state energy¹⁵² and has a real space connotation. In fact, the lone drawback of using the electron density as a visualization tool is its lack of sensitivity, which makes capturing weak non-covalent interactions problematic. In particular, the electron density is less sensitive than the energy to van der Waals interactions,⁹⁹ even though noticeable modifications have been observed for systems with high polarizability density or low dimensionality—see, e.g., Ref. 153. Alternatively, one might also consider molecular electrostatic potential (MEP) maps¹⁵⁴ as a treatment for visually examining non-covalent interactions. Akin to directly observing the electron density, while MEPs ably describe strong electrostatic interactions such as those present in halogen and hydrogen bonds,^155–158 they lack the sensitivity (or can even be misleading) for interpreting weaker interactions (see, e.g., Ref. 159). Another possibility for visualizing molecular interactions is the Electron Localization Function (ELF),^160,161 which is a scalar field associated with same-spin pair densities of electrons that is a popular descriptor for accessing both the bonding and atomic shell structure. Unfortunately, while extremely useful for characterizing covalent bonding, ELF only recognizes non-covalent interactions arising from strong electrostatic forces, mainly those occurring in hydrogen bonds,^162,163 halogen bonds,^164,165 and pnicogen bonds.¹⁶⁶ 

Given that each of the above tools fails to properly describe weak non-covalent interactions, it becomes necessary to search for an alternative approach elsewhere. If one concedes that the electron density does contain all the relevant information about a system, then an attractive possibility is look at the density topology, i.e., the density’s gradient field and Hessian. This is the idea behind Bader’s Quantum Theory of Atoms in Molecules (QTAIM),^3,167 which exploits the topology of the density to redefine atoms and bonds in molecules in a rigorous mathematical framework, consistent with the principles of quantum physics. Two other tools also capably detect and permit the visualization of weak interactions that are partially inspired by QTAIM: the Non-Covalent Interaction index⁵ and the Density Overlap Region Indicator.⁴ Each of these three approaches will be discussed in more detail below.

QTAIM defines atoms as regions of space with the maximum of density located on the nucleus and bounded by a surface orthogonal to the density gradient. The “bond path,” a line of locally maximum density connecting two maxima placed on the nuclei (“nuclear critical points”), replaces the notion of a bond. Unconventionally, these bond paths are often not straight (see Fig. 7) and they may appear for both covalent and non-covalent patterns.^168–171 Apart from visualizing bond paths and bond critical points (BCPs, maxima of the density along bond paths), QTAIM facilitates the analysis and interpretation of interactions through parameters such as the Laplacian of the electron density and its eigenvalues, which are computed at the BCPs, or the curvature of the bond paths.

QTAIM translates the traditional chemistry language (see Table I) into mathematical constructs. As a result, there is no fundamental discrimination, for instance, between a covalent bond and a non-covalent interaction, which leads to non-intuitive depictions of van der Waals complexes (see, e.g., the thiophene dimer in Fig. 7). At the extreme, are situations such as that presented in Fig. 8, which depicts the large number of London dispersion interactions between a noble gas and a fullerene (He@C₆₀). This system contains 60 bond paths, which is, perhaps, physically realistic but contrary to normal chemical intuition. The bond paths are better thought of as representing interacting atoms rather than any actual chemical bond.¹⁷² QTAIM identifies all atomic interactions as “local pairings of the densities of opposite spin electrons”¹⁷³ which, in fact, covers all non-covalent interactions. However, when the density is computed from DFT, which only roughly accounts for dispersion, it becomes unclear whether the van der Waals interactions remain correctly depicted by QTAIM.

QTAIM is essentially known as a visualization tool (thanks to its availability in various software) but an energy decomposition analysis tool has also been derived from its framework. The Interacting Quantum Atom (IQA) approach^174,175 partitions the energy components into intraatomic and interatomic terms by associating the terms of the one- and two-particle reduced density matrices (RDMs) with atomic basins. Moreover, the two-particle reduced density matrices are partitioned into Coulomb, exchange, and correlation terms, the latter of which is introduced through the use of a Kohn-Sham,^175,176 coupled-cluster,^177,178 CASSCF, or MP2 wavefunction.¹⁷⁴ Even though the energy partitioning in IQA is less clear-cut than SAPT (the physical interpretation of the correlation term is not as straightforward as the dispersion or the charge-transfer components), IQA has the advantage of relying upon a clear definition of atoms and fragments (as it is built on the QTAIM framework). This becomes particularly important for intramolecular interactions when orbital-based atom definitions become blurred by the existence of covalent bonds. On the other hand, this real-space partitioning also has a downside—the computational inefficiency associated with performing numerous integrations on a 3D grid, although new faster algorithms are continuously being developed (see, e.g., Refs. 179 and 180).

The previously discussed QTAIM is a mathematically rigorous, intrinsically consistent theory that facilitates the visualization of bonds and other interactions. However, it is easy to envision cases where a simpler and more intuitive tool would be desirable, something akin to ELF but that allows the visualization of van der Waals interactions. A noteworthy method that accomplishes exactly this is the Non-Covalent Interaction (NCI) index.^5,182–184 NCI adapts ideas from both DFT and QTAIM. The actual NCI index is the reduced gradient s of the density,

which, within the DFT framework, describes to what degree the system differs from a homogeneous electron gas. An isosurface of s determines the spatial area of the interaction, while the sign of the second eigenvalue of the density Hessian determines whether the interaction is attractive or repulsive. In the case of a thiophene dimer (see Fig. 7), NCI identifies two types of interacting regions: a steric repulsion (drawn in red) in the ring centers and a weakly attractive (in green) π-stacking between rings.

Interestingly, NCI is capable of revealing weak interactions, including pure London dispersion, even when the density is derived from a method that poorly describes van der Waals forces.¹⁸⁵ Usually,^5,182 although not always,¹⁸⁶ NCI correctly identifies these interactions even when promolecular densities (i.e., sums of atomic densities) are used. In turn, that allows the treatment of very large systems such as DNA fragments^187,188 and proteins.¹⁸⁹ The choice of the computational method, however, does have a quantitative impact on the NCI index,¹⁸⁵ but many users of NCI are more interested in qualitative trends than quantitative values.

NCI is able to handle also the covalent interactions,^5,182,190 but to obtain a clear and intuitive picture of both strong and weak interactions at the same time the combination of NCI with ELF is needed.¹⁹¹ This is a significant nuisance as a huge advantage of the original NCI is its simple, nearly “black-box” use.

Alternatively, for the simultaneous visualization of both covalent and non-covalent interactions in real space, one can use the Density Overlap Region Indicator (DORI),⁴ which is also a scalar field based on the density and the density gradient. DORI descends from the aptly named Single Exponential Decay Detector (SEDD),¹⁹² which finds regions corresponding to localized electrons, based on the observation that the electron density in those regions should be described by exponential functions (as they correspond to single orbitals). DORI takes this philosophy further and associates all the density overlap regions (i.e., bonds, molecular interactions, and steric clashes) with their deviations from exponentiality. The detected interactions are then identified as attractive or repulsive using, as in NCI, the sign of the second eigenvalue of the density Hessian. In non-covalent cases, DORI produces pictures similar to NCI (compare Figs. 7(b) and 7(c)) while also capturing the bonding patterns of ionic and covalent bonds without any special adjustments.

An appealing albeit not extensively explored alternative to considering the derivatives of the density is examining the density response. A suitable object is the density-density response function χ(r,r’,ω) that connects the external perturbation potential V(r’,ω) with a system’s linear density response

Even though $ρ r e s p ( 𝐫 , ω )$ is a multidimensional object, it is much less complex than the wavefunction. More importantly, it describes the (local) susceptibility of the system to an oscillating electrical field, which, in turn, obviously links to charge fluctuations and the physics of the dispersion interaction. Preliminary efforts to exploit this link have been applied to molecular crystals^193,194 and weakly interacting complexes.¹⁹⁵ 

As mentioned earlier, QTAIM is a widely applicable and rigorous theory, but the delivered insight is sometimes counterintuitive or cryptic. QTAIM and IQA have helped decipher the nature of halogen^{186,196–198} and hydrogen^{168,199–202} bonds, but discussion on their use within the context of London dispersion interactions is still ongoing.²⁰³ Even if QTAIM can be used for visualizing van der Waals interactions, interpreting the results is by no means “black-box.” Additionally, QTAIM analysis is generally computationally expensive, thereby preventing its application to large systems. IQA appears promising and may serve as an appealing alternative to traditional variational EDAs in the decomposition of intramolecular interactions. The relevance of IQA for treating weak interactions will expand once a better understanding of the electronic correlation terms is achieved, for which progress is currently being made.^177,178,204

Visualization tools like NCI and DORI have a different objective than QTAIM and IQA: they were designed to be pragmatic, easy to use, and computationally efficient. As a result, they capably identify the spatial extent and the underlying nature of interactions present in molecular systems and complexes. The relevance of NCI for studying large systems has been demonstrated on DNA, proteins, as well as nanotubes,^205,206 polymers,²⁰⁷ and crystals under high pressure.²⁰⁸ Moreover, this tool has also been used to explore reactivity,^{114,184,209,210} dynamics,¹⁹¹ and solvation²¹¹ effects. Likewise, DORI has been used to explain the mechanisms of alkylperoxyl radical termination,²¹² to investigate the conformational behavior of peptides in different media,^213,214 to design new materials,²¹⁵ and has also served as a component of new density functionals.^216,217

As both NCI and DORI rely only upon the electronic density and its derivatives, they can be applied to densities obtained experimentally. Similarly, they can analyze excited states,^218,219 whenever an excited state density is available. Provided that one is interested in qualitative pictures, the range of applicability of scalar fields is therefore quite immense. The easy interpretation and low sensitivity to the underlying computational method make these visualization tools particularly attractive to experimentalists and to others not wishing to invest considerable computational resources.

This perspective overviewed the tools used to translate the highly mathematized quantum chemical description of non-covalent interactions into more intuitive concepts and graphical depictions. They function by extracting information from well-defined objects of the quantum physics world: the wavefunction, the energy operator (i.e., Hamiltonian), and the electron density. To a large degree, the strengths and weaknesses of an individual tool strongly depend upon which physical object is used in their construction. For example, wavefunction-based approaches are unlikely to break down for strong interactions, but, in contrast to Hamiltonian-based schemes (e.g., SAPT), they may lack the sensitivity to capture weak dispersion interactions. Density-based tools (i.e., QTAIM, DORI, and NCI), on the other hand, are often employed for visualization purposes. Of course, the wavefunction, the Hamiltonian, and the electron density are all interrelated, so categorizing the tools on this basis is somewhat artificial. Nevertheless, both potential applications and insights obtained from a specific class of approaches should be analyzed by keeping their limited purpose in mind. These inherent differences create a modus operandi for each tool, where some phenomena and notions would have multiple interpretations while others do not have a true counterpart. Illustrative examples of such “idioms” are SAPT’s dispersion-exchange and induction-exchange energy components and the representation of vdW interactions in QTAIM. The complete removal of the ambiguity built into the tools is unlikely, but would be highly desirable. Toward achieving that objective, two principle developmental directions emerge: unifying existing approaches to avoid ill-defined terms^{122–124,128} and expanding the field capabilities of existing approaches.^{56,139,124,125,179,219,220} The first avenue has already lead to a more stable definition of charge transfer^77,96 or London dispersion terms.³⁶ Alternatively, SAPT has recently been expanded to treat intramolecular interactions^123,124,128 and open-shell systems.¹³⁹ The latter achievement could be a major step toward the development of a multi-reference (i.e., able to treat bond-breaking phenomena) or even an excited-state variant of SAPT. Additional efforts for reducing the computational cost are also on-going.^142,147 A rudimentary judgment of recent developments following each of these strategies leads to the conclusion that this is an exciting time in the domain of computational tools for the interaction analysis.

Funding from EPFL and from Swiss NSF Grant No. “200021_156001” is gratefully acknowledged. We thank Mr. Laurent Vannay, Dr. Ganna Gryn’ova, and Dr. Matthew Wodrich for insightful discussions.