<jupyter_start><jupyter_text>Symmetry-Adapted Perturbation Theory (SAPT0)Symmetry-adapted perturbation theory (SAPT) is a perturbation theory aimed specifically at calculating the interaction energy between two molecules. Compared to the more conventional supermolecular approach where the interaction energy is computed as the difference between the electronic energy of the complex and the sum of electronic energies for the individual molecules (monomers), $E_{\rm int}=E_{\rm AB}-E_{\rm A}-E_{\rm B}$, SAPT obtains the interaction energy directly - no subtraction of similar terms is needed. Even more important, the result is obtained as a sum of separate corrections accounting for the electrostatic, induction, dispersion, and exchange contributions to interaction energy, so the SAPT decomposition facilitates the understanding and physical interpretation of results. In the wavefunction-based variant presented here [Jeziorski:1994], SAPT is actually a triple perturbation theory. The zeroth-order Hamiltonian is the sum of the monomer Fock operators, $H_0=F_{\rm A}+F_{\rm B}$, and the perturbation correction $E^{(nkl)}$ corresponds to $n$th, $k$th, and $l$th order effects, respectively, of the intermolecular interaction operator $V$, the monomer-A Moller-Plesset fluctuation potential $W_{\rm A}=H_{\rm A}-F_{\rm A}$, and an analogous monomer-B potential $W_{\rm B}$. Thus, the SAPT correction $E^{(nkl)}$ is of the $n$th order in the *intermolecular interaction* and of the $(k+l)$th order in the *intramolecular correlation*.In this example, we will calculate the interaction energy between two molecules at the simplest, SAPT0 level of theory [Parker:2014]. In SAPT0, intramolecular correlation is neglected, and intermolecular interaction is included through second order. Specifically,\begin{equation}E_{\rm int}^{\rm SAPT0}=E^{(100)}_{\rm elst}+E^{(100)}_{\rm exch}+E^{(200)}_{\rm ind,resp}+E^{(200)}_{\rm exch-ind,resp}+E^{(200)}_{\rm disp}+E^{(200)}_{\rm exch-disp}\end{equation}In this equation, the consecutive corrections account for the electrostatic, first-order exchange, induction, exchange induction, dispersion, and exchange dispersion effects, respectively. The additional subscript ``resp'' denotes that these corrections are computed including the monomer relaxation (response) effects at the coupled-perturbed Hartree-Fock (CPHF) level of theory.Before we proceed to the computation of the individual SAPT0 corrections, let us make two comments on the specifics of the calculation of the exchange corrections. The exchange terms stem from the symmetry adaptation, specifically, from the presence of the $(N_{\rm A}+N_{\rm B})$-electron antisymmetrizer ${\cal A}$ that enforces the antisymmetry of the wavefunction upon an interchange of a pair of electrons between the monomers. Typically, the full operator ${\cal A}$ in SAPT is approximated as ${\cal A}=1+{\cal P}$, where the *single-exchange operator* ${\cal P}=\sum_{a\in {\rm A}}\sum_{b\in {\rm B}}P_{ab}$ collects all transpositions of a single pair of electrons between the interacting molecules. This approach is known as the *single exchange approximation* or the *$S^2$ approximation* --- the latter name refers to keeping terms that are quadratic in the intermolecular overlap integrals $S$ and neglecting terms that vanish like $S^4$, $S^6$, $\ldots$. In Psi4,the $E^{(100)}_{\rm exch}$ correction can be computed without the $S^2$ approximation, and the nonapproximated formulas for $E^{(200)}_{\rm exch-ind,resp}$ and $E^{(200)}_{\rm exch-disp}$ have also been derived [Schaffer:2013]. Nevertheless, in this example we will employ the $S^2$ approximation in all exchange corrections. Second, there exist two formalisms for the derivation of SAPT exchange corrections: the second-quantization approach [Moszynski:1994a] and the density matrix formalism [Moszynski:1994b]. The two methodologies lead to completely different SAPT expressions which, however, lead to identical results as long as the full dimer basis set is employed. Below, we will adopt the density formalism that is more general (valid in dimer and monomer basis sets) and exhibits more favorable computational scaling (however, more different types of two-electron integrals are required).  1. Preparation of the matrix elementsThe expressions for SAPT0 corrections contain similar quantities as the ones for other correlated electronic structure theories: one- and two-electron integrals over molecular orbitals (MOs), Hartree-Fock (HF) orbital energies, and various amplitudes and intermediates. The feature unique to SAPT is that one has two sets of occupied and virtual (unoccupied) MOs, one for molecule A and one for molecule B (the MOs for the two molecules are not mutually orthogonal, and they may span the same one-electron space but do not have to do so). The most direct consequence of having two sets of MOs is a large number of different MO-basis two-electron integrals $(xy\mid zw)$: each of the four indices can be an occupied orbital of A, a virtual orbital of A, an occupied orbital of B, or a virtual orbital of B. Even when we account for all possible index symmetries, a few dozen types of MO integrals are possible, and we need a code for the integral transformation from atomic orbitals (AOs) to MOs that can produce all of these types. This transformation, and a number of other useful routines, is present in the `helper_SAPT` module that one has to load at the beginning of the SAPT run.<jupyter_code># A simple Psi 4 input script to compute SAPT interaction energies
#
# Created by: <NAME>
# Date: 12/1/14
# License: GPL v3.0
#

