Physical mechanisms of intermolecular interactions from symmetry-adapted perturbation theory

Abstract

Symmetry-adapted perturbation theory (SAPT) is a method for computational studies of noncovalent interactions between molecules. This method will be discussed here from the perspective of establishing the paradigm for understanding mechanisms of intermolecular interactions. SAPT interaction energies are obtained as sums of several contributions. Each contribution possesses a clear physical interpretation as it results from some specific physical process. It also exhibits a specific dependence on the intermolecular separation R. The four major contributions are the electrostatic, induction, dispersion, and exchange energies, each due to a different mechanism, valid at any R. In addition, at large R, SAPT interaction energies are seamlessly connected with the corresponding terms in the asymptotic multipole expansion of interaction energy in inverse powers of R. Since such expansion explicitly depends on monomers’ multipole moments and polarizabilities, this connection provides additional insights by rigorously relating interaction energies to monomers’ properties.

Appendices

Appendix: Third-order induction energy

In contrast to the second-order case, separation of the induction energy in higher orders is not straightforward. A definition of the infinite-order induction energy was proposed in ref. [53] (see also a more extensive discussion of this approach in ref. [107]). This definition leads to the formula for $E^{(3)}_{\text {ind}}$ given by Eq. (60) in ref. [107] (or equivalently by Eq. (20) in ref. [53] where the reduced-resolvent notation is used). In both cases, $E^{(3)}_{\text {ind}}$ is expressed via the zeroth- and first-order functions only. To get an alternative physical interpretation of this quantity, we derive here an expression for $E^{(3)}_{\text {ind}}$ involving only the second-order induction functions. The infinite-order induction energy can be defined [53] as the minimum of the expectation value of H₀ + V with the trial function of the form $\tilde {\Psi }^{A}_{\text {ind}} \tilde {\Psi }^{B}_{\text {ind}}$

$$J_{\text{ind}}\left[\tilde{\Psi}^{A}_{\text{ind}},\tilde{\Psi}^{B}_{\text{ind}}\right] \equiv \frac{\langle \tilde{\Psi}^{A}_{\text{ind}} \tilde{\Psi}^{B}_{\text{ind}}| H_{0} +V| \tilde{\Psi}^{A}_{\text{ind}} \tilde{\Psi}^{B}_{\text{ind}} \rangle} {\langle \tilde{\Psi}^{A}_{\text{ind}} \tilde{\Psi}^{B}_{\text{ind}}| \tilde{\Psi}^{A}_{\text{ind}} \tilde{\Psi}^{B}_{\text{ind}} \rangle} .$$

(A.1)

A simple way to derive the equations for the induction functions is to first assume that the exact induction wave function for monomer B is known. Equation (A.1) can then be written as

$$J_{\text{ind}}\left[\tilde{\Psi}^{A}_{\text{ind}},{\Psi}^{B}_{\text{ind}}\right] = \frac{\langle \tilde{\Psi}^{A}_{\text{ind}} | H_{A} + \bar{\Omega}_{B}| \tilde{\Psi}^{A}_{\text{ind}} \rangle} {\langle \tilde{\Psi}^{A}_{\text{ind}}| \tilde{\Psi}^{A}_{\text{ind}} \rangle} ,$$

(A.2)

where

$$\begin{aligned} \bar{\Omega}_{B} = \frac{\langle {\Psi}^{B}_{\text{ind}}| V {\Psi}^{B}_{\text{ind}} \rangle + \langle {\Psi}^{B}_{\text{ind}}|H_{B} |{\Psi}^{B}_{\text{ind}} \rangle} {\langle {\Psi}^{B}_{\text{ind}}|{\Psi}^{B}_{\text{ind}} \rangle} . \end{aligned}$$

(A.3)

The standard application of the variational principle leads then to the following equation for the induction wave function

$$\begin{array}{@{}rcl@{}} \left[ H_{A} + \bar{\Omega}_{B} \right] {\Psi}^{A}_{\text{ind}} = {\mathcal E} {\Psi}^{A}_{\text{ind}} , \end{array}$$

(A.4)

where ${\mathcal E} = J_{\text {ind}}[{\Psi }^{A}_{\text {ind}}, {\Psi }^{B}_{\text {ind}}]$. The analogous equation for monomer B is

$$\begin{array}{@{}rcl@{}} \left[ H_{B} + \bar{\Omega}_{A} \right] {\Psi}^{B}_{\text{ind}} = {\mathcal E} {\Psi}^{B}_{\text{ind}} . \end{array}$$

(A.5)

If all quantities in Eqs. (A.4) and (A.5) are expanded in powers of V, one gets (using a self-explanatory short-hand notation)

$$\begin{array}{@{}rcl@{}} {\mathcal E} = {\mathcal E}^{0} + {\mathcal E}^{1} + {\mathcal E}^{2} + {\mathcal E}^{3} + {\dots} = {E_{0}^{A}} + {E_{0}^{B}} + E^{(1)}_{\text{elst}} + E^{(2)}_{\text{ind}} + E^{(3)}_{\text{ind}} + {\dots} , \end{array}$$

(A.6)

$$\begin{array}{@{}rcl@{}} |A\rangle = |A^{0}\rangle + |A^{1}\rangle + |A^{2}\rangle + {\dots} , \end{array}$$

(A.7)

and similarly for monomer B, where all the terms in the sums except for the last ones have been defined before. To derive formulas for the latter quantities, we start from the following expansion of $\bar {\Omega }_{B}$:

$$\begin{array}{@{}rcl@{}} \bar{\Omega}_{B} = {{\Omega}^{0}_{B}} + {{\Omega}^{1}_{B}} + {{\Omega}^{2}_{B}} + {{\Omega}^{3}_{B}} + {\dots} , \end{array}$$

(A.8)

where

$$\begin{array}{@{}rcl@{}} {{\Omega}^{0}_{B}} &=& {E_{0}^{B}}; \quad {{\Omega}^{1}_{B}} = \langle B^{0}|V|B^{0} \rangle = {\Omega}_{B} + \langle B^{0}|V^{A}|B^{0} \rangle + V_{0}; \end{array}$$

(A.9)

$$\begin{array}{@{}rcl@{}} {{\Omega}^{2}_{B}} &= &2\langle B^{0}|V|B^{1} \rangle - \langle B^{0}|{\Omega}_{A}|B^{1} \rangle \end{array}$$

(A.10)

and

