Symmetry-adapted perturbation theory (SAPT) is a theory of intermolecular
interactions. When computing intermolecular interaction energies one typically
computes the energy of two molecules infinitely separated and in contact, then
computes the interaction energy by subtraction. SAPT, in contrast, is a
perturbative expression for the interaction energy itself. The various terms
in the perturbation series are physically meaningful, and this decomposition of
the interaction energy can aid in the interpretation of the results. A brief
overview of the theory is given below; for additional technical details, the
reader is referred to Jeziorski *et al.*.^{424, 425}
Additional context can be found in a pair of more recent review
articles.^{897, 385}

In SAPT, the Hamiltonian for the $\mathrm{A}\mathrm{\cdots}\mathrm{B}$ dimer is written as

$$\widehat{H}={\widehat{F}}^{A}+{\widehat{F}}^{B}+\xi {\widehat{W}}^{A}+\eta {\widehat{W}}^{B}+\zeta \widehat{V},$$ | (13.37) |

where ${\widehat{W}}^{A}$ and ${\widehat{W}}^{B}$ are Møller-Plesset fluctuation operators for fragments $A$ and $B$, whereas $\widehat{V}$ consists of the intermolecular Coulomb operators. This part of the perturbation is conveniently expressed as

$$\widehat{V}=\sum _{i\in A}\sum _{j\in B}\widehat{v}(ij)$$ | (13.38) |

with

$$\widehat{v}(ij)=\frac{1}{\left|{\overrightarrow{r}}_{i}-{\overrightarrow{r}}_{j}\right|}+\frac{{\widehat{v}}_{A}(j)}{{N}_{A}}+\frac{{\widehat{v}}_{B}(i)}{{N}_{B}}+\frac{{V}_{0}}{{N}_{A}{N}_{B}}.$$ | (13.39) |

The quantity ${V}_{0}$ is the nuclear interaction energy between the two fragments and

$${\widehat{v}}_{A}(j)=-\sum _{I\in A}\frac{{Z}_{I}}{\left|{\overrightarrow{r}}_{j}-{\overrightarrow{R}}_{I}\right|}$$ | (13.40) |

describes the interaction of electron $j\in B$ with nucleus $I\in A$.

Starting from a zeroth-order Hamiltonian ${\widehat{H}}_{0}={\widehat{F}}^{A}+{\widehat{F}}^{B}$ and zeroth-order
wave functions that are direct products of monomer wave functions,
$|{\mathrm{\Psi}}_{0}\u27e9=|{\mathrm{\Psi}}_{A}\u27e9|{\mathrm{\Psi}}_{B}\u27e9$, the SAPT
approach is based on a symmetrized Rayleigh-Schrödinger perturbation
expansion^{424, 425} with respect to the perturbation parameters
$\xi $, $\eta $, and $\zeta $ in Eq. (13.37). The resulting interaction energy can
be expressed as^{424, 425}

$${E}_{\mathrm{int}}=\sum _{i=1}^{\mathrm{\infty}}\sum _{j=0}^{\mathrm{\infty}}\left({E}_{\mathrm{pol}}^{(ij)}+{E}_{\mathrm{exch}}^{(ij)}\right).$$ | (13.41) |

Because it makes no sense to treat ${\widehat{W}}^{A}$ and ${\widehat{W}}^{B}$ at different
orders of perturbation theory, there are only two indices in this expansion:
$j$ for the monomer fluctuations potentials and $i$ for the intermolecular
perturbation. The terms ${E}_{\mathrm{pol}}^{(ij)}$ are known collectively as the
*polarization expansion*, and these are precisely the same terms that would
appear in ordinary Rayleigh-Schrödinger perturbation theory, which is valid
when the monomers are well-separated. The polarization expansion contains
electrostatic, induction and dispersion interactions, but in the
*symmetrized* Rayleigh-Schrödinger expansion, each term ${E}_{\mathrm{pol}}^{(ij)}$
has a corresponding exchange term, ${E}_{\mathrm{exch}}^{(ij)}$, that arises from an
anti-symmetrizer ${\widehat{\mathcal{A}}}_{AB}$ that is introduced in order to project
away the Pauli-forbidden components of the interaction energy that would
otherwise appear.^{425}