import time
import numpy as np
from helper_SAPT import *
np.set_printoptions(precision=5, linewidth=200, threshold=2000, suppress=True)
import psi4

# Set Psi4 & NumPy Memory Options
psi4.set_memory('2 GB')
psi4.core.set_output_file('output.dat', False)

numpy_memory = 2<jupyter_output><empty_output><jupyter_text>Next, we specify the geometry of the complex (in this example, it will be the water dimer). Note that we have to let Psi4 know which atoms belong to molecule A and which ones are molecule B. We then call the `helper_SAPT` function to initialize all quantities that will be needed for the SAPT corrections. In particular, the HF calculations will be performed for molecules A and B separately, and the two sets of orbital energies and MO coefficients will be waiting for SAPT to peruse.<jupyter_code># Set molecule to dimer
dimer = psi4.geometry("""
O   -0.066999140   0.000000000   1.494354740
H    0.815734270   0.000000000   1.865866390
H    0.068855100   0.000000000   0.539142770
--
O    0.062547750   0.000000000  -1.422632080
H   -0.406965400  -0.760178410  -1.771744500
H   -0.406965400   0.760178410  -1.771744500
symmetry c1
""")

psi4.set_options({'basis': 'jun-cc-pVDZ',
                  'e_convergence': 1e-8,
                  'd_convergence': 1e-8})

sapt = helper_SAPT(dimer, memory=8)<jupyter_output>Initializing SAPT object...

RHF for monomer A finished in 0.86 seconds.
RHF for monomer B finished in 0.86 seconds.
Building ERI tensor...
...built ERI tensor in 2.015 seconds.
Size of the ERI tensor is 0.08 GB, 56 basis functions.