$$\begin{aligned} {{\Omega}^{3}_{B}} &= \langle B^{1}|V|B^{1} \rangle + 2 \langle B^{0}|V|B^{2} \rangle\\& - 2 \langle B^{0}|{\Omega}_{A}|B^{2} \rangle - \langle B^{0}|V|B^{0} \rangle\langle B^{1}|B^{1} \rangle . \end{aligned}$$

(A.11)

The second-order equation defining $|B^{2} \rangle$ is

$$\begin{aligned} (H_{B} - {E_{0}^{B}}) |B^{2}\rangle& = \left(\langle B^{0}|{\Omega}_{A}|B^{0} \rangle-{\Omega}_{A}\right) |B^{1}\rangle\\ &+ \left(\langle B^{0}|{\Omega}_{A}|B^{1} \rangle -2 \langle A^{0}|V|A^{1} \rangle \right) |B^{0}\rangle , \end{aligned}$$

(A.12)

with an analogous equation for $|A^{1} \rangle$. Finally, the third-order equation is

$$\begin{array}{@{}rcl@{}} (H_{A} - {E_{0}^{A}}) \big|A^{3}\rangle = \sum\limits_{i=1}^{3} \left({\mathcal E}^{i} - {{\Omega}_{B}^{i}} \right) |A^{3-i}\rangle , \end{array}$$

(A.13)

resulting in the following expression for the third-order induction energy

$$\begin{aligned} {\mathcal E}^{3} &= \langle A^{0}|{\Omega}_{B}|A^{2} \rangle + 2 \langle A^{0} B^{0} | V | A^{1} B^{1} \rangle \\&+ \langle B^{1}|{\Omega}_{A}|B^{1} \rangle - \langle B^{0}|{\Omega}_{A}|B^{0} \rangle\langle B^{1}|B^{1} \rangle . \end{aligned}$$

(A.14)

Since

$$\begin{aligned}\langle B^{1}|{\Omega}_{A}|B^{1} \rangle & = \langle B^{0}|{\Omega}_{A}|B^{2} \rangle + \langle B^{0}|{\Omega}_{A}|B^{0} \rangle\langle B^{1}|B^{1} \rangle\\ & - 2 \langle A^{0} B^{0} | V | A^{1} B^{1} \rangle , \end{aligned}$$

(A.15)

we finally obtain

$$E^{(3)}_{\text{ind}} \equiv {\mathcal E}^{3} = \langle A^{0}|{\Omega}_{B}|A^{2} \rangle + \langle B^{0}|{\Omega}_{A}|B^{2} \rangle .$$

(A.16)

Returning to the explicit notation, Eq. (A.16) becomes

$$E^{(3)}_{\text{ind}} = \langle {{\Psi}_{0}^{A}} |{\Omega}_{B}| {\Psi}_{\text{ind}}^{(2)A} \rangle + \langle {{\Psi}_{0}^{B}} |{\Omega}_{A}| {\Psi}_{\text{ind}}^{(2)B} \rangle .$$

(A.17)

To find an explicit expression for ${\Psi }_{\text {ind}}^{(2)A}$, let us first write the equivalent of Eq. (A.12) for monomer A:

$$\begin{aligned} (H_{A} - {E_{0}^{A}}) |A^{2}\rangle& = \left(\langle A^{0}|{\Omega}_{B}|A^{0} \rangle-{\Omega}_{B}\right) |A^{1}\rangle\\ & + \left(\langle A^{0}|{\Omega}_{B}|A^{1} \rangle -2 \langle B^{0}|V|B^{1} \rangle \right) |A^{0}\rangle , \end{aligned}$$

(A.18)

which leads to the following formula for the second-order induction function of monomer A expressed as a spectral sum

$$\begin{aligned} {\Psi}_{\text{ind}}^{(2)A}& = \sum\limits_{k\ne 0} \frac{\langle {{\Psi}_{k}^{A}} |\langle {{\Psi}_{0}^{A}} | {\Omega}_{B} {{\Psi}_{0}^{A}} \rangle -{\Omega}_{B}| {\Psi}_{\text{ind}}^{(1)A} \rangle} {{E_{k}^{A}} - {E_{0}^{A}}} {{\Psi}_{k}^{A}}\\ &- 2 \sum\limits_{k\ne 0} \frac{ \langle {{\Psi}_{k}^{A}} {{\Psi}_{0}^{B}} | V_{\text{ee}} | {{\Psi}_{0}^{A}} {\Psi}_{\text{ind}}^{(1)B} \rangle } {{E_{k}^{A}} - {E_{0}^{A}}} {{\Psi}_{k}^{A}} . \end{aligned}$$

(A.19)

In practice, $E^{(3)}_{\text {ind}}$ is much easier to compute using the expression depending only on the zeroth- and first-order functions, see Eq. (A.22) below, but the use of the second-order function allows a transparent interpretation of the third-order induction energy. Of course, the importance of this function stems also from the fact that it can be used to compute the fourth- and fifth-order induction energy contributions.

Expression (A.17) can be written in terms of densities as

$$\begin{array}{@{}rcl@{}} E^{(3)}_{\text{ind}} &= \frac{1}{2}\int \rho^{(02)A}_{\text{ind}}({\boldsymbol r}) \omega^{B}({\boldsymbol r}) d^{3}\boldsymbol{r} + \frac{1}{2}\int \rho^{(02)B}_{\text{ind}}({\boldsymbol r}) \omega^{A}({\boldsymbol r}) d^{3}\boldsymbol{r} \end{array}$$

(A.20)

where $\rho ^{(02)A}_{\text {ind}}({\boldsymbol r})$, a component of the total second-order induction density $\rho ^{(2)A}_{\text {ind}}({\boldsymbol r})$, is defined by a procedure similar to that leading to Eq. (21), but with ${\Psi }_{\text {ind}}^{(2)A}$ added to Ψ^A

$$\begin{aligned}\rho^{(2)A}_{\text{ind}}({\boldsymbol r}) &= 2 N_{A} \sum\limits_{s,s_{2},\dots,s_{N}}\int {{\Psi}_{0}^{A}}(\boldsymbol{x},\boldsymbol{x}_{2}, \dots, \boldsymbol{x}_{N_{A}}) {\Psi}_{\text{ind}}^{(2)A}(\boldsymbol{x},\boldsymbol{x}_{2}, \dots, \boldsymbol{x}_{N_{A}}) d^{3}\boldsymbol{r}_{2} {\dots} d^{3}\boldsymbol{r}_{N_{A}} \\ &\quad +N_{A} \sum\limits_{s,s_{2},\dots,s_{N}}\int {\Psi}_{\text{ind}}^{(1)A}(\boldsymbol{x},\boldsymbol{x}_{2}, \dots, \boldsymbol{x}_{N_{A}}) {\Psi}_{\text{ind}}^{(1)A}(\boldsymbol{x},\boldsymbol{x}_{2}, \dots, \boldsymbol{x}_{N_{A}}) d^{3}\boldsymbol{r}_{2} {\dots} d^{3}\boldsymbol{r}_{N_{A}} \\ & = \rho^{(02)A}_{\text{ind}}({\boldsymbol r}) + \rho^{(11)A}_{\text{ind}}({\boldsymbol r}), \end{aligned}$$