The version of SAPT that is implemented in Q-Chem assumes that $\xi =\eta =0$,
an approach that is usually called SAPT0.^{385} Within the
SAPT0 formalism, the interaction energy is formally expressed by the following
symmetrized Rayleigh-Schrödinger expansion:^{424, 425}

$${E}_{\mathrm{int}}(\zeta )=\frac{\u27e8{\mathrm{\Psi}}_{0}|\zeta \widehat{V}{\widehat{\mathcal{A}}}_{AB}|\mathrm{\Psi}(\zeta )\u27e9}{\u27e8{\mathrm{\Psi}}_{0}|{\widehat{\mathcal{A}}}_{AB}|\mathrm{\Psi}(\zeta )\u27e9},$$ | (13.42) |

The anti-symmetrizer ${\widehat{\mathcal{A}}}_{AB}$ in this expression can be written as

$${\widehat{\mathcal{A}}}_{AB}=\frac{{N}_{A}!{N}_{B}!}{({N}_{A}+{N}_{B})!}{\widehat{\mathcal{A}}}_{A}{\widehat{\mathcal{A}}}_{B}\left(\widehat{1}+{\widehat{\mathcal{P}}}^{AB}+{\widehat{\mathcal{P}}}^{\prime}\right),$$ | (13.43) |

where ${\widehat{\mathcal{A}}}_{A}$ and ${\widehat{\mathcal{A}}}_{B}$ are anti-symmetrizers for the two
monomers and ${\widehat{\mathcal{P}}}^{AB}$ is a sum of all one-electron exchange
operators between the two monomers. The operator ${\widehat{\mathcal{P}}}^{\prime}$ in
Eq. (13.43) denotes all of the three-electron and
higher-order exchanges. This operator is neglected in what is known as the
“single-exchange” approximation,^{424, 425} which
is expected to be quite accurate at typical van der Waals and larger
intermolecular separations, but sometimes breaks down at smaller intermolecular
separations.^{523}

Only terms up to $\zeta =2$ in Eq. (13.42)—that is, second order in the intermolecular interaction—have been implemented in Q-Chem. It is common to relabel these low-order terms in the following way [cf. Eq. (13.41)]:

$${E}_{\mathrm{int}}^{\mathrm{SAPT0}}={E}_{\mathrm{elst}}^{(1)}+{E}_{\mathrm{exch}}^{(1)}+{E}_{\mathrm{pol}}^{(2)}+{E}_{\mathrm{disp}}^{(2)}.$$ | (13.44) |

The electrostatic part of the first-order energy correction is denoted
${E}_{\mathrm{elst}}^{(1)}$ and represents the Coulomb interaction between the two
monomer electron densities.^{425} The quantity
${E}_{\mathrm{exch}}^{(1)}$ is the corresponding first-order (*i.e.*, Hartree-Fock)
exchange correction. Explicit formulas for these corrections can be found in
Ref. 424. The second-order term from the polarization
expansion, denoted ${E}_{\mathrm{pol}}^{(2)}$ in Eq. (13.44),
consists of a dispersion contribution (which arises for the first time at
second order) as well as a second-order correction for induction. The latter can be written

$${E}_{\mathrm{ind}}^{(2)}={E}_{\mathrm{ind}}^{(2)}(A\leftarrow B)+{E}_{\mathrm{ind}}^{(2)}(B\leftarrow A),$$ | (13.45) |

where the notation $A\leftarrow B$, for example, indicates that the frozen charge density of $B$ polarizes the density of $A$. In detail,

$${E}_{\mathrm{ind}}^{(2)}(A\leftarrow B)=2\sum _{ar}{t}_{ar}{({w}_{B})}_{ra}$$ | (13.46) |