...finished initializing SAPT object in  3.90 seconds.<jupyter_text>Before we start computing the SAPT0 corrections, we still need to specify the pertinent notation and define the matrix elements that we will be requesting from `helper_SAPT`. In the classic SAPT papers [Rybak:1991], orbital indices $a,a',a'',\ldots$ and $b,b',b'',\ldots$ denote occupied orbitals of monomers A and B, respectively. The virtual orbitals of monomers A and B are denoted by $r,r',r'',\ldots$ and $s,s',s'',\ldots$, respectively. The overlap integral $S^x_y=\langle x|\rangle y$ reduces to a Kronecker delta when two orbitals from the same monomer are involved, for example, $S^a_{a'}=\delta_{aa'}$, $S^a_r=0$, however, the intermolecular overlap integrals cannot be simplified in any general fashion. Any kind of overlap integral can be requested by calling `sapt.s`, for example, `sapt.s('ab')` gives the $S^a_b$ matrix. For the convenience of implementation, the one-electron (nuclear attraction) $(v_{\rm X})^x_y$ (X = A or B) and nuclear repulsion $V_0$ contributions are usually folded into the two-electron integrals $v^{xy}_{zw}\equiv (xz|yw)$ forming the *dressed* integrals $\tilde{v}$:\begin{equation}\tilde{v}^{xy}_{zw}=v^{xy}_{zw}+(v_{\rm A})^{y}_{w}S^{x}_{z}/N_{\rm A}+(v_{\rm B})^{x}_{z}S^{y}_{w}/N_{\rm B}+V_0S^{x}_{z}S^{y}_{w}/N_{\rm A}N_{\rm B},\end{equation}where $N_{\rm X}$, X=A,B, is the number of electrons in monomer X. An arbitrary *dressed* integral $\tilde{v}^{xy}_{zw}$ can be requested by calling `sapt.vt('xyzw')`. Finally, the HF orbital energy for either monomer can be obtained by calling `sapt.eps`; for example, `sapt.eps('r')` returns a 1D array of virtual orbital energies for monomer A.  2. Electrostatic energyThe SAPT0 electrostatic energy $E^{(100)}_{\rm elst}$ is simply the expectation value of the intermolecular interaction operator $V$ over the zeroth-order wavefunction which is the product of HF determinants for monomers A and B. For the interaction of two closed-shell systems, this energy is obtained by a simple summation of *dressed* two-electron integrals over occupied orbitals of A and B:\begin{equation}E^{(100)}_{\rm elst}=4\tilde{v}^{ab}_{ab}.\end{equation}<jupyter_code>### Start E100 Electrostatics
elst_timer = sapt_timer('electrostatics')
Elst10 = 4 * np.einsum('abab', sapt.vt('abab'), optimize=True)
elst_timer.stop()
### End E100 Electrostatics<jupyter_output>Starting electrostatics...
...electrostatics took a total of  0.09 seconds.<jupyter_text>3. First-order exchange energyThe SAPT0 first-order exchange energy $E^{(100)}_{\rm exch}$ within the $S^2$ approximation and the density matrix formalism is given by Eq. (40) of [Moszynski:1994b]:\begin{align}E^{(100)}_{\rm exch}=&-2\left[\tilde{v}^{ba}_{ab}+S^b_{a'}\left(2\tilde{v}^{aa'}_{ab}-\tilde{v}^{a'a}_{ab}\right)+S^{a}_{b'}\left(2\tilde{v}^{b'b}_{ab}-\tilde{v}^{bb'}_{ab}\right)\right.\\ &\left.-2S^b_{a'}S^{a'}_{b'}\tilde{v}^{ab'}_{ab}-2S^{b'}_{a'}S^{a}_{b'}\tilde{v}^{a'b}_{ab}+S^b_{a'}S^{a}_{b'}\tilde{v}^{a'b'}_{ab}\right]\end{align}and involves several different types of *dressed* MO integrals as well as some intermolecular overlap integrals (not that all indices still pertain to occupied orbitals in this formalism). In Psi4NumPy, each tensor contraction in the above expression can be performed with a single `np.einsum` call:<jupyter_code>### Start E100 Exchange
exch_timer = sapt_timer('exchange')
vt_abba = sapt.vt('abba')
vt_abaa = sapt.vt('abaa')
vt_abbb = sapt.vt('abbb')
vt_abab = sapt.vt('abab')
s_ab = sapt.s('ab')

Exch100 = np.einsum('abba', vt_abba, optimize=True)

tmp = 2 * vt_abaa - vt_abaa.swapaxes(2, 3)
Exch100 += np.einsum('Ab,abaA', s_ab, tmp, optimize=True)

tmp = 2 * vt_abbb - vt_abbb.swapaxes(2, 3)
Exch100 += np.einsum('Ba,abBb', s_ab.T, tmp, optimize=True)

Exch100 -= 2 * np.einsum('Ab,BA,abaB', s_ab, s_ab.T, vt_abab, optimize=True)
Exch100 -= 2 * np.einsum('AB,Ba,abAb', s_ab, s_ab.T, vt_abab, optimize=True)
Exch100 += np.einsum('Ab,Ba,abAB', s_ab, s_ab.T, vt_abab, optimize=True)