(A.21)

with the second-order induction function given by Eq. (A.19) and the two terms in the last expression corresponding to the consecutive terms in the preceding expression. Thus, the second-order induction densities have two components: one originating from the second-order induction function and one originating from the product of first-order induction functions, and only the former term enters the expression for $E^{(3)}_{\text {ind}}$. Equation (A.20) shows that the third-order induction interaction can be interpreted as the interaction of the second-order induction density components resulting from ${\Psi }_{\text {ind}}^{(2)X}$ with the unperturbed electrostatic potentials. Note again the factors of 1/2 multiplying the Coulomb interactions of $\rho ^{(02)A}_{\text {ind}}({\boldsymbol r})$ with ω^B(r) in the first term, and analogously in the second.

Using Eq. (A.15), formula (A.16), expressed in terms of the unperturbed functions and of the second-order induction functions, can be transformed into formula (60) from ref. [107] expressed in terms of the unperturbed functions and of the first-order induction functions only:

$$\begin{aligned} E^{(3)}_{\text{ind}} &= \langle A^{1}|{\Omega}_{B}|A^{1} \rangle + \langle B^{1}|{\Omega}_{A}|B^{1} \rangle - \langle A^{0}|{\Omega}_{B}|A^{0} \rangle\langle A^{1}|A^{1} \rangle\\ &- \langle B^{0}|{\Omega}_{A}|B^{0} \rangle\langle B^{1}|B^{1} \rangle + 4 \langle A^{0} B^{0} | V | A^{1} B^{1} \rangle , \end{aligned}$$

(A.22)

which written in explicit notation is the same as Eq. (60) in ref. [107]. This formula can be expressed in terms of densities as

$$\begin{aligned} E^{(3)}_{\text{ind}}& = \int \rho^{(11)A}_{\text{ind}}({\boldsymbol r}) \tilde{\omega}^{B}({\boldsymbol r}) d^{3}\boldsymbol{r} + \int \rho^{(11)B}_{\text{ind}}({\boldsymbol r}) \tilde{\omega}^{A}({\boldsymbol r}) d^{3}\boldsymbol{r} \\ & +\iint \rho^{(1)A}_{\text{ind}}({\boldsymbol r}_{1}) \frac{1}{r_{12}} \rho^{(1)B}_{\text{ind}}({\boldsymbol r}_{2}) d^{3}\boldsymbol{r}_{1}d^{3}\boldsymbol{r}_{2} ,\end{aligned}$$

(A.23)

where

$$\tilde{\omega}^{B}({\boldsymbol r}) = \omega^{B}({\boldsymbol r}) - \frac{1}{N_{A}} \int {\rho^{A}_{0}}({\boldsymbol r}) \omega^{B}({\boldsymbol r}) d^{3}\boldsymbol{r} .$$

The third term is the same as in Eq. (25). The first term is numerically the same as the first term in Eq. (25), but the density is now the component of the second-order induction density resulting from the product of the first-order induction functions and it interacts with the shifted electric potential of molecule B, similarly for the second term.

One more interpretation of $E^{(3)}_{\text {ind}}$ can be obtained via regrouping the terms in formula (A.17) and defining an alternative second-order induction function, different from the one of Eq. (A.19). The resulting expression for the third-order induction energy is algorithmically different from Eq. (A.17), but gives the same numerical value of this quantity. The former definition just neglects the second expression on the right-hand side of Eq. (A.19)

$$\begin{aligned}\bar{\Psi}_{\text{ind}}^{(2)A} = \sum\limits_{k\ne 0} \frac{\langle {{\Psi}_{k}^{A}} |\langle {{\Psi}_{0}^{A}} | {\Omega}^{B} {{\Psi}_{0}^{A}} \rangle -{\Omega}^{B}| {\Psi}_{\text{ind}}^{(1)A} \rangle} {{E_{k}^{A}} - {E_{0}^{A}}} {{\Psi}_{k}^{A}} . \end{aligned}$$

(A.24)

The alternative second-order equation is

$$\begin{array}{@{}rcl@{}} (H_{A} - {E_{0}^{A}}) |\bar{A}^{2}\rangle = \left(\langle A^{0}|{\Omega}_{B}|A^{0} \rangle-{\Omega}_{B}\right) |A^{1}\rangle . \end{array}$$

(A.25)

Notice that $|\bar {A}^{2}\rangle$ cannot be called the second-order induction function since the function of Eq. (A.18) is the unique and only such function. We can now design an energy expression

$$E^{(3)}_{\text{ind}} = \langle A^{0}|{\Omega}_{B}|\bar{A}^{2} \rangle + \langle B^{0}|{\Omega}_{A}|\bar{B}^{2} \rangle + 4 \langle A^{0} B^{0} | V | A^{1} B^{1} \rangle$$

(A.26)

or in explicit notation

$$E^{(3)}_{\text{ind}} = \langle {{\Psi}_{0}^{A}} |{\Omega}_{B}| \bar{\Psi}_{\text{ind}}^{(2)A} \rangle + \langle {{\Psi}_{0}^{B}} |{\Omega}_{A}| \bar{\Psi}_{\text{ind}}^{(2)B} \rangle + 4\langle {{\Psi}_{0}^{A}} {{\Psi}_{0}^{B}} |V| {\Psi}_{\text{ind}}^{(1)A} {\Psi}_{\text{ind}}^{(1)B} \rangle ,$$

(A.27)

which produces the same numerical values as given by Eq. (A.17). To show this, transform Eq. (A.26) into Eq. (A.22). To this end, write Eq. (A.25) as

$$|\bar{A}^{2}\rangle = {R_{0}^{A}} \left(\langle A^{0}|{\Omega}_{B}|A^{0} \rangle-{\Omega}_{B}\right) |A^{1}\rangle ,$$

(A.28)

where ${R_{0}^{A}}$ is the reduced resolvent of monomer A, and use it in Eq. (A.26). We get

