13.12 Symmetry-Adapted Perturbation Theory (SAPT)

13.12.1 Theory

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..427, 428 Additional context can be found in a pair of more recent review articles.906, 388

In SAPT, the Hamiltonian for the AB dimer is written as

H^=F^A+F^B+ξW^A+ηW^B+ζV^, (13.37)

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

V^=iAjBv^(ij) (13.38)

with

v^(ij)=1|ri-rj|+v^A(j)NA+v^B(i)NB+V0NANB. (13.39)

The quantity V0 is the nuclear interaction energy between the two fragments and

v^A(j)=-IAZI|rj-RI| (13.40)

describes the interaction of electron jB with nucleus IA.

Starting from a zeroth-order Hamiltonian H^0=F^A+F^B and zeroth-order wave functions that are direct products of monomer wave functions, |Ψ0=|ΨA|ΨB, the SAPT approach is based on a symmetrized Rayleigh-Schrödinger perturbation expansion427, 428 with respect to the perturbation parameters ξ, η, and ζ in Eq. (13.37). The resulting interaction energy can be expressed as427, 428

Eint=i=1j=0(Epol(ij)+Eexch(ij)). (13.41)

Because it makes no sense to treat W^A and 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 Epol(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 Epol(ij) has a corresponding exchange term, Eexch(ij), that arises from an anti-symmetrizer 𝒜^AB that is introduced in order to project away the Pauli-forbidden components of the interaction energy that would otherwise appear.428

The version of SAPT that is implemented in Q-Chem assumes that ξ=η=0, an approach that is usually called SAPT0.388 Within the SAPT0 formalism, the interaction energy is formally expressed by the following symmetrized Rayleigh-Schrödinger expansion:427, 428

Eint(ζ)=Ψ0|ζV^𝒜^AB|Ψ(ζ)Ψ0|𝒜^AB|Ψ(ζ), (13.42)

The anti-symmetrizer 𝒜^AB in this expression can be written as

𝒜^AB=NA!NB!(NA+NB)!𝒜^A𝒜^B(1^+𝒫^AB+𝒫^), (13.43)

where 𝒜^A and 𝒜^B are anti-symmetrizers for the two monomers and 𝒫^AB is a sum of all one-electron exchange operators between the two monomers. The operator 𝒫^ 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,427, 428 which is expected to be quite accurate at typical van der Waals and larger intermolecular separations, but sometimes breaks down at smaller intermolecular separations.527

Only terms up to ζ=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)]:

EintSAPT0=Eelst(1)+Eexch(1)+Epol(2)+Edisp(2). (13.44)

The electrostatic part of the first-order energy correction is denoted Eelst(1) and represents the Coulomb interaction between the two monomer electron densities.428 The quantity Eexch(1) is the corresponding first-order (i.e., Hartree-Fock) exchange correction. Explicit formulas for these corrections can be found in Ref. 427. The second-order term from the polarization expansion, denoted Epol(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

Eind(2)=Eind(2)(AB)+Eind(2)(BA), (13.45)

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

Eind(2)(AB)=2artar(wB)ra (13.46)

where

(wB)ar=(v^B)ar+b(ar|bb) (13.47)

and tar=(wB)ar/(ϵa-ϵr). The second term in Eq. (13.45), in which A polarizes B, is obtained by interchanging labels.417 Finally, the second-order dispersion correction has a form reminiscent of the MP2 correlation energy:

Edisp(2)=4abrs(ar|bs)(ra|sb)ϵa+ϵb-ϵr-ϵs. (13.48)

The induction and dispersion corrections both have accompanying exchange corrections (exchange-induction and exchange-dispersion).427, 428

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 Eind(2) and Eexch-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.427, 428 Operationally, this substitution involves replacing the second-order induction amplitudes, tar in Eq. (13.46), with amplitudes obtained from solution of the coupled-perturbed Hartree-Fock equations.749 (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 form428, 388

δEintHF=EintHF-(Eelst(1)+Eexch(1)+Eind,resp(2)+Eexch-ind,resp(2)) (13.49)

to the interaction energy. Here, EintHF is the counterpoise-corrected Hartree-Fock binding energy for AB. Both the response corrections and the δEintHF 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),999 and does offer an improvement on SAPT0 for some strongly hydrogen-bonded systems.372 Unfortunately, SAPT(KS) results are generally in poor agreement with benchmark dispersion energies,372 owing to incorrect asymptotic behavior of approximate exchange-correlation potentials.643 The dispersion energies can be greatly improved through the use of long-range corrected (LRC) functionals in which the range-separation parameter, ω, is “tuned” so as to satisfy the condition ϵHOMO=-IE, where ϵHOMO is the HOMO energy and “IE" represents the ionization energy.530 Monomer-specific values of ω, tuned using the individual monomer IEs, substantially improve SAPT(KS) dispersion energies, though the results are still not of benchmark quality.530 Other components of the interaction energy, however, can be described quite accurately SAPT(KS) in conjunction with a tuned version of LRC-ωPBE.530 Use of monomer-specific ω values is controlled by the variable DFT-LRC in the $xpol section, as described in Section 13.11.3, with monomer-specific values of ω that must be specified in a $lrc_omega input section. As an example, some values of ω for various monomers obtained using the IE-tuning condition and the LRC-ωPBE functional are listed in Table 13.1. Clearly there is a non-trivial variation amongst the optimally-tuned values of ω.

S22 monomers S66 monomers ions
Monomer ω / a0-1 Monomer ω / a0-1 Monomer ω / a0-1
adenine 0.271 MeNH2 0.397 F- 0.480
2-aminopyridine 0.293 MeOH 0.438 Cl- 0.372
benzene 0.280 AcNH2 0.453 SO2-4 0.344
ethyne 0.397 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
Table 13.1: Tuned values of the range separation parameter (ω) for monomers in the S22447 and S66785 data sets, along with a few ions.

Finally, some discussion of basis sets is warranted. Typically, SAPT calculations are performed in the so-called dimer-centered basis set (DCBS),1000 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 |ΨA 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).1000 However, MCBS results are generally of poorer quality. As an efficient alternative to the DCBS, Jacobson and Herbert417 introduced a projected (“proj”) or “pseudocanonicalized” basis set, borrowing an idea from dual-basis MP2 calculations.874 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.