Exch100 *= -2
exch_timer.stop()
### End E100 (S^2) Exchange<jupyter_output>Starting exchange...
...exchange took a total of  0.32 seconds.<jupyter_text>4. Dispersion energyThe SAPT0 dispersion energy $E^{(200)}_{\rm disp}$ is given by the formula\begin{equation}E^{(200)}_{\rm disp}=4t^{rs}_{ab}v^{ab}_{rs}\end{equation}where the *dispersion amplitude* $t^{rs}_{ab}$, representing a single excitation on A and a single excitation on B, involves a two-electron integral and an excitation energy denominator:\begin{equation}t^{rs}_{ab}=\frac{v_{ab}^{rs}}{\epsilon_a+\epsilon_b-\epsilon_r-\epsilon_s}\end{equation}Note that for this particular type of integral $\tilde{v}^{ab}_{rs}=v^{ab}_{rs}$: therefore, `sapt.v` instead of `sapt.vt` is used to prepare this tensor.<jupyter_code>### Start E200 Disp
disp_timer = sapt_timer('dispersion')
v_abrs = sapt.v('abrs')
v_rsab = sapt.v('rsab')
e_rsab = 1/(-sapt.eps('r', dim=4) - sapt.eps('s', dim=3) + sapt.eps('a', dim=2) + sapt.eps('b'))

Disp200 = 4 * np.einsum('rsab,rsab,abrs->', e_rsab, v_rsab, v_abrs, optimize=True)
### End E200 Disp<jupyter_output>Starting dispersion...<jupyter_text>5. Exchange dispersion energySome of the formulas for the SAPT0 exchange-dispersion energy $E^{(200)}_{\rm exch-disp}$ in the original papers contained errors. The corrected formula for this term is given by e.g. Eq. (10) of [Patkowski:2007]:\begin{align}E^{(200)}_{\rm exch-disp}=&-2t^{ab}_{rs}\left[\tilde{v}^{sr}_{ab}+S^s_a (2\tilde{v}^{a'r}_{a'b}-\tilde{v}^{ra'}_{a'b})+ S^s_{a'} (2\tilde{v}^{ra'}_{ab}-\tilde{v}^{a'r}_{ab})\right.\\ &+ S^r_b (2\tilde{v}^{sb'}_{ab'}-\tilde{v}^{b's}_{ab'})+ S^r_{b'} (2\tilde{v}^{b's}_{ab}-\tilde{v}^{sb'}_{ab}) \\ &+S^{r}_{b} S^{b'}_{a'} \tilde{v}^{a's}_{ab'}-2 S^{r}_{b'} S^{b'}_{a'} \tilde{v}^{a's}_{ab}-2 S^{r}_{b} S^{b'}_{a} \tilde{v}^{a's}_{a'b'}+4 S^{r}_{b'} S^{b'}_{a} \tilde{v}^{a's}_{a'b} \\ &-2 S^{s}_{a} S^{a'}_{b} \tilde{v}^{rb'}_{a'b'}+4 S^{s}_{a'} S^{a'}_{b} \tilde{v}^{rb'}_{ab'}+ S^{s}_{a} S^{a'}_{b'} \tilde{v}^{rb'}_{a'b}-2 S^{s}_{a'} S^{a'}_{b'} \tilde{v}^{rb'}_{ab} \\ &+ S^{r}_{b'} S^{s}_{a'} \tilde{v}^{a'b'}_{ab}-2 S^{r}_{b} S^{s}_{a'} \tilde{v}^{a'b'}_{ab'}-2 S^{r}_{b'} S^{s}_{a} \tilde{v}^{a'b'}_{a'b} \\ &\left. + S^{a'}_{b} S^{b'}_{a} \tilde{v}^{rs}_{a'b'}-2 S^{a'}_{b} S^{b'}_{a'} \tilde{v}^{rs}_{ab'}-2 S^{a'}_{b'} S^{b'}_{a} \tilde{v}^{rs}_{a'b}\right]\end{align}The corresponding Psi4NumPy code first recreates the dispersion amplitudes $t^{rs}_{ab}$ and then prepares the tensor `xd_absr` that is equal to the entire expression in brackets. The additional two intermediates `h_abrs` and `q_abrs` collect terms involving one and two overlap integrals, respectively.<jupyter_code>### Start E200 Exchange-Dispersion

# Build t_rsab
t_rsab = np.einsum('rsab,rsab->rsab', v_rsab, e_rsab, optimize=True)

# Build h_abrs
vt_abar = sapt.vt('abar')
vt_abra = sapt.vt('abra')
vt_absb = sapt.vt('absb')
vt_abbs = sapt.vt('abbs')