$$\begin{aligned}\langle A^{0}|{\Omega}_{B}|\bar{A}^{2} \rangle &= \langle A^{0}|{\Omega}_{B} {R_{0}^{A}} \left(\langle A^{0}|{\Omega}_{B}|A^{0} \rangle-{\Omega}_{B}\right) |A^{1}\rangle \\&= \langle A^{1}|{\Omega}_{B}|A^{1} \rangle - \langle A^{0}|{\Omega}_{B}|A^{0} \rangle\langle A^{1}|A^{1} \rangle \end{aligned}$$

(A.29)

and similarly for the second term, which indeed gives the desired formula. Equation (A.27) expressed in terms of densities yields Eq. (25).

Appendix: Conversations

Mo et al. commented:

The authors focus on the elegance of the SRS formulation of SAPT, but perhaps ignore chemistry, and at the same time they are dismissing all other EDA methods as nonunique. They make statements about the uniqueness and precision of SAPT in providing energy partition that sums to the total interaction energy between the molecules/fragments A and B, and dismisses effects like charge transfer, covalency, and charge delocalization which arise from other EDA methods (and also from NBO, VB, and BLW). Our general comment is that there are multiple perspectives in a matter of fact, which are in no way unphysical. What is unscientific is to claim uniqueness and truth for one of these choices, namely SAPT, and to dismiss on this ground all other approaches. This is done without providing the reader with a single example that compares SAPT (e.g., what about BrHBr⁻?) to other EDA methods. In a nutshell the paper is a blizzard of equations without any example. This is a major problem for most chemists, who would like to see examples with numerical data, as proofs of correctness of statements.

Reply:

We do not believe we ignore chemistry, but it depends what one has in mind by “chemistry”. Noncovalent intermolecular interactions are a part of chemistry and our whole paper is devoted to such interactions. Thus, in this sense we cannot agree that we ignore chemistry. On the other hand, we do not discuss making and breaking of chemical bonds since these are not processes that SAPT was designed for (although please see a discussion of this issue later on). In fact, the first sentence in the abstract states “Symmetry-adapted perturbation theory (SAPT) is a method for computational studies of noncovalent interactions between molecules.”

We also do not believe we are dismissing all other EDA methods as nonunique. While this is a plain fact that the EDA methods are highly nonunique, the physical components in some of these methods come reasonably close to the corresponding SAPT components. Such methods are, in our opinion, important since their application on top of some supermolecular calculations of interaction energies does give sufficiently precise physical insights.

Concerning the dismissal of “effects like charge-transfer, covalency, charge delocalization,” indeed, SAPT is not designed to investigate covalency effects as it is practically limited to noncovalent interactions. We make no statements on EDA methods applied to chemically reactive systems. One the other hand, we do not dismiss charge-delocalization effects (which in our terminology are equivalent to charge-transfer effects) and we discuss these effects in “SAPT contributions and EDA methods.” To summarize this discussion: SAPT does include charge-delocalization effects, but it appears there is no unique ways to separate them. However, very reasonable ways to perform such separation approximately have been designed by Misquitta and Stone [140, 142, 143] and analyzed from the point of view of applications in developments of force fields [77].

While we agree that different viewpoints are useful in science, it does not mean that all viewpoints are correct. In particular, if method X states that the dispersion interaction in a given dimer is zero, while method Y states that it is one of major attractive forces, only one of these viewpoints can be correct. Thus, if the authors of method X believe their results is correct, they should explain why method Y makes wrong predictions.

There are so many examples of SAPT interaction energy decompositions in literature that we did not think examples are needed in the present paper. Nevertheless, in the revised version we added Table 1 with such examples. The geometries of the dimers included are available in literature, so these data can perhaps serve as a useful reference point for authors of EDA approaches. We have not included the BrHBr⁻ complex in Table 1. While one can trivially compute SAPT components for this system, it is a system which is to a large effect covalently bound, see the recent work on FHF⁻ dimer [144]. While both systems would be interesting cases to study by SAPT, their special character does not make them appropriate as examples of SAPT analysis of noncovalent interactions. One may add that a somewhat similar system, H₂O-F⁻, is included in ref. [77].

Most published SAPT calculations listed the components of interaction energies, at least for some selected geometries. Here is a selection of such papers for readers who would like to study more examples: water dimer (refs. [20, 22, 145, 146]), helium dimer (refs. [21, 147]), water-uracil (ref. [148]), He-F⁻ (ref. [149]), Ar–H₂ (ref. [150]), Ar–HF (ref. [84]), He-H₂O (ref. [151]), water trimer (ref. [152]), and many other.

Mo et al. commented:

On pages 13–14 the authors state that “one cannot define uniquely the charge transferred from one monomer.” But the fact is that charge transfer accompanies nevertheless many reactions. How does SAPT handles the CT in e.g., S_N2 reaction? As the authors say, CT is handled as damping of other terms of SAPT. It is hard to buy such a statement.

Reply:

We never say that there is no charge delocalization (or transfer) in noncovalent intermolecular interactions, we only say that the amount of charge that was delocalized cannot be uniquely determined. These are two different statements. Again, we make no statements about chemical reactions.

We never say that charge-delocalization energy is due to damping of asymptotic expansion. We just point out that in the method of determining this energy developed by Misquitta and Stone in ref. [140], the numerical value of this quantity depends on assumptions concerning damping.

Mo et al. commented:

On page 14, the authors speak about CT and ask: “what would be the use of such information besides a physical insight?” Let us ask the authors: is there any science without insight?

Reply:

We were not trying to belittle the importance of physical insight, but the major goal of science are predictions about nature. We added an additional sentence at this place to make our views clear.

Mo et al. commented:

On page 14, the authors do recognize charge-transfer-delocalization “and one may assume that this is due to the difficulties of such a basis set to model charge delocalization”. It is good that Stone found a way to add CT to SAPT. There is not much chemistry without CT.

Reply:

We recognize the existence of charge delocalization all the time. In the sentence quoted, we discuss the assumptions of the method proposed by Stone and Misquitta [142] in 2009 and this sentence is not a statement of recognition of charge delocalization. The work of these authors shows how difficult it is to define charge-delocalization energies within SAPT as the three proposed approaches [140, 142, 143] produce quite different numerical values of these quantities. This is one more confirmation of the point we make: charge delocalization is a fact, a unique determination of its energetic effect is not possible, but reasonable approximate definitions of this effect may be useful.

Mo et al. commented:

In fact, MO (or MO-CI – any level) theory tells us that when two molecules approach, there are orbital interactions dominated by one HOMO in one side and one LUMO on the other side. These orbital interactions stabilize the complex by CT interaction that creates covalency. Where is this in SAPT?

Reply:

Again, SAPT is not designed to treat covalent bonds.

Mo et al. commented:

In addition, when two molecules approach one another, their individual molecular orbitals will be perturbed and thus reshuffled. This physical effect is the polarization effect. Thus, the obvious question would be how the SAPT method quantifies the electron (or charge) transfer and polarization energies. The author did address this question in passing in the end of the paper by mentioning the approaches by Stone and Misquitta.

Reply:

The polarization effect, which we call the induction effect, is discussed in detail in “Induction interaction”. We can only repeat one more time that a unique separation of induction energies into parts due to “fixed” (but deformed) and delocalized charge is not possible. Nonunique, approximate separations such as those proposed by Misquitta and Stone [140, 142, 143] can be useful.

Mo et al. commented:

In chemistry, not only inter- but also intramolecule charge transfer plays a significant role and has been well recognized. For example, in benzene (C₆H₆), the delocalization of the six p electrons has a profound influence on molecular structures and properties. It will be helpful for readers if some data from the SAPT computations in this regard can be presented.

Reply:

In SAPT, monomers can be described at various levels of theory, see “Levels of intramonomer electron correlation in SAPT”. Even the lowest possible level, the HF method, should take account of these effects. This is demonstrated by the fact that the SAPT potential for the benzene dimer [153] is the most accurate published one for this system, capable to produce predictions of spectroscopic accuracy [154].

Mo et al. commented:

We are also curious about the basis set dependency of the SAPT method. Taking the simple example of a two-body complex H₃N⋯BH₃, can the author show the SAPT results with the basis sets from 6-31G to 6-311+G(d) to 6-311++G(d,p) to aug-cc-pvtz for this very simple complex? Of course, the comparison of different correlated methods would also be helpful. We believe that the case of H₃N⋯BH₃ with various basis sets would be illuminating if the author is willing to share the computation results.

Reply:

Multiple basis set convergence studies for SAPT components including the basis sets listed above and several other basis sets have been published, see in particular refs. [22] and [129]. We believe no further studies of basis set convergence are needed.

Mo et al. commented:

Getting back to the “Introduction,” the author wrote that “there is no place for terms not present in SAPT since SAPT’s contributions sum up to an accurate value of the interaction energy”. This is quite confusing as all other EDA approaches also sum up all terms to an accurate value of the interaction energy. We do not see disagreements with the interaction energy values, and all controversies come from the interpretation of the energy terms. The accuracy of the SAPT towards the final interaction energies cannot be used as evidence for the correctness of its physical interpretations or energy partition schemes. Again, data can speak better.

Reply:

This statement means the following: virtually all EDA methods identify terms labeled in the same way as in SAPT: electrostatic, first-order exchange, induction, and dispersion energies. In SAPT, these terms sum up to an accurate total interaction energies. Thus, if these four components are defined in an EDA in such a way that they are close to the SAPT values, they also sum up to total interaction energies. Then, no other terms with significant magnitude can be added. We have modified the quoted text to make our point more clear.

Mo et al. commented:

Besides, any theoretical results need be justified by experimental evidence, directly or indirectly. On page 5, the author wrote “In fact, there is no resemblance between SAPT(DFT) and supermolecular DFT interaction energies for majority of dimers”. This is again quite confusing. DFT interaction energies rely on the DFT methods themselves not any particular EDA method. It seems that the author is comparing orange with apple here as SAPT(DFT) and DFT are not at the same theoretical levels.

Reply:

SAPT results are fully confirmed by experimental evidence. The most convincing confirmation are vibration-rotation-tunneling (VRT) spectra of dimers and trimers: spectra computed from SAPT potentials agree very well with experiment [84, 154,155,156,157]. Another example are crystal-structure predictions from SAPT-based force fields, which correctly identify the experimental crystal as one of the top-ranked polymorphs [158,159,160,161,162,163].

A different question is if the individual energy components predicted by SAPT can be related to experimental data. The evidence is less direct here, but it does exist. In the long-range region, the interaction energy of a dimer made from polar monomers is dominated by the electrostatic energy. For such systems, scattering experiments can sometimes identify the so-called long-range entrance-channel states which are located in the regions of strongest electrostatic interactions. Furthermore, the electrostatic energies are guiding crystallographers in designing crystals of polar molecules. The induction energy dominates the long-range total potential in interactions of ions with rare-gas atoms. Therefore scattering experiments on such systems directly probe the induction components of PESs. Similarly, interactions of rare-gas atoms are dominated at long range by the dispersion energies. Thus, the measured s-wave scattering lengths probe dispersion interactions. Finally, the exchange component is related to van der Waals radii of elements. This component determines the repulsive wall of potentials and the repulsive wall in turn determines the van der Waals radii.

The second question is orthogonal to the first. Yes, SAPT(DFT) and supermolecular DFT are different levels of theory. Still, different approaches can produce similarly accurate interaction energies. For example, interaction energies from high-order SAPT based on the FCI description of monomers agree to several digits with supermolecular FCI energies [16, 18, 66, 127]. On the other hand, while SAPT(DFT) gives accurate interaction energies, those from supermolecular DFT, as it is well known, have in general dramatically large errors, mainly due to the fact that semilocal DFT methods do not reproduce dispersion energies at physically relevant intermolecular separations [118, 164]. The underlying reason for this problem is the “shortsightedness” of interelectron interactions in semilocal DFT approximations. However, DFT describes monomers reasonably well. This is why Kohn-Sham monomer determinants and TD-DFT monomer density-density response functions can be used to construct SAPT(DFT) components. The dispersion energy is obtained in SAPT(DFT) from wave-function-type expressions and therefore the DFT shortsightedness problem does not matter. References [26, 113] provide theoretical justifications for high accuracy of SAPT(DFT). It is easy to understand this in the case of electrostatic energy. Most variants of semilocal generalized-gradient approximation (GGA) DFTs give quite accurate electron densities, except at large separations from nuclei. Because of the latter problem, the first SAPT(KS) calculations [109] gave poor electrostatic energies. This problem can be fixed by applying the asymptotic correction as done in refs. [110, 111], leading to densities accurate at all separations. Since the electrostatic energy is just an integral of electron densities of unperturbed monomers, if the densities are accurate, so is the electrostatic energy.

Mo et al. commented:

On Page 2 the authors state that SAPT in higher order is accurate even for diatomics, e.g., LiH. Can he show an example or two? There is no chemistry without delocalization, there is no chemistry without covalency. What about H₂? Where is the covalency in SAPT? What about resonance?

Reply:

The question presumably concerns chemically bound diatomics (for NCIs in diatomics, SAPT is accurate already in the low order). For two diatomics, H₂ and LiH, SAPT was applied to chemically bonded ground states of these systems [15, 16, 18, 54, 87, 165]. Table I in ref. [15] shows that in 60th order SRS recovers the energy of the singlet state of H₂ at R = 2.0 bohr to within 0.0005%. Table III in ref. [18] shows that the best working variant of SAPT recovers the LiH binding energy at the chemical minimum to within 0.0001%. Figure 8 in the same paper shows potential energy curves in the region of chemical minimum computed using 4th-order SAPT. So clearly, SAPT can recover covalent interactions. Yet, because of the necessity to apply a high-order treatment, we do not recommend SAPT to study strong covalent interactions. However, this excellent convergence says nothing about presence or absence of charge delocalization. As already stated several times, SAPT does account for delocalization effects, but cannot separate them from polarization effects. Finally, resonances do appear in intermolecular interactions when the interacting systems are degenerate. It is possible to construct SAPT applicable to such systems [11].

Mo et al. commented:

Some technical questions:

Are the orbitals in A and B orthogonal? Presumably they are not since Pauli repulsion is accounted for during the anti-symmetrization procedure. But Pauli repulsion necessarily bring about electronic effects like CT. Where are these in the SAPT picture? Any example?

Reply:

Orbitals within system A are orthogonal to each other, and similarly for systems B. However, orbitals of A are not orthogonal to orbitals of B. This is independent of antisymmetrization as such non-orthogonality exists already at the level of RSPT. The non-orthogonality is treated in SAPT exactly, i.e., proper orbital overlap integrals appear in all formulas. Indeed, antisymmetrization leads to a distortion of charge density, so the exchange components of SAPT do contain some delocalization, it is not only the induction energy which contributes to charge delocalization. Since delocalization effects are not separable from charge distortion effects that do not involve any shift of charge, no examples can be given.

Mo et al. commented:

Doesn’t SAPT miss one electronic effect like CT because its perturbation Hamiltonian includes only bielectronic Coulombic terms [$H_{0} = H_{A}+H_{B}+ H_{C} + V_{AB}+V_{AC}+V_{BC} + \dots$]?

Reply:

The terms on the right-hand side define the exact Hamiltonian of Schrödinger’s quantum mechanics for atoms and molecules. No three-electrons interactions are present in such Hamiltonians. [BTW, since this is the total Hamiltonian of a cluster, it should be denoted as H rather than H₀ since in the customary notation H₀ = H_A + H_B + H_C.]

Brink and Borrfors commented:

It is very reassuring that SAPT to high orders is formally an exact theory and can handle covalent bond formation. However, the most common SAPT variants are truncated at the second order or possibly the third order and better suited for weaker interactions. In addition, there are other approximations that commonly are employed, such as $\delta E_{\text {int, resp}}^{\text {HF}}$ (eq. 7) and the S² approximation. Is there any approach for estimating the accuracy of the employed SAPT level for a given problem? Can the $\delta E_{\text {int, resp}}^{\text {HF}}$ value be used as such an indicator? For example, does a value lower than a certain number or lower than a certain fraction of the total SAPT energy indicate that the SAPT level is sufficiently accurate for the problem at hand?

Reply:

As discussed above, SAPT is not designed for interactions leading to formation of covalent bonds. For NCIs, second-order treatment works well in practice. In fact, in most SAPT calculations basis set incompleteness errors are larger than SAPT theory-level errors. In a number of papers, SAPT interaction energies were compared to CCSD(T) energies at complete basis set (CBS) limits. Perhaps the most thorough comparison was performed in ref. [115] on 10 dimers and about 100 configurations total. The median unsigned percentage error computed for all dimers in an augmented triple-zeta basis relative to CCSD(T)/CBS was 2.6%. This should be compared to the same error for CCSD(T) in the same basis set amounting to 1.2%.

Estimates of SAPT accuracy by comparisons to other accurate methods are the only reliable ones. The magnitude of $\delta E_{\text {int, resp}}^{\text {HF}}$ is not an indication of the overall error of SAPT. In fact, SAPT performs very well for interactions of strongly polar systems, while $\delta E_{\text {int, resp}}^{\text {HF}}$ is always large for such systems. Thus, other than the average errors such as those found in ref. [115], there are no a priori estimates of the size of SAPT error for a given dimer. In practice, one usually performs CCSD(T) calculation for a couple of points on a potential energy surface to estimates the uncertainties of SAPT, as well as performs a few calculations at the CBS limits to estimate basis set incompleteness errors.

It appears that the errors due to the use of the S² approximation and due to the addition of $\delta E_{\text {int, resp}}^{\text {HF}}$ are smaller than the errors due to the truncation of SAPT expansion at the second order, although there is no study showing this unequivocally. The S² approximation can now be eliminated [124,125,126]; however, it affects the results significantly only at very short separations and the S² errors are mostly removed by the use of $\delta E_{\text {int, resp}}^{\text {HF}}$. The physical reason for including $\delta E_{\text {int, resp}}^{\text {HF}}$ is to account for the induction and exchange-induction effects of the third (or fourth) and higher orders. For polar systems, the advantages of adding the third and higher-order induction effects much outweigh the small inaccuracies [94, 106] that this addition introduces in the first and second order.

Brink and Borrfors commented:

What are the main reasons for the higher accuracy of SAPT(DFT) compared to supermolecular DFT? Is the difference in interaction energy dominated by the more accurate description of the dispersion energy in the former approach or is SAPT(DFT) able to describe other energy contributions more accurately, as well? Is it still possible to determine third-order and higher terms by a similar equation to eq. 7 ($\delta E_{\text {int, resp}}^{\text {HF}}$ approximation) in SAPT(DFT)? A related question concerns the functional dependence of SAPT(DFT). Is SAPT(DFT) less dependent than supermolecular DFT on the choice of DFT functional? Furthermore, which energy term in SAPT is most functional dependent?

Reply:

The reasons that SAPT(DFT) interaction energies are more accurate than the supermolecular DFT ones have already been discussed in the reply to one of the Mo et al. questions. Indeed, the inability of semilocal GGA approaches to recover dispersion interactions is one of the reasons. However, it is not the only reason. An extensive discussion [118] of the other reasons based on analysis of numerical results for several dimers led to the conclusion that inaccuracies originating from DFT components unrelated to the dispersion energy are of similar magnitude. This work analyzed the quantity