where

$${({w}_{B})}_{ar}={({\widehat{v}}_{B})}_{ar}+\sum _{b}(ar|bb)$$ | (13.47) |

and ${t}_{ar}={({w}_{B})}_{ar}/({\u03f5}_{a}-{\u03f5}_{r})$. The second term in
Eq. (13.45), in which $A$ polarizes $B$, is obtained by
interchanging labels.^{414} Finally, the second-order dispersion
correction has a form reminiscent of the MP2 correlation energy:

$${E}_{\mathrm{disp}}^{(2)}=4\sum _{abrs}\frac{(ar|bs)(ra|sb)}{{\u03f5}_{a}+{\u03f5}_{b}-{\u03f5}_{r}-{\u03f5}_{s}}.$$ | (13.48) |

The induction and dispersion corrections both have accompanying exchange
corrections (exchange-induction and
exchange-dispersion).^{424, 425}

The similarity between Eq. (13.48) and the MP2 correlation energy means that SAPT jobs, like MP2 calculations, can be greatly accelerated using resolution-of-identity (RI) techniques, and an RI version of SAPT is available in Q-Chem. To use it, one must specify an auxiliary basis set. The same ones used for RI-MP2 work equally well for RI-SAPT, but one should always select the auxiliary basis set that is tailored for use with the primary basis of interest, as in the RI-MP2 examples in Section 6.6.1.

It is common to replace ${E}_{\mathrm{ind}}^{(2)}$ and ${E}_{\mathrm{exch}\text{-}\mathrm{ind}}^{(2)}$ in
Eq. (13.44) with their “response” (“resp”) analogues,
which are the infinite-order correction for polarization arising from a frozen partner
density.^{424, 425} Operationally, this
substitution involves replacing the second-order induction amplitudes, ${t}_{ar}$
in Eq. (13.46),
with amplitudes obtained from solution of the coupled-perturbed Hartree-Fock equations.^{741}
(The perturbation is simply the electrostatic potential of the other monomer.) In addition,
it is common to correct the SAPT0 binding energy for higher-order polarization effects
by adding a correction term of the form^{425, 385}

$$\delta {E}_{\mathrm{int}}^{\mathrm{HF}}={E}_{\mathrm{int}}^{\mathrm{HF}}-\left({E}_{\mathrm{elst}}^{(1)}+{E}_{\mathrm{exch}}^{(1)}+{E}_{\mathrm{ind},\mathrm{resp}}^{(2)}+{E}_{\mathrm{exch}\text{-}\mathrm{ind},\mathrm{resp}}^{(2)}\right)$$ | (13.49) |

to the interaction energy. Here, ${E}_{\mathrm{int}}^{\mathrm{HF}}$ is the counterpoise-corrected Hartree-Fock binding energy for $\mathrm{A}\mathrm{\cdots}\mathrm{B}$. Both the response corrections and the $\delta {E}_{\mathrm{int}}^{\mathrm{HF}}$ correction are optionally available in Q-Chem’s implementation of SAPT.