tmp = 2 * vt_abar - vt_abra.swapaxes(2, 3)
h_abrs = np.einsum('as,AbAr->abrs', sapt.s('as'), tmp, optimize=True)

tmp = 2 * vt_abra - vt_abar.swapaxes(2, 3)
h_abrs += np.einsum('As,abrA->abrs', sapt.s('as'), tmp, optimize=True)

tmp = 2 * vt_absb - vt_abbs.swapaxes(2, 3)
h_abrs += np.einsum('br,aBsB->abrs', sapt.s('br'), tmp, optimize=True)

tmp = 2 * vt_abbs - vt_absb.swapaxes(2, 3)
h_abrs += np.einsum('Br,abBs->abrs', sapt.s('br'), tmp, optimize=True)

# Build q_abrs
vt_abas = sapt.vt('abas')
q_abrs =      np.einsum('br,AB,aBAs->abrs', sapt.s('br'), sapt.s('ab'), vt_abas, optimize=True)
q_abrs -= 2 * np.einsum('Br,AB,abAs->abrs', sapt.s('br'), sapt.s('ab'), vt_abas, optimize=True)
q_abrs -= 2 * np.einsum('br,aB,ABAs->abrs', sapt.s('br'), sapt.s('ab'), vt_abas, optimize=True)
q_abrs += 4 * np.einsum('Br,aB,AbAs->abrs', sapt.s('br'), sapt.s('ab'), vt_abas, optimize=True)

vt_abrb = sapt.vt('abrb')
q_abrs -= 2 * np.einsum('as,bA,ABrB->abrs', sapt.s('as'), sapt.s('ba'), vt_abrb, optimize=True)
q_abrs += 4 * np.einsum('As,bA,aBrB->abrs', sapt.s('as'), sapt.s('ba'), vt_abrb, optimize=True)
q_abrs +=     np.einsum('as,BA,AbrB->abrs', sapt.s('as'), sapt.s('ba'), vt_abrb, optimize=True)
q_abrs -= 2 * np.einsum('As,BA,abrB->abrs', sapt.s('as'), sapt.s('ba'), vt_abrb, optimize=True)

vt_abab = sapt.vt('abab')
q_abrs +=     np.einsum('Br,As,abAB->abrs', sapt.s('br'), sapt.s('as'), vt_abab, optimize=True)
q_abrs -= 2 * np.einsum('br,As,aBAB->abrs', sapt.s('br'), sapt.s('as'), vt_abab, optimize=True)
q_abrs -= 2 * np.einsum('Br,as,AbAB->abrs', sapt.s('br'), sapt.s('as'), vt_abab, optimize=True)

vt_abrs = sapt.vt('abrs')
q_abrs +=     np.einsum('bA,aB,ABrs->abrs', sapt.s('ba'), sapt.s('ab'), vt_abrs, optimize=True)
q_abrs -= 2 * np.einsum('bA,AB,aBrs->abrs', sapt.s('ba'), sapt.s('ab'), vt_abrs, optimize=True)
q_abrs -= 2 * np.einsum('BA,aB,Abrs->abrs', sapt.s('ba'), sapt.s('ab'), vt_abrs, optimize=True)

# Sum it all together
xd_absr = sapt.vt('absr')
xd_absr += h_abrs.swapaxes(2, 3)
xd_absr += q_abrs.swapaxes(2, 3)
ExchDisp20 = -2 * np.einsum('absr,rsab->', xd_absr, t_rsab, optimize=True)