$$E_{\text{int}}^{\text{extra}} = E_{\text{int}}^{\text{DFT}} - E_{\text{int}}^{\text{dispersionless}},$$

where $E_{\text {int}}^{\text {DFT}}$ is the supermolecular DFT interaction energy and $E_{\text {int}}^{\text {dispersionless}}$ is a near-exact interaction energy minus the dispersion and exchange-dispersion contribution. If the dispersion energy was the only problem of DFT, it should recover $E_{\text {int}}^{\text {dispersionless}}$ well, i.e., $E_{\text {int}}^{\text {extra}}$ should be small (except possibly at very small intermonomer separations where the electrons of the interacting monomers get into the “visibility” region of DFT). Figures 4 and 6 in ref. [118] show that this is not the case, in fact, the recovery of $E_{\text {int}}^{\text {dispersionless}}$ is poor. So the answer is confirmative: the accurate description of the dispersion interaction in SAPT(DFT) compared to essentially no description in supermolecular DFT is one reason for SAPT(DFT) being so much more accurate, but SAPT(DFT)’s ability to describe the other interaction energy contributions more accurately than does supermolecular DFT is another, perhaps equally important reason.

Yes, the addition of $\delta E_{\text {int, resp}}^{\text {HF}}$ is as rigorous in SAPT(DFT) as in SAPT based on wave-function description of monomers. The reason is that the orders in V in each version of SAPT are rigorously separated from each other. In particular, SAPT(DFT) in its current version includes only first- and second-order terms, while $\delta E_{\text {int, resp}}^{\text {HF}}$ includes only the third- and higher-order terms (plus a small “contamination” in lower orders [94, 106] which it the reason the addition of $\delta E_{\text {int, resp}}^{\text {HF}}$ is an approximation, as discussed above).

Yes, SAPT(DFT) interaction energies change insignificantly when different variants of GGAs are used (provided an asymptotic correction is applied) compared to dramatically different interaction energies from different variants of supermolecular DFT. This issue was investigated in a number of papers [26, 113, 166,167,168,169]. Interestingly, the PBE0 functional [170, 171] shows consistently the best performance in SAPT(DFT) calculations.

The SAPT component most dependent on the choice of the density functional depends on type of interactions. For dispersion-bonded systems like rare-gas dimers, the effect is the largest in absolute terms for the dispersion energies, see for example Table IV in ref. [26]. For dimers of polar monomers, the largest effects come from the first-order and induction energies, see Table V in ref. [26] and Table IV in ref. [113].

Brink and Borrfors commented:

A limitation of SAPT for the analysis of larger systems seems to be the lack of an efficient procedure for structure optimization of molecular complexes. In particular, strong interactions often lead to conformational changes and changes to intramolecular geometry parameters, e.g., intramolecular bond lengths. How are the structure optimizations of such systems best handled? When employing SAPT(DFT) it does not seem advisable to use a supermolecular approach for structure optimization as supermolecular DFT is much less accurate than SAPT(DFT) for intermolecular interactions. Would it be possible to use a mixed approach where supermolecular calculations are used to determine binding conformations and intramolecular parameters and where SAPT is used for refining intermolecular distances? Would such a procedure be sufficiently accurate and can it be automated? A related question concerns the best approach for computing vibrational corrections (zero point and thermal corrections) to complexation enthalpies and free energies?

Reply:

Actually, our programs provide one of the most efficient approaches to structure optimization for molecular complexes. Let us focus first on rigid monomers. In this case, one needs monomers’ structures. To get them, one can use standard electronic structure programs. For not too large monomers (containing up to couple dozen atoms), one can usually determine a small number of starting monomer’s configurations based on chemical intuition. Then, local optimization algorithms (i.e., algorithms finding the minimum closest from the starting point) will reliably find global and local minima (conformers) of each monomer. If one wanted to proceed in this way to find minimum structures of the dimer, this approach would frequently fail since the locations of minima on the potential energy surface of the dimer are often in very nonintuitive places. Thus, many starting points would have to be tried, which makes such optimizations very expensive even at the DFT level. Our approach is to first fit a potential energy surface and then use the fit function to search for minima. Both tasks are performed completely automatically by the autoPES programs [123, 141]. This protocol is the mixed approach mentioned in the question (and yes, the procedure is accurate, robust, and fully automated). For optimization of monomer geometries any method can be used, e.g., MP2, not necessarily DFT (and if DFT is used, it should be a dispersion-corrected DFT approach). Since the latest version of autoPES can develop flexible-monomer potential energy surfaces, one can now optimize full-dimensional dimer structures. This allows investigations of effects of intermolecular interactions on monomer conformations.

As mentioned earlier, SAPT potentials for smaller clusters have been often used to compute VRT spectra of these clusters. Such calculations produce very accurate zero-point energies and give energy levels allowing computations of thermodynamic quantities (see for example ref. [172]). For larger clusters, the potential energy surface can be used to compute the Hessian and proceed in the standard way to obtain thermodynamic quantities in the harmonic approximation.

Popelier commented:

Question 1: Let us take a single molecule, such as Br(CH₂)₁₀Br, and curl back its chain so that the two Br atoms end up in close contact but without being bonded. Surely there is a dispersion-like interaction between the two Br atoms but can SAPT, as presented in this article, calculate its energy value? SAPT’s basic assumption is the partitioning of the total Hamiltonian into a sum of Hamiltonians of separated monomers. However, Br(CH₂)₁₀Br is a single molecule and cannot be separated into monomers. Is there a conceptual challenge in partitioning the Hamiltonian for this covalently bound system, in the absence of monomers? Is a way forward the “atomic SAPT partition” or A-SAPT (J. Chem. Phys. 2014, 141, 044115)? However, A-SAPT struggles to produce chemically useful partitions of the electrostatic energy, caused by the buildup of oscillating partial charges on adjacent functional groups. This is why “functional-group SAPT” or F-SAPT (J. Chem. Theor. Comp. 2014, 10, 4417) was proposed. But then F-SAPT is formulated entirely in terms of fragments with integral charge (including zero), which may suit this molecule but which is not realistic in general.

Reply:

Yes, SAPT is a theory starting with the assumption that the system separates into a set of monomers at infinite distances from each other and these monomers are well-defined molecules (not necessarily closed-shell). As discussed above, this separation should not involve any breaking of chemical bonds. Thus, standard SAPT cannot be applied to Br(CH₂)₁₀Br.