It is tempting to replace Hartree-Fock MOs and eigenvalues in the SAPT0
formulas with their Kohn-Sham counterparts, as a low-cost means of introducing
monomer electron correlation. The resulting procedure is known as
SAPT(KS),^{988} and does offer an improvement on SAPT0 for some
strongly hydrogen-bonded systems.^{369} Unfortunately, SAPT(KS)
results are generally in poor agreement with benchmark dispersion
energies,^{369} owing to incorrect asymptotic behavior of
approximate exchange-correlation potentials.^{638} The
dispersion energies can be greatly improved through the use of long-range
corrected (LRC) functionals in which the range-separation parameter, $\omega $,
is “tuned” so as to satisfy the condition ${\u03f5}_{\mathrm{HOMO}}=-\text{IE}$, where ${\u03f5}_{\mathrm{HOMO}}$ is the HOMO energy and “IE"
represents the ionization energy.^{526} Monomer-specific values
of $\omega $, tuned using the individual monomer IEs, substantially improve
SAPT(KS) dispersion energies, though the results are still not of benchmark
quality.^{526} Other components of the interaction energy, however,
can be described quite accurately SAPT(KS) in conjunction with a tuned version
of LRC-$\omega $PBE.^{526} Use of monomer-specific $\omega $ values is
controlled by the variable DFT-LRC in the *$xpol* section, as described in
Section 13.11.3, with monomer-specific values of $\omega $
that must be specified in a *$lrc_omega* input section.
As an example, some values of $\omega $ for various monomers obtained using the
IE-tuning condition and the LRC-$\omega $PBE functional are listed in
Table 13.1. Clearly there is a non-trivial variation amongst the
optimally-tuned values of $\omega $.

S22 monomers | S66 monomers | ions | |||
---|---|---|---|---|---|

Monomer | $\omega $ / ${a}_{0}^{-1}$ | Monomer | $\omega $ / ${a}_{0}^{-1}$ | Monomer | $\omega $ / ${a}_{0}^{-1}$ |

adenine | 0.271 | ${\mathrm{MeNH}}_{2}$ | 0.397 | F${}^{-}$ | 0.480 |

2-aminopyridine | 0.293 | MeOH | 0.438 | Cl${}^{-}$ | 0.372 |

benzene | 0.280 | ${\mathrm{AcNH}}_{2}$ | 0.453 | SO${}_{4}^{2-}$ | 0.344 |

ethyne | 0.397 | $\mathrm{AcOH}$ | 0.381 | Li${}^{+}$ | 2.006 |

ethene | 0.359 | cyclopentane | 0.420 | Na${}^{+}$ | 1.049 |

methane | 0.454 | neopentane | 0.287 | K${}^{+}$ | 0.755 |

formamide | 0.460 | pentane | 0.365 | ||

formic acid | 0.412 | peptide | 0.341 | ||

water | 0.502 | pyridine | 0.316 | ||

HCN | 0.452 | ||||

indole | 0.267 | ||||

ammonia | 0.440 | ||||

phenol | 0.292 | ||||

pyrazine | 0.367 | ||||

2-pyridoxine | 0.294 | ||||

thymine | 0.284 | ||||

uracil | 0.295 |

Finally, some discussion of basis sets is warranted. Typically, SAPT calculations are performed in the
so-called dimer-centered basis set (DCBS),^{989} which means that the combined $A+B$
basis set is used to calculate the zeroth-order wave functions for both $A$ and $B$.
This leads to the unusual situation that there are more MOs than basis functions: one set of occupied
and virtual MOs for each monomer, both expanded in the same (dimer) AO basis. As an alternative to
the DCBS, one might calculate $|{\mathrm{\Psi}}_{A}\u27e9$ using only $A$’s basis functions (similarly
for $B$), in which case the SAPT calculation is said to employ the monomer-centered basis set
(MCBS).^{989} However, MCBS results are generally of poorer quality.
As an efficient alternative to the DCBS, Jacobson and Herbert^{414}
introduced a projected (“proj”) or “pseudocanonicalized” basis set, borrowing
an idea from dual-basis MP2 calculations.^{866} In this approach, the SCF iterations
are performed in the MCBS but then Fock matrices for
fragments $A$ and $B$ are constructed in the dimer ($A+B$) basis set and then
pseudocanonicalized, meaning that the occupied-occupied and
virtual-virtual blocks of these matrices are diagonalized. This procedure does not mix occupied
and virtual orbitals, and thus
leaves the fragment densities and zeroth-order fragment energies unchanged. However, it
does provide a larger set
of virtual orbitals that extend over the partner fragment. This larger virtual space is then
used to evaluate the
perturbative corrections. All three of these basis options (MCBS, DCBS, and projected
basis) are available in Q-Chem.