disp_timer.stop()
### End E200 Exchange-Dispersion<jupyter_output>...dispersion took a total of  2.86 seconds.<jupyter_text>6. CPHF coefficients and induction energyAs already mentioned, the induction and exchange-induction contributions to SAPT0 are calculated including the relaxation of one molecule's HF orbitals in the electrostatic potential generated by the other molecule. Mathematically, this relaxation is taken into account by computing the CPHF coefficients $C^a_r$ for monomer A [Caves:1969] that specify the linear response of the HF orbitals of A to the electrostatic potential $\omega_{\rm B}$ generated by the nuclei and electrons of the (unperturbed) monomer B and the analogous coefficients $C^b_s$ that describe the response of B to the electrostatic potential of A. The CPHF coefficients are computed by solving the system of equations\begin{equation}(\epsilon_r-\epsilon_a)C^a_r+(2v^{ar'}_{ra'}-v^{r'a}_{ra'})C^{a'}_{r'}+(2v^{aa'}_{rr'}-v^{a'a}_{rr'})C^{r'}_{a'}=-2\tilde{v}^{ab}_{rb}. \end{equation}and similarly for monomer B. Once the CPHF coefficients are ready, the SAPT0 induction energy $E^{(200)}_{\rm ind,resp}$ can be computed very easily:\begin{equation}E^{(200)}_{\rm ind,resp}=4\tilde{v}^{rb}_{ab}C^a_r+4\tilde{v}^{as}_{ab}C^b_s\end{equation}The call to the `helper_SAPT` function `sapt.chf` generates the corresponding contribution to $E^{(200)}_{\rm ind,resp}$ as a byproduct of the calculation of the CPHF coefficients $C^a_r$/$C^b_s$.<jupyter_code>### Start E200 Induction and Exchange-Induction

# E200Induction and CPHF orbitals
ind_timer = sapt_timer('induction')

CPHF_ra, Ind20_ba = sapt.chf('B', ind=True)
sapt_printer('Ind20,r (A<-B)', Ind20_ba)

CPHF_sb, Ind20_ab = sapt.chf('A', ind=True)
sapt_printer('Ind20,r (A->B)', Ind20_ab)

Ind20r = Ind20_ba + Ind20_ab<jupyter_output>Starting induction...
Ind20,r (A<-B)           -1.39782430 mH       -0.87714733 kcal/mol
Ind20,r (A->B)           -2.86985793 mH       -1.80086168 kcal/mol<jupyter_text>7. Exchange induction energyJust like for induction energy, the SAPT0 exchange-induction energy $E^{(200)}_{\rm exch-ind,resp}$ decomposes into two parts describing the exchange quenching of the polarization of A by B and of the polarization of B by A:\begin{equation}E^{(200)}_{\rm exch-ind,resp}=E^{(200)}_{\rm exch-ind,resp}({\rm A}\leftarrow{\rm B})+E^{(200)}_{\rm exch-ind,resp}({\rm B}\leftarrow{\rm A})\end{equation}Now, the formula for the A$\leftarrow$B part is given e.g. by Eq. (5) of [Patkowski:2007]:\begin{align}E^{(200)}_{\rm exch-ind,resp}({\rm A}\leftarrow {\rm B})=&-2 C^a_r \left[\tilde{v}^{br}_{ab}+2S^b_a\tilde{v}^{a'r}_{a'b}+2S^b_{a'}\tilde{v}^{ra'}_{ab}-S^b_a\tilde{v}^{ra'}_{a'b}-S^b_{a'}\tilde{v}^{a'r}_{ab}+2S^r_{b'}\tilde{v}^{b'b}_{ab}\right.\\ &-S^r_{b'}\tilde{v}^{bb'}_{ab}-2S^b_a S^r_{b'}\tilde{v}^{a'b'}_{a'b}-2S^b_{a'}S^{a'}_{b'}\tilde{v}^{rb'}_{ab}-2S^{b'}_{a'}S^r_{b'}\tilde{v}^{a'b}_{ab}-2S^{b'}_a S^{a'}_{b'}\tilde{v}^{rb}_{a'b}\\ & \left.+S^b_{a'}S^r_{b'}\tilde{v}^{a'b'}_{ab}+S^b_a S^{a'}_{b'}\tilde{v}^{rb'}_{a'b}\right]\end{align}and the corresponding formula for the B$\leftarrow$A part is obtained by interchanging the symbols pertaining to A with those of B $(a\leftrightarrow b,r\leftrightarrow s)$ in the above expression. In this example, the CPHF coefficients $C^a_r$ and $C^b_s$ obtained in the previous section are combined with *dressed* two-electron integrals and overlap integrals to compute the $E^{(200)}_{\rm exch-ind,resp}$ expression term by term.<jupyter_code># Exchange-Induction