A-SAPT and F-SAPT have been proposed with the goal to approximately assign interaction energy contributions to atoms or groups of atoms in the standard SAPT approach involving dimers of two closed-shell monomers. An extension of SAPT to intramonomer NCIs, named ISAPT, was proposed later [173]. Another approach of this type was developed in ref. [174]. One should emphasize that all intramolecular SAPT applications require one to make several assumptions and that different but equally reasonable assumptions can lead to very different predictions for a given system.

Popelier commented:

Question 2: It is stated that SAPT defines energy contributions each of which results from a differential equation that has an exact solution. Please give examples of such differential equations as they do not seem to appear in the standard SAPT literature.

Reply:

The differential equations for the wave function corrections are the foundations of RSPT. Such equations appear in many papers developing SAPT, see for example Eqs. (5), (10), and (27) in ref. [21] and Eqs. (18), (29), (A.4), and (A.12) in the present work.

Popelier commented:

Question 3: Does the author agree with the opinion of Konrad Patkowski who writes in reference [58] that “my personal least favorite SAPT term is the $\delta E^{(2)}_{HF}$ correction of Equation (20) (...). I consider its presence as an admission that pure SAPT has a difficulty that cannot be fully resolved from within, and it requires outside help in a form of supermolecular HF”? Note that his Eq. (20) is the same as Eq. (7) in the current article and hence $\delta E^{(2)}_{HF} = \delta E_{\text {int, resp}}^{\text {HF}}$.

Reply:

As already discussed above, from the practical point of view the addition of $\delta E_{\text {int, resp}}^{\text {HF}}$ poses no problems. It does account for higher-order induction and exchange-induction effects and the errors introduced by this addition are very small. Furthermore, for nonpolar and mildly polar systems the addition of $\delta E_{\text {int, resp}}^{\text {HF}}$ is not needed if the third-order SAPT interaction energies are computed [108]. Thus, the issue is more of aesthetic than practical nature. It may be possible to increase the range of systems that do not need $\delta E_{\text {int, resp}}^{\text {HF}}$ by removing some approximations in the present set of third-order terms. Another possible step in this direction is to apply the formula for the second-order induction wave function derived in the present work in computations of the fourth- and fifth-order induction energy corrections. One more possible avenue is to make the regularized SAPT [18, 88] applicable to general monomers (the regularized SAPT exhibits a faster convergence of induction energies).

Popelier commented:

Question 4: It is stated that “One cannot uniquely determine the total charge transferred from monomer A to monomer B since this requires choosing an arbitrary boundary between the monomers.” The claim that charge cannot be uniquely assigned to atoms or even molecules is typically perpetuated, yet there is a great need to do so both in terms of interpretative (quantum) chemistry and force field construction. Which criterion (or criteria) give(s) rise to uniqueness if it is not experimental arbitration? According to this article, present-day SAPT (“the theory of intermolecular forces (...) providing the ‘standard model’ for EDA methods”) is declared unique because the symmetrized Rayleigh-Schrödinger method is the only one used in practice. With such perhaps relaxed uniqueness criterion, can Occam’s razor not be used to propose the topological partitioning as a satisfactory method to settle the debate on how to quantify charge transfer (even at the level of tens of millielectrons)? Moreover, an extensive and thorough comparison between fuzzy (interpenetrating) and non-fuzzy (space-filling, e.g., QTAIM) partitioning methods (J. Comp. Chem. 2007, 28, 161) showed that the latter “may be preferred from the chemical consistency point of view” as they also ”better preserve the atomic or fragment identity from the energetic point of view”.

Reply:

First, we changed the terminology to “reference model”. Nevertheless, we maintain that SAPT components are the quantities that EDA methods should compare to (and mostly do). As discussed in the reply to Mo et al., there is, actually, reasonably convincing experimental evidence for the physical character of SAPT components.

Also in the reply to Mo et al., we have extensively discussed our position on charge delocalization. Of course, AIM-atoms give a possible definition of charge delocalization, but still do not allow to determine its energetic effect using SAPT.

Popelier commented:

Question 5: Can SAPT match the quantification of steric effects that the Interacting Quantum Atoms (IQA) method is able to achieve (J. Phys. Chem. A 2016, 120, 9647; ChemPhysChem 2021, 22, 775; Chemistry Open 2019, 8, 560)?

Reply:

As stated in the quote from Richard Feynman in the “Introduction,” molecules are “repelling upon being squeezed into one another.” A simple interpretation relates this repulsion to steric effects: since atoms are approximately spherical, a molecule can be viewed as a shape resulting from superposition of such spheres, and if two interacting molecules come close enough to each other, the penetration of the spheres results in a repulsion. The more the spheres overlap, the larger the repulsive component becomes. In force fields, the repulsive effects are approximated by exponential terms or by large inverse powers of R (mostly 1/R¹²) with positive coefficients. SAPT, of course, accounts for steric effects, as demonstrated by the fact that SAPT PESs are accurate in the repulsive regions. The steric effects computed using SAPT cannot be compared with those computed using IQA as the two quantities are defined in different and unrelated ways. We would like to reiterate here that, as it is extensively discussed in “Exchange interactions”, SAPT provides physical interpretation of the steric interactions grounded in quantum mechanics. The main repulsive effect comes from the first-order exchange energy, Eq. (43). Although, as discussed in “Exchange interactions”, the exchange interactions are due to electron tunneling between monomers and not to overlap of wave functions, this tunneling is proportional to such overlap and can be modelled using overlap integrals [175,176,177]. The exchange interactions appear also in the second-order of SAPT, see Eqs. (46) and (47). These exchange interactions are repulsive as well, but significantly smaller than the first-order ones. They are also proportional to overlap integrals. Another group of effects related to the wave-function overlap are damping effects [53, 79,80,81,82]. In most cases, damping leads to positive contributions to interaction energies; however, as discussed briefly in the “Introduction,” it can also produce negative contributions, which then violate the simple steric interactions picture. Since the first-order exchange energy dominates the repulsive interactions, the values of this component presented in Table 1 provide a reasonable approximation of the total steric effect.

An example of unusual steric effects on shapes of PESs is provided in a recent study of the ammonia dimer [178]. This PES is very unusual in that it has a very narrow canyon-like valley where two equivalent minima are located, with a very small barrier between them at the lowest-energy saddle point. An overlap-driven variation of the components in the direction perpendicular to the interconversion path through this saddle, i.e., across the valley, explains the narrowness of this valley. One has to consider both the exchange and damping effects to get a complete explanation.