# A <- B
vt_abra = sapt.vt('abra')
vt_abar = sapt.vt('abar')
ExchInd20_ab  =     np.einsum('ra,abbr', CPHF_ra, sapt.vt('abbr'), optimize=True)
ExchInd20_ab += 2 * np.einsum('rA,Ab,abar', CPHF_ra, sapt.s('ab'), vt_abar, optimize=True)
ExchInd20_ab += 2 * np.einsum('ra,Ab,abrA', CPHF_ra, sapt.s('ab'), vt_abra, optimize=True)
ExchInd20_ab -=     np.einsum('rA,Ab,abra', CPHF_ra, sapt.s('ab'), vt_abra, optimize=True)

vt_abbb = sapt.vt('abbb')
vt_abab = sapt.vt('abab')
ExchInd20_ab -=     np.einsum('ra,Ab,abAr', CPHF_ra, sapt.s('ab'), vt_abar, optimize=True)
ExchInd20_ab += 2 * np.einsum('ra,Br,abBb', CPHF_ra, sapt.s('br'), vt_abbb, optimize=True)
ExchInd20_ab -=     np.einsum('ra,Br,abbB', CPHF_ra, sapt.s('br'), vt_abbb, optimize=True)
ExchInd20_ab -= 2 * np.einsum('rA,Ab,Br,abaB', CPHF_ra, sapt.s('ab'), sapt.s('br'), vt_abab, optimize=True)

vt_abrb = sapt.vt('abrb')
ExchInd20_ab -= 2 * np.einsum('ra,Ab,BA,abrB', CPHF_ra, sapt.s('ab'), sapt.s('ba'), vt_abrb, optimize=True)
ExchInd20_ab -= 2 * np.einsum('ra,AB,Br,abAb', CPHF_ra, sapt.s('ab'), sapt.s('br'), vt_abab, optimize=True)
ExchInd20_ab -= 2 * np.einsum('rA,AB,Ba,abrb', CPHF_ra, sapt.s('ab'), sapt.s('ba'), vt_abrb, optimize=True)

ExchInd20_ab +=     np.einsum('ra,Ab,Br,abAB', CPHF_ra, sapt.s('ab'), sapt.s('br'), vt_abab, optimize=True)
ExchInd20_ab +=     np.einsum('rA,Ab,Ba,abrB', CPHF_ra, sapt.s('ab'), sapt.s('ba'), vt_abrb, optimize=True)

ExchInd20_ab *= -2
sapt_printer('Exch-Ind20,r (A<-B)', ExchInd20_ab)

# B <- A
vt_abbs = sapt.vt('abbs')
vt_absb = sapt.vt('absb')
ExchInd20_ba  =     np.einsum('sb,absa', CPHF_sb, sapt.vt('absa'), optimize=True)
ExchInd20_ba += 2 * np.einsum('sB,Ba,absb', CPHF_sb, sapt.s('ba'), vt_absb, optimize=True)
ExchInd20_ba += 2 * np.einsum('sb,Ba,abBs', CPHF_sb, sapt.s('ba'), vt_abbs, optimize=True)
ExchInd20_ba -=     np.einsum('sB,Ba,abbs', CPHF_sb, sapt.s('ba'), vt_abbs, optimize=True)

vt_abaa = sapt.vt('abaa')
vt_abab = sapt.vt('abab')
ExchInd20_ba -=     np.einsum('sb,Ba,absB', CPHF_sb, sapt.s('ba'), vt_absb, optimize=True)
ExchInd20_ba += 2 * np.einsum('sb,As,abaA', CPHF_sb, sapt.s('as'), vt_abaa, optimize=True)
ExchInd20_ba -=     np.einsum('sb,As,abAa', CPHF_sb, sapt.s('as'), vt_abaa, optimize=True)
ExchInd20_ba -= 2 * np.einsum('sB,Ba,As,abAb', CPHF_sb, sapt.s('ba'), sapt.s('as'), vt_abab, optimize=True)

vt_abas = sapt.vt('abas')
ExchInd20_ba -= 2 * np.einsum('sb,Ba,AB,abAs', CPHF_sb, sapt.s('ba'), sapt.s('ab'), vt_abas, optimize=True)
ExchInd20_ba -= 2 * np.einsum('sb,BA,As,abaB', CPHF_sb, sapt.s('ba'), sapt.s('as'), vt_abab, optimize=True)
ExchInd20_ba -= 2 * np.einsum('sB,BA,Ab,abas', CPHF_sb, sapt.s('ba'), sapt.s('ab'), vt_abas, optimize=True)

ExchInd20_ba +=     np.einsum('sb,Ba,As,abAB', CPHF_sb, sapt.s('ba'), sapt.s('as'), vt_abab, optimize=True)
ExchInd20_ba +=     np.einsum('sB,Ba,Ab,abAs', CPHF_sb, sapt.s('ba'), sapt.s('ab'), vt_abas, optimize=True)

ExchInd20_ba *= -2
sapt_printer('Exch-Ind20,r (A->B)', ExchInd20_ba)
ExchInd20r = ExchInd20_ba + ExchInd20_ab

ind_timer.stop()
### End E200 Induction and Exchange-Induction<jupyter_output>Exch-Ind20,r (A<-B)       0.88275409 mH        0.55393614 kcal/mol
Exch-Ind20,r (A->B)       1.45460411 mH        0.91277717 kcal/mol
...induction took a total of  5.56 seconds.<jupyter_text>8. Summary tableAll the SAPT0 interaction energy contributions have been calculated. All that is left to do is to print out the contributions and the total energy, and to compare the results with the SAPT0 corrections calculated directly by Psi4.<jupyter_code>print('SAPT0 Results')
print('-' * 70)
sapt_printer('Exch10 (S^2)', Exch100)
sapt_printer('Elst10', Elst10)
sapt_printer('Disp20', Disp200)
sapt_printer('Exch-Disp20', ExchDisp20)
sapt_printer('Ind20,r', Ind20r)
sapt_printer('Exch-Ind20,r', ExchInd20r)

print('-' * 70)
sapt0 = Exch100 + Elst10 + Disp200 + ExchDisp20 + Ind20r + ExchInd20r
sapt_printer('Total SAPT0', sapt0)

# ==> Compare to Psi4 <==
psi4.set_options({'df_basis_sapt':'aug-cc-pvtz-ri'})
psi4.energy('sapt0')
Eelst = psi4.variable('SAPT ELST ENERGY')
Eexch = psi4.variable('SAPT EXCH10(S^2) ENERGY')
Eind  = psi4.variable('SAPT IND20,R ENERGY')
Eexind  = psi4.variable('SAPT EXCH-IND20,R ENERGY')
Edisp  = psi4.variable('SAPT DISP20 ENERGY')
Eexdisp  = psi4.variable('SAPT EXCH-DISP20 ENERGY')
psi4.compare_values(Eelst, Elst10, 6, 'Elst100')
psi4.compare_values(Eexch, Exch100, 6, 'Exch100(S^2)')
psi4.compare_values(Edisp, Disp200, 6, 'Disp200')
psi4.compare_values(Eexdisp, ExchDisp20, 6, 'Exch-Disp200')
psi4.compare_values(Eind, Ind20r, 6, 'Ind200,r')
psi4.compare_values(Eexind, ExchInd20r, 6, 'Exch-Ind200,r')<jupyter_output>SAPT0 Results
----------------------------------------------------------------------
Exch10 (S^2)             10.67144721 mH        6.69642917 kcal/mol
Elst10                  -13.76363621 mH       -8.63680559 kcal/mol
Disp20                   -2.49082954 mH       -1.56301796 kcal/mol
Exch-Disp20               0.52133730 mH        0.32714385 kcal/mol
Ind20,r                  -4.26768223 mH       -2.67800901 kcal/mol
Exch-Ind20,r              2.33735821 mH        1.46671331 kcal/mol
----------------------------------------------------------------------
Total SAPT0              -6.99200527 mH       -4.38754623 kcal/mol
    Elst100...........................................................PASSED
    Exch100(S^2)......................................................PASSED
    Disp200...........................................................PASSED
    Exch-Disp200......................................................PASSED
    Ind200,r.......................................................[...]
{'fraction_non_alphanumeric': 0.17942771794670231, 'fraction_numerical': 0.03934439116421665, 'mean_word_length': 4.842250287026406, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '237396', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 10197, 'n_tokens_neox': 9177, 'n_words': 2527, 'path': 'Example/Psi4Numpy/07-SymmetryAdaptedPerturbationTheory/7a_sapt0_mo.ipynb', 'repo_name': 'yychuang/109-2-compchem-lite', 'size': '25454'}