††thanks: These two authors contributed equally to this work.††thanks: These
two authors contributed equally to this work.

Present address: ]Google Quantum AI, Mountain View, CA, 94043

# Interaction Energies on Noisy Intermediate-Scale Quantum Computers

Matthias Loipersberger  0000-0002-3648-0101 QC Ware Corporation, Palo Alto,
CA, 94301, USA  Fionn D. Malone  0000-0001-9239-0162 [  QC Ware Corporation,
Palo Alto, CA, 94301, USA  Alicia R. Welden  0000-0002-2238-9825 QC Ware
Corporation, Palo Alto, CA, 94301, USA  Robert M. Parrish  0000-0002-2406-4741
<EMAIL_ADDRESS>QC Ware Corporation, Palo Alto, CA, 94301, USA  Thomas
Fox  0000-0002-1054-4701 Medicinal Chemistry, Boehringer Ingelheim Pharma GmbH
& Co. KG, Birkendorfer Straße 65, 88397 Biberach an der Riß, Germany  Matthias
Degroote  0000-0002-8850-7708 Elica Kyoseva  0000-0002-9154-0293 Nikolaj Moll
0000-0001-5645-4667<EMAIL_ADDRESS>Raffaele Santagati
0000-0001-9645-0580 Michael Streif  0000-0002-7509-4748 Quantum Lab,
Boehringer Ingelheim, 55218 Ingelheim am Rhein, Germany

###### Abstract

The computation of interaction energies on noisy intermediate-scale quantum
(NISQ) computers appears to be challenging with straightforward application of
existing quantum algorithms. For example, use of the standard supermolecular
method with the variational quantum eigensolver (VQE) would require extremely
precise resolution of the total energies of the fragments to provide for
accurate subtraction to the interaction energy. Here we present a symmetry-
adapted perturbation theory (SAPT) method that may provide interaction
energies with high quantum resource efficiency. Of particular note, we present
a quantum extended random-phase approximation (ERPA) treatment of the SAPT
second-order induction and dispersion terms, including exchange counterparts.
Together with previous work on first-order terms, this provides a recipe for
complete SAPT(VQE) interaction energies up to second order. The SAPT
interaction energy terms are computed as first-level observables with no
subtraction of monomer energies invoked, and the only quantum observations
needed are the the VQE one- and two-particle density matrices. We find
empirically that SAPT(VQE) can provide accurate interaction energies even with
coarsely optimized, low circuit depth wavefunctions from the quantum computer,
simulated through ideal statevectors. The errors on the total interaction
energy are orders of magnitude lower than the corresponding VQE total energy
errors of the monomer wavefunctions.

Keywords

Quantum Computing, Intermolecular Interactions, Quantum Chemistry, Drug
Discovery, SAPT

##  I Introduction

Quantum computing has emerged as a promising platform to approach classically
challenging problems in chemistry Cao  _et al._ (2018); Bauer  _et al._
(2020). The most interesting near-term application is the simulation of
strongly-correlated systems for which the electronic structure cannot be
described with a single Slater determinant. For such systems Kohn-Sham density
functional theory Hohenberg and Kohn (1964); Kohn and Sham (1965) (KS-DFT) may
fail to describe the electronic structure correctly; popular examples are the
Fe-S clusters or the FeMo-cofactor Reiher  _et al._ (2017); Li  _et al._
(2019).

Classically, a proper treatment of these strongly correlated systems is
achieved with multi-reference methods where the naïve combinatorial scaling of
the wavefunction ansatz limits its applications. Note however that much
progress has been made on classical heuristics for wavefunction methods that
exhibit less than combinatorial scaling and that may be highly accurate for a
broad range of problems White and Martin (1999); Chan and Sharma (2011);
Lüchow (2011); Garniron  _et al._ (2018); Levine  _et al._ (2020); Holmes  _et
al._ (2016); Sharma  _et al._ (2017). Alternatively, quantum algorithms
Aspuru-Guzik  _et al._ (2005); McArdle  _et al._ (2020); Cao  _et al._ (2019);
Bauer  _et al._ (2020) might be used to solve the Schrödinger equation with a
resource cost that scales polynomially with the number of qubits.
Unfortunately, the currently available noisy intermediate-scale quantum (NISQ)
hardware Preskill (2018) suffers from relatively poor gate fidelity and low
qubit count Cohn  _et al._ (2021) which poses two key challenges. First, it is
important for NISQ-tailored quantum algorithms Bharti  _et al._ (2022) to
minimize quantum resources. The most prominent NISQ methods are hybrid
quantum-classical algorithms like the variational quantum eigensolver (VQE)
Peruzzo  _et al._ (2014); McClean  _et al._ (2016), quantum Krylov methods
Parrish and McMahon (2019); Stair  _et al._ (2020); Huggins  _et al._ (2020);
Cohn  _et al._ (2021); Klymko  _et al._ (2022) or the fermionic quantum Monte
Carlo method Huggins  _et al._ (2022). The second challenge is to find
specific applications that could harness quantum computing Elfving  _et al._
(2020); Liu  _et al._ (2022). Many application studies in chemistry use either
reduced model systems or molecules with a simple electronic structure Kirsopp
_et al._ (2021); Rice  _et al._ (2021); Kim  _et al._ (2022); Greene-Diniz
_et al._ (2022). There are several promising application areas for quantum
chemistry in computer aided drug design Heifetz (2020) namely, exploring
potential energy surfaces O’Brien  _et al._ (2021), simulating metalloenzymes
Goings  _et al._ (2022) and computing protein-ligand interaction energies
Malone  _et al._ (2022), the last of which we consider in this work.

The computation of non-covalent interaction energies is a routine task in
classical quantum chemistry Rezac and Hobza (2016); Kodrycka and Patkowski
(2019) and the standard procedure is the supermolecular approach: the
interaction energy is calculated as the difference between the ground state
energies of the dimer and two monomers separated to infinity Boys and Bernardi
(1970a). However, transferring this approach to a NISQ type quantum computer
is difficult for several reasons: First, the VQE total energies (on the order
of thousands of kcal/mol) need to be tightly converged to resolve interaction
energies on the order a few kcal/mol with the supermolecular approach. This is
a disadvantage on NISQ hardware because the total energy expectation value is
obtained statistically and high-precision expectation values require a high
number of measurements. Furthermore, converging the VQE total energy to high
accuracy requires deep circuits associated with a large sets of parameters,
where it becomes increasingly difficult to reach the global minimum on the
parameter surface McClean  _et al._ (2018). Second, accounting for the basis
set superposition error (BSSE) in the supermolecular approach Boys and
Bernardi (1970a) is commonly achieved by expanding the basis in the monomer
calculations to the size of the dimer basis. This unnecessarily increases the
qubit count requirements for the individual monomers and can potentially lead
to convergence issues for VQE Peruzzo  _et al._ (2014); McClean  _et al._
(2016).

To this end, this work provides an alternative pathway towards interaction
energies with high accuracy and low quantum resource requirements by using
symmetry-adapted perturbation theory (SAPT) Jeziorski  _et al._ (1976, 1994);
Patkowski (2020) in combination with VQE monomer wave functions [SAPT(VQE)].
This approach directly computes the interaction energy as a sum of small
expectation values; in contrast, the supermolecular approach computes the
small interaction energy (several kcal/mol) as a difference of large total
energies (thousands of kcal/mol). This work builds on our previous work Malone
_et al._ (2022) where we presented the implementation of the first-order
SAPT(VQE) terms of electrostatics and exchange. However, the first-order
SAPT(VQE) terms alone are not capable of computing accurate interaction
energies - standard levels of SAPT also require the computation of the second-
order induction and dispersion terms Jeziorski  _et al._ (1994). Up to now,
the absence of a SAPT(VQE) recipe for the complete second-order SAPT terms has
been a major potential weakness of the approach - indeed, other authorsKirsopp
_et al._ (2021) have noted that “[first-order SAPT(VQE)] computed interaction
energies did not reproduce ligand rankings yielded by more accurate 2nd order
SAPT calculations, due to the missing induction and dispersion components in
the 1st order approximation,” and conclude that “[the first-order SAPT(VQE)]
workflow is limited by truncation of the SAPT expansion at 1st order and it is
not clear how their method can be effectively extended to higher orders.”

In this work, we ameliorate this crucial deficiency by direct implementation
of the second-order induction and dispersion terms, together with their
exchange counterparts. Our approach follows the SAPT(FCI) approach of Korona
Korona  _et al._ (1997) but with the naively-exponential-scaling FCI piece
replaced by an active-space VQE wavefunction which is intended to be
implemented on a forthcoming NISQ computer (in this work we use ideal
statevector simulators for the numerical tests). To implement the second-order
terms, we follow the extended random phase approximation (ERPA) formalism for
SAPT(CASSCF) (complete active space self-consistent field (CASSCF)) developed
by Hapka et al. Hapka  _et al._ (2021a), with VQE standing in for the FCI
solver in CASSCF. In the approach, the response equations for the coupled
polarization propagators are carried out in a truncated hole-particle basis
reminiscent of the Casida expansion Casida  _et al._ (1998) in TD-DFT or the
quantum subspace expansion (QSE) for VQE excited states McClean  _et al._
(2017). While this treatment necessarily does not include all Hamiltonian
states even in the FCI limit, it does use a set of coupled hole-particle
states that span the full range of energetic scales of the Hilbert space, and
which are empirically known to provide highly-converged results for the
induction and dispersion energies. Notably, the use of the ERPA formalism in
SAPT(VQE) allows for the computation of the second-order induction and
dispersion terms, together with exchange counterparts (here in the $S^{2}$
approximation) with the active-space one- and two-particle density matrices of
the VQE ground state wavefunctions appearing as the only “new” quantum
observables. These observables are polynomial scaling and typically readily
available as a byproduct of the VQE optimization procedure. The subsequent
ERPA equations, induction and dispersion contractions, and exchange
counterparts are polynomial scaling classical operations. They are
significantly more complicated than the first-order terms - roughly 200
equations are needed to describe the implementation (see supplemental material
for full details), and the naive CuPy implementation of the equations
implemented here is restricted to smaller systems than our previous paper due
to classical postprocessing overhead (though optimizations such as
hole/active/particle separations and density fitting might significantly
reduce this overhead). Complexity notwithstanding, the approach presented
provides a recipe for a SAPT(VQE) doppelganger of SAPT(FCI) complete through
all second order terms. SAPT methods complete through second order are well
known to produce interaction energies with high accuracy even with modest
basis sets Korona (2013); Parker  _et al._ (2014).

Below, we first lay out the motivation and high-level theory for the ERPA
treatment of the second-order SAPT(VQE) approach. The ERPA equations,
contractions to induction and dispersion terms, and exchange counterparts are
straightforward but extremely verbose, so much of their explicit presentation
is deferred to the supplemental material. We then demonstrate the numerical
performance of all four SAPT(VQE) terms and total interaction energies for
several small multireference dimers and a model heme-nitrosyl hydrogen bonding
complex, using classical ideal statevector simulators to emulate the VQE.

##  II Theory

The most direct “supermolecular” route to the interaction energy of two
monomers $A$ and $B$ is to simply compute the total energy of the combined
system $E_{AB}$ and subtract the total energies of its non-interacting
constituents $E_{A}+E_{B}$, i.e.,

$E_{\mathrm{int}}=E_{AB}-E_{A}-E_{B},$ (1)

ideally using the exact wave function for each system.

An alternative approach to computing intermolecular interaction energies is
symmetry adapted perturbation theory (SAPT), which obtains the interaction
energy with a different approach. In particular, we can write the Hamiltonian
of the combined system as

$\hat{H}=\hat{H}_{A}+\hat{H}_{B}+\hat{V},$ (2)

where we assume $\hat{H}_{X}\lvert\Psi_{X}\rangle=E_{X}\lvert\Psi_{X}\rangle$,
where $\lvert\Psi_{X}\rangle$ is the ground state wavefunction of monomer $X$
and $\hat{V}$ contains only the Coulombic interactions between monomer $A$ and
$B$. With this partitioning of the Hamiltonian we can build a perturbation
theory for the intermolecular interaction energy directly, thus avoiding
computing potentially very large total energies. More explicitly we have

$E_{\mathrm{int}}=\sum_{n}(E_{\mathrm{pol}}^{(n)}+E_{\mathrm{exch}}^{(n)}),$
(3)

where $E_{\mathrm{pol}}^{(n)}$ and $E_{\mathrm{exch}}^{(n)}$ are $n$th-order
polarization and exchange energies respectively. The combination of both first
and second order terms yields the following expression for the interaction
energy (see SI sections 1 to 3 for a detailed derivation):

$E_{\mathrm{int}}\approx
E_{\mathrm{SAPT}}=E_{\mathrm{elst}}+E_{\mathrm{exch}}+E_{\mathrm{ind,u}}+E_{\mathrm{disp}}$
(4)

where $E_{\mathrm{elst}}$ denotes the electrostatic contribution to the total
interaction energy energy, $E_{\mathrm{exch}}$ the exchange energy,
$E_{\mathrm{ind,u}}$ the induction energy and $E_{\mathrm{disp}}$ the
dispersion energy. If we use high quality monomer wavefunctions approaching
full configuration interaction (FCI), we conceptually approch the SAPT(FCI)
method of Korona and co-workers. In Korona’s original work, the computation of
the second-order induction and dispersion terms is performed by a direct
response property treatment of the static (induction) and frequency dependent
dynamic susceptibility tensors (dispersion), followed by contractions of these
property tensors to form the polarization and exchange SAPT
contributions.Korona  _et al._ (1997); Korona (2008a, b); Korona and Jeziorski
(2008); Korona (2009a) In this work, we instead employ the extended random
phase approximation (ERPA) Chatterjee and Pernal (2012) as pioneered within
SAPT(CASSCF) by Hapka Hapka  _et al._ (2018, 2019, 2021b, 2021a) to avoid
computing excited state properties on the quantum computer as they often
require a significant measurement overhead on NISQ-era quantum computers
Parrish and McMahon (2019); Parrish  _et al._ (2019); Cai  _et al._ (2020);
Ollitrault  _et al._ (2020). We note that the use of ERPA as a proxy for
explicit response properties of VQE resembles the quantum subspace expansion
(QSE) method McClean  _et al._ (2017), wherein a basis of single and double
excitations out of a VQE reference is used to provide an excited state ansatz
that is truncated in character (but not in excitation energy). Further details
of SAPT and the ERPA procedure are provided in the supporting information (SI
sections 1–3). We note that SAPT interaction energies using this second-order
truncation of SAPT are typically highly accurate, even with Hartree-Fock
wavefunctions (for the case of single-reference systems), and well established
to produce accurate interaction energies in many common cases Parker  _et al._
(2014). Third- and higher-order extensions are likely possible along similar
response property or ERPA lines as used here, but are typically found to not
improve the SAPT interaction energy significantly past the second-order level.

The current NISQ-era hardware is limited to tens of qubits (spin-orbitals);
therefore an active space formalism is necessary to describe realistic
chemical systems. In the active space approach, we partition the one-electron
orbital set into core orbitals, active orbitals and virtual orbitals. Ideally,
the active orbitals contain the orbitals required to describe the entangled
electrons properly. The active space of the wave function is then is
calculated on a quantum computer (see SI section 1.2 for more information). In
the SAPT(VQE) approach, one or both of the monomer active space wavefunctions
are generated by VQE-type quantum circuits:

$|\Psi_{\mathrm{VQE}}\rangle\equiv\hat{U}_{\mathrm{VQE}}|\Phi_{I}\rangle$ (5)

where $|\Phi_{I}\rangle$ is some initial state (typically the Hartree–Fock
determinant). From these wavefunctions we obtain single-particle and two-
particle reduced density matrices that go into the computation of the SAPT
interaction energy. In this work we use a modified version of the unitary
cluster Jastrow wavefunction ansatz (Matsuzawa and Kurashige, 2020) (VQE)
which takes the form

$|\Psi_{\mathrm{VQE}}\rangle=\prod_{k}\exp(-\hat{K}^{(k)})\exp(\hat{T}^{(k)})\exp(+\hat{K}^{(k)})|\Phi_{I}\rangle,$
(6)

where $\hat{K}^{(k)}$ and $\hat{T}^{(k)}$ are one- and two-body operators, and
$k$ is a parameter that controls the depth of the circuit and as a result its
variational freedom. We use a slightly modified $k$-uCJ ansatz from Ref. 61,
which we denote as $k$-muCJ for clarity, with the ‘m’ standing for modified
(see SI section 1.3 for more details). The SAPT(VQE) workflow is outlined at a
high level in Fig. 1. We note that the SAPT(VQE) method as formulated within
the ERPA picture is independent of the quantum algorithm used to determine the
density matrices and thus can likely be readily adopted to any quantum
algorithm of choice.

##  III Results

In a first step, we test SAPT(VQE) with two classic intermolecular interaction
motifs (water dimer and t-shaped benzene dimer). In a second step, we apply
SAPT(VQE) to a heme-nitrosyl model complex; these systems are highly relevant
in both biological Maia and Moura (2014) and pharmaceutical Moncada (1991)
chemistry. The SAPT(VQE) results presented in this section are the result of
ideal statevector VQE simulations (see SI section 4 for for more details). In
all examples, we benchmark the accuracy of the VQE/SAPT(VQE) results by
comparing to classical SAPT(CAS-CI) energies using the same orbitals and
active space (see Figs. 2–5). The complete active space configuration
interaction (CAS-CI) wavefunction represents the exact wavefunction within the
active space approximation and thus, SAPT(CAS-CI) results represent the best
possible interaction energy within the SAPT approximation but not the exact
interaction energy (see Fig. S6 for a comparison of VQE and CAS-CI wave
functions). We note that this comparison is only possible for small active
space sizes due to the combinatorial scaling of the CAS-CI wavefunction
ansatz.

###  III.1 Multi-Reference Benchmark Systems

The chemistry of non-covalent interactions governs a wide range of interaction
motifs such as hydrogen bonds or dispersion bound systems. However, the
electronic structure of these simple systems is often well described by
classical single reference methods. Therefore, we modified two of the classic
systems in our previous study Malone  _et al._ (2022), namely the water dimer
and the t-shaped benzene dimer to make the electronic structure strongly
correlated and thus challenging to compute accurately for conventional single
reference methods.

The first test case is a hydrogen bonding motif: the stretched water dimer
complex, which is depicted in Fig. 2 (a). The two partially broken single
bonds make this system strongly correlated and require a multi reference
treatment to accurately describe the electronic structure. We included all
eight valence electrons of the stretched monomer and eight spatial orbitals
(8e, 8o) in the active space (for a detailed procedure on how the active
orbitals were selected for the CASSCF calculations see SI section 1). The
CASSCF natural orbital occupation numbers (NOON, see Fig. S2) exhibit
deviations from integer values, which is an indicator of strong correlation.
Consequently, the single reference RHF method fails to describe this system as
apparent by the large deviation of the absolute energy of the monomer (see
Fig. 2 (b) black dotted line). This system also provides a challenge for the
quantum algorithm as the $k$-muCJ ansatz needs a repetition factor of $k>10$
to converge the absolute energy to the stretched monomer below 1 kcal/mol.

In sharp contrast to this, the errors of the total interaction energy as well
as the errors of each individual SAPT energy term are multiple orders of
magnitude lower than the absolute VQE errors (see Fig. 2 (b)). In fact, the
very shallow $k=1$ circuit is accurate enough to provide interaction energies
in comparison to the SAPT(CAS-CI) results. In a next step, we probed the bond
dissociation of the water complex (along the intermolecular O-H bond labeled
$r$ in Fig. 2 (a)). We find a similar behavior: the errors in the interaction
energies and each energy component are below 1 kcal/mol for all intermolecular
distances and several orders of magnitude lower than the error in absolute
energies. For $k=4$, each interaction energy is below the 1 kcal/mol
threshold; for $k=1$ some errors are slightly larger at small intermolecular
distances.

The second test case is a dispersion bound complex: the T-shaped benzene
p-benzyne dimer, depicted in Fig. S4 (a). The p-benzyne monomer has a
biradical ground state, which is difficult to describe with classical single
reference methods Crawford  _et al._ (2001). The key findings are identical to
the previous test case. Thus, these findings hold for very different types of
intermolecular interaction motifs, intermolecular distances, different active
spaces and different type difficult electronic structures.

###  III.2 Hydrogen Bonding to Heme-Nitrosyl Model Complexes

As an application example, we study hydrogen bonding to a Manganese nitrosyl
complex. Nitric oxide (NO) is a small molecule with important biological
implications such as signal transduction Bruckdorfer (2005); Moncada (1991) or
as a key intermediate in the nitrogen cycle Maia and Moura (2014); Lehnert
_et al._ (2018). At the center of these processes are metalloporpyhrins Kadish
_et al._ (2000); Lehnert  _et al._ (2013); Hunt and Lehnert (2015), where NO
binds to the metal center as a nitrosyl ligand Santolini (2011); Walker
(2005).

In order to understand and control these biological processes, the chemistry
around the Metal-NO bonds must be elucidated in terms of electronic structure
and reactivity Hayton  _et al._ (2002) as illustrated by theoretical Radon and
Pierloot (2008), experimental Lehnert  _et al._ (2021) and medicinal Tfouni
_et al._ (2010); Hickok and Thomas (2010); Serafim  _et al._ (2012) work.
Unfortunately, the metal-NO bond in nitrosyl complexes poses a challenge for
many quantum chemistry methods due to the redox active nature of NO ligand
Ampßler  _et al._ (2020). There are three possible oxidation states for the NO
moiety: NO-, NO∙ and NO+, which is illustrated for a generic M(II)-NO complex
in Figure 4. In many cases, the bond is best described in terms of a
superposition of these states. This strong correlation makes this a
challenging system for many single reference methods such as DFT Boguslawski
_et al._ (2011). This results in a large variety of recommended functionals
depending on the specific nitrosyl complex studied Radon and Pierloot (2008);
Siegbahn  _et al._ (2010); Goodrich and Lehnert (2013); Lehnert  _et al._
(2013).

In this work, we study the hydrogen bonding to a heme-model manganese-nitrosyl
complex. This model can serve as a proxy of how a metal-heme bound NO
interacts with a protein environment as depicted in Fig. 4. The metal-
coordinating cyano and ammonia ligands resemble a porphyrin coordination
environment in terms of ligand field, total charge and $\pi$-acidity (see Fig.
4). As hydrogen bond donors we chose HF, H2O, NH3 and CH4, which cover a wide
range of donor strengths similarly to active sites in proteins. The resulting
four hydrogen complexes are depicted in Fig. 4 and are abbreviated as MnNO…HF,
MnNO…HOH, MnNO…HNH2 and MnNO…HCH3.

As a first step, we analyzed the electronic structure of [Mn(CN)2(NH3)3NO]^0
and found the system to be strongly correlated (see SI section 5 for a more
detailed discussion). We used a (6e, 6o) active space for subsequent CASSCF
and VQE calculations to include all 6 electrons of the Mn-NO bond. The key six
active orbitals are centered around the Mn-NO moiety and are very similar to
the active orbitals in real heme nitrosyl complexes Radon  _et al._ (2010)
(see Fig. S2 and S5). The VQE as the $k$-uCJ ansatz required up to $k=8$
layers to converge to the CASSCF energy within 1 kcal/mol (the variations are
caused by the additional ghost-basis functions from the different hydrogen
bonding donors). In contrast, the SAPT interaction energy is already
significantly below that threshold even for $k=1$. It is noteworthy that this
finding holds true for the wide range of interaction energies in this series
($-$0.3 to $-$9.3 kcal/mol, see Fig. 5 (a)). Thus, we confirm the core finding
of this work for different interaction motifs, a wide range of interaction
energies, different strongly correlated electronic structure systems and both
equilibrium and non-equilibrium bond distances.

In addition to the interaction energy, SAPT provides a decomposition into
physical meaningful terms helping to unravel the origin of the interaction.
Fig. 6 (b) plots each component of the SAPT(VQE) ($k=4$) calculation of the
series of hydrogen bonded complexes plus the stretched water dimer as a
reference of a typical hydrogen bond (see see SI section 5.1, Table S1 and
Fig. S7 for more details including bond distances). The main driving force for
binding is the electrostatic term as expected for hydrogen bonds. The
strongest contrast between the water dimer and the Mn-NO hydrogen bonds is
observed in the exchange energy and is the main driving force for the
difference in interaction energies. This may be rationalized by difference in
the diffuseness of the lone pairs: the bound NO becomes (partly) NO^+, which
makes the lone pair more compact in space than the lone pair in the water
dimer, thus resulting in less exchange repulsion.

At last, we compare the SAPT(VQE) interaction energy to DFT based
supermolecular (BSSE corrected Boys and Bernardi (1970b)) interaction
energies, the standard approach on classical hardware. The exact comparison is
difficult as there is no standard procedure to obtain accurate interaction
energies for strongly correlated systems. In addition, we use a small basis
set [due to technical limits in the current CuPy classical implementation of
SAPT(VQE)]. However, nitrosyl complexes are an example of the non-universality
problem of approximate density functionals as the hydrogen bonding moiety and
the nitrosyl moiety prefer different approximate density functionals Radon and
Pierloot (2008); Radon  _et al._ (2010); Mardirossian and Head-Gordon (2014)
and thus reliable predictions are only possible with careful system specific
benchmarking when experimental data is available Lehnert  _et al._ (2013). In
contrast, SAPT is expected to robustly give accurate results for hydrogen
bonds given proper monomer wavefunctions, e.g. via a quantum algorithm in
SAPT(VQE). To illustrate this point, Fig. 6 plots the interaction energies of
SAPT(VQE) , SAPT(CAS-CI) and several popular DFT functionals (through the
supermolecular approach). We included many popular functionals as well as
several top performing functionals for non-covalent interactions Mardirossian
and Head-Gordon (2017). We see in Fig. 6 that the SAPT(VQE) ($k=4$) is almost
identical to the SAPT(CAS-CI) in all four cases. The DFT functionals exhibit a
significant spread for each complex. The B97-D functional predicted the
smallest interaction energy in all four cases, but the highest interaction
energy is predicted in each case by a different functional. Furthermore, we
see the relative ordering of the functionals change for each system (color
sequence in each plot). This illustrates the non-universality problem for
approximate exchange correlation functionals even for very similar nitrosyl
complexes (this also holds true for larger basis set as illustrated in Fig.
S8). Note that the SAPT(CAS-CI) results are the reference for the SAPT(VQE)
calculations and do _not_ represent the true interaction energy, thus, only
the SAPT(VQE) and not the DFT interaction energies should be compared against
this reference.

To demonstrate that the erroneous behavior of DFT is related with the system
studied here, we calculate nitric acid hydrogen complexes with HF, H2O (see
Fig. S9). These hydrogen complexes are the main group analogues of the
nitrosyl complex where we replace the Mn-NO with a H-NO bond. This results in
a much simpler electronic structure without strong correlation where we can
generate reference energies using coupled cluster wave function methods. We
find that many DFT functionals perform within 1 kcal/mol accuracy.
Interestingly, we see that the relative ordering of the functionals changes
notably from the MnNO to the HNO systems. Furthermore, we note that SAPT
predicts the interaction energies with $<$0.5 kcal/mol error in both cases
using the optimal basis set (see Fig. S9 (a) and (b); see SI S5.2 for details
on the reference energies). We can expect SAPT(VQE) ansätze to yield similarly
accurate results for strongly correlated examples as presented above when the
optimal basis is used. Therefore, the SAPT formalism, presented in this work,
is able to provide accurate interaction energies both for simple and difficult
electronic structures, while the accuracy of DFT deteriorates for the latter.

##  IV Conclusion

With the developments of the present manuscript, we have what we believe
represents a minimally complete path to accurate determination of
intermolecular interaction energies on a NISQ-type computer. Our previous
study Malone  _et al._ (2022) established the theoretical framework of
SAPT(VQE) but was limited to only first order terms of electrostatics and
exchange for a proof-of-concept demonstration. This work obviates this
limitation by including the second order terms of induction and dispersion,
including their exchange counterparts, which results in a level of SAPT well-
established to produce accurate interaction energies with chemical accuracy
Parker  _et al._ (2014). In this hybrid quantum-classical procedure, we obtain
the monomer wavefunctions on the quantum computer via the VQE algorithm and
measure the one- and two-particle reduced density matrices of the monomers
(simulated through ideal statevectors in the present work). On the classical
computer we compute the first and second order SAPT contributions based on the
reduced density matrices from VQE calculation. The direct computation of
excited states for the second order terms is avoided via an extended random
phase approximation (ERPA) formalism.

We find empirically that SAPT(VQE) can provide accurate interaction energies
even with coarsely optimized, low circuit depth wavefunctions from the quantum
computer. The resulting errors of first and second order contributions, in
addition to the the total interaction energies, are orders of magnitude lower
than the corresponding VQE total energies of the monomer wavefunctions. Our
empirical findings are based on the application of the SAPT(VQE) method to
several systems with strongly correlated electronic structures: two classic
intermolecular interaction motifs and several hydrogen bonding complexes of a
heme-nitrosyl model complex, a class of biological highly relevant
metalloenzymes where classic quantum chemistry methods such as DFT struggle to
obtain accurate interaction energies. Thus, this works paves the way to obtain
accurate interaction energies on a NISQ-era quantum computer with few quantum
resources. It is a first step to alleviate one of the major challenges in
quantum chemistry where in-depth knowledge of both the method and system is
required a priori to reliably generate accurate interaction energies.

While a basic path to NISQ-type computations of interaction energies is now
reasonably clear, much remains to improve the details of the operational
concept. One basic direction that needs improvement is the classical
acceleration of the ERPA and SAPT terms - our naïve CuPy code for this
consideration was severely limited in system size due to a non-optimal
treatment of core/active/virtual simpliciation and a lack of density fitting
of the response functions. Korona and Jeziorski (2008) In another instance, it
may be that the ERPA formalism for the second order terms could be improved by
direct treatment of the response of the monomers to external perturbations,
including coupled response of the quantum and classical monomer wavefunction
parameters. It might also be the case that more-advanced treatment of the
required observables along the lines of double factorization Motta  _et al._
(2021); Huggins  _et al._ (2021) could significantly reduce the number of
required measurements when the time comes to deploy this method in the
presence of shot noise. Another highly interesting question is how the present
approach might map (if at all) to the fault tolerant regime where statistical
evaluation of expectation values is definitionally prohibitive. Finally,
beyond the concept of interaction energies, this work represents our general
hypothesis of the level of specialization needed to converge various important
chemical observables - we believe that other properties such as gradients,
polarizabilities, spectroscopies, etc may require similar quantum adaption of
rather verbose classical methods like SAPT to provide good convergence of
observables on NISQ devices.

Figure 1: The workflow of SAPT(VQE) to obtain accurate interaction energies
using a NISQ quantum algorithm: first, we define two monomers and the
interaction of interest (exemplary here: two water molecules); second, we
define the active space of the electrons which are simulated on the quantum
computer; third, the quantum computer is used to compute the electronic
structure via a quantum algorithm such as VQE. The converged wavefunction
yields the reduced one- and two-particle density matrices (1- and 2-RDM);
fourth, the classical computer computes the interaction energy as a post
processing step using the RDMs via a sum of electrostatic, exchange, induction
and dispersion energies. The former two terms were published by some of us
previously Malone  _et al._ (2022), the latter two are presented in this study
and emplpy an extended random phase approximation (ERPA) formalism. They
require solving response equations and are necessary to obtain accurate
interaction energies; fifth, the SAPT energy components and interaction
energies provide an in-depth understanding of intermolecular interaction of
interest (note that in this work the quantum computing results were obtained
on a simulator).

(a)

(b)

(c)

Figure 2: (a) Structure of the stretched water dimer; (b) absolute errors of
the VQE total energy (stretched monomer) and each SAPT(VQE) energy term as a
function of the repetition factor $k$ at $r$ = 2.0 Å. This shows that accurate
interaction energies can already be obtained with coarsely optimized, low
circuit depth wavefunctions from the quantum computer (errors relative to the
CAS-CI and SAPT(CAS-CI) energies, the dotted gray line represents the chemical
accuracy 1 kcal/mol threshold); (c) absolute error of the VQE total energy
(stretched water monomer) and each SAPT(VQE) ($k=1$) and SAPT(VQE) ($k=4$)
energy term as a function of the intermolecular distance $r$. This shows that
this empirical finding also holds for non-equilibrium bond distances.

(a)

(b)

(c)

Figure 3: (a) Structure of the benzene p-benzyne dimer; (b) absolute errors of
the VQE total energy (p-benzyne monomer) and each SAPT(VQE) energy term as a
function of the repetition factor $k$ at r = 3.9 Å (error relative to the CAS-
CI and SAPT(CAS-CI) energies, the dotted gray line represents the chemical
accuracy 1 kcal/mol threshold); (c) absolute error of of the VQE total energy
(p-benzyne monomer) and each SAPT(VQE) ($k=1$) and SAPT(VQE) ($k=4$) energy
term as a function of the intermolecular distance r. This shows that our
empirical finding about accurate interaction energies also holds for different
interaction motifs and different strongly correlated electronic structures.
Figure 4: Upper left panel: different oxidation states of a transition metal
bound NO; bottom left panel: structure of the [Mn(CN)2(NH3)3NO]^0 complex and
the natural orbitals (NOON) with the largest deviation from integer value
occupation ($\pi$ type metal to ligand backbonding); right panel:
simplification of the heme ligand framework and schematic representation of
the four manganese nitrosyl hydrogen bond complexes.

(a)

(b)

Figure 5: (a) Top figure: absolute errors of the VQE total energy for the
[Mn(CN)_2(NH3)_3NO]^0 complex and the “stretched” water; bottom figure: the
SAPT(VQE) interaction energy as a function of the repetition factor $k$ for
the water dimer and each model heme-nitrosyl hydrogen bond complex showing
that our empirical finding about accurate interaction energies also holds for
a wide range of interaction energies (error relative to the CAS-CI and
SAPT(CAS-CI) energies, the dotted gray line represents the 1 kcal/mol
threshold); (b) term-by-term decomposition of the SAPT(VQE) ($k=4$)
interaction energies of each model heme-nitrosyl hydrogen complex (see Tab. S1
for details).

Figure 6: Interaction energies using several approximate DFT exchange
correlation functionals (supramolecular approach) in comparison to the
SAPT(CAS-CI) and SAPT(VQE) ($k=4$) interaction energies showing the
sensitivity of the results for a given hydrogen bonding complex - exchange
correlation functional pair (see Table S3 for more details).

## Data availability

XYZ structures and NPZ files for the MnNO…HF example. We also include python
code to obtain and save the active space Hamiltonian via PySCF.

## Acknowledgments

QC Ware Corp. acknowledges generous funding from Boehringer Ingelheim for this
research project. We thank Clemens Utschig-Utschig for insightful discussions
and Daria Miroshko for the graphics design of Fig. 1.

## Author contributions

All authors discussed and designed the formalism of an ERPA-based extension of
SAPT(VQE) for induction and dispersion. F. D. M., A. R. W., and R. M. P.
carried out the initial derivation of the equations. F.D.M. performed
extensive simplifications and spin-adaption of the equations. F. D. M. wrote
the CuPy implementation of the method. M. L., T. F., M. D, and N. M. prepared
the geometries and performed the classical reference computations. M.L. and
F.D.M. performed the simulated SAPT(VQE) computations. All authors assisted
with the production and analysis of the applications results and with the
writing of the manuscript. Correspondence should be addressed to R.M.P. or
N.M.

## Competing interest

M. L., A. R. W. and R. M. P. own stock/options in QC Ware Corp.

##  Appendix A Theory

###  A.1 Indices and Labels

We adopt the following notation for orbital sets used in this work:

  * •

$A/B$ \- monomer nuclear indices.

  * •

$\mu/\nu$ \- nonorthogonal atomic spatial orbital basis indices (i.e.,
Gaussian basis indices).

  * •

$p/q$ \- orthogonal molecular spatial orbital basis indices.

  * •

$i/j$ \- orthogonal occupied spatial orbital basis indices.

  * •

$t/u$ \- orthogonal active spatial orbital basis indices.

  * •

$a/b$ \- orthogonal virtual spatial orbital basis indices.

Repeated indices within a monomer will be denoted with primes, e.g.,
$p,p^{\prime},p^{\prime\prime},p^{\prime\prime\prime}$. When dealing with
spin-orbital quantities, we use the context specific notation of an “unbarred”
orbital index to denote $\alpha$ and a “barred” orbital index to denote
$\beta$, i.e., $p^{\dagger}$ is an $\alpha$ spin-orbital creation operator on
spatial orbital index $p$, while $\bar{p}^{\dagger}$ is a $\beta$ spin-orbital
creation operator on spatial orbital index $p$.

###  A.2 Symmetry Adapted Perturbation Theory

The interaction energy between two monomers is defined as the the total energy
of the combined system and subtract the total energies:

$E_{\mathrm{int}}=E_{AB}-E_{A}-E_{B},$ (7)

where the total energies represent the full CI solution at the basis set
limit. In practice, approximate methods such as density functional theory or
coupled cluster methods are used to compute accurate enough total energies to
resolve the binding energy accurately.

An alternative approach to computing intermolecular interaction energies is
symmetry adapted perturbation theory (SAPT) which is valid for non-covalent
interactions. Instead of computing total energies, SAPT assumes that the
intermolecular interactions are weak and thus we can compute the interaction
energy via perturbation theory. In particular, we can write the Hamiltonian of
the combined system as

$\hat{H}=\hat{H}_{A}+\hat{H}_{B}+\hat{V},$ (8)

where we assume $\hat{H}_{X}\lvert\Psi_{X}\rangle=E_{X}\lvert\Psi_{X}\rangle$,
where $\lvert\Psi_{X}\rangle$ is the ground state wavefunction of monomer $X$
and $\hat{V}$ contains only the Coulombic interactions between monomer $A$ and
$B$. With this partitioning of the Hamiltonian we can build a perturbation
theory for the intermolecular interaction energy directly thus avoiding
computing potentially very large total energies. More explicitly we have

$E_{\mathrm{int}}=\sum_{n}(E_{\mathrm{pol}}^{(n)}+E_{\mathrm{exch}}^{(n)}),$
(9)

where $E_{\mathrm{pol}}^{(n)}$ and $E_{\mathrm{exch}}^{(n)}$ are $n$th-order
polarization and exchange energies respectively. An added benefit of SAPT is
that we obtain an intuitive breakdown of the interaction energy components,
into electrostatic, induction, dispersion and exchange components which can be
used to provide chemical insight into the binding process. At the lowest order
this gives rise to SAPT0Jeziorski  _et al._ (1994) yielding:

$E_{\mathrm{int}}\approx
E_{SAP0}=E^{(1)}_{\mathrm{elst}}+E^{(1)}_{\mathrm{exch}}+E^{(2)}_{\mathrm{ind,u}}+E^{(2)}_{\mathrm{disp}}+E^{(2)}_{\mathrm{exch-
disp}}+E^{(2)}_{\mathrm{exch-ind,u}},$ (10)

where $E^{(1)}_{\mathrm{{elst}}}$ corresponds to the first order
electrostatics term, $E^{(1)}_{\mathrm{{exch}}}$ to the first order exchange
term, $E^{(1)}_{\mathrm{{ind,u}}}$ to the first order induction term,
$E^{(2)}_{\mathrm{{disp}}}$ to the second order dispersion term,
$E^{(2)}_{\mathrm{exch-disp}}$ to the second order exchange-dispersion term,
$E^{(2)}_{\mathrm{exch-ind,u}}$ to the second order exchange-induction term
(see below for a detailed derivation). In the main manuscript we combine both
dispersion and inductions terms and drop the superscripts yielding:

$E_{\mathrm{int}}\approx
E_{\mathrm{elst}}+E_{\mathrm{exch}}+E_{\mathrm{ind,u}}+E_{\mathrm{disp}}.$
(11)

Ref. 49 shows that this SAPT0 using Hartee-Fock wavefunctions and a medium
size jun-cc-pvdz basis can yield highly accurate results for a broad range of
non-covalent interactions.

As in-depth discussed in our previous workMalone  _et al._ (2022), we use the
density matrix formulation of SAPT Moszynski  _et al._ (1994); Korona (2008a,
b, 2009b) as recently fully implemented for complete active space self
consistent field (CASSCF) wavefunctions by Hapka and et al. Hapka  _et al._
(2021b) This formalism allows for the evaluation of the terms appearing in Eq.
10 using just the ground state one- and two-particle reduced density matrices
of the monomers with additional response terms for the second order terms.
Instead of Hartee-Fock density matrices, we use a quantum computer via a VQE
ansatz to determine an accurate ground state wave function of the system (for
more details on those algorithms see the next sections). The detailed
derivation of the first order terms ($E^{(1)}_{\mathrm{{elst}}}$,
$E^{(1)}_{\mathrm{{exch}}}$) be found in Ref. Malone  _et al._ (2022).

The second order SAPT terms for induction and dispersion energies (as well as
their exchange counterparts) requires the calculation of excited state
properties on a quantum computer (see below for a more detailed derivation).
Although several approaches have been suggested in the literature to compute
excited state properties on NISQ-era quantum computers, they often require a
significant measurement overhead Parrish and McMahon (2019); Parrish  _et al._
(2019); Cai  _et al._ (2020); Ollitrault  _et al._ (2020). To reduce this
burden we employ the extended random phase approximation (ERPA) Chatterjee and
Pernal (2012) which requires only the one- and two-body reduced density
matrices to be evaluated on the quantum computer McClean  _et al._ (2017).
This approximation has previously been shown by others to produce quite
accurate interaction energies when employed in SAPT based on CASSCF
wavefunctions and we will show that this carries over for VQE wavefunctions.

Given that NISQ-era devices are currently limited to tens of qubits (spin-
orbitals) we will use an active space approach analogous to CAS-CI methods. In
the active space approach we partition the one-electron orbital set into
$N_{c}$ core orbitals, $N_{a}$ active orbitals and $N_{i}$ virtual orbitals.
This partitioning gives rise to modified monomer Hamiltonian given by (for
example for monomer $A$)

$\begin{split}\hat{H_{A}}^{\prime}=&\sum_{tt^{\prime}\sigma}\tilde{h}_{tt^{\prime}}a^{\dagger}_{t\sigma}a_{t^{\prime}\sigma}+\\\
&\frac{1}{2}\sum_{\sigma\sigma^{\prime}}\sum_{tt^{\prime}t^{\prime\prime}t^{\prime\prime\prime}}(tt^{\prime\prime}|t^{\prime}t^{\prime\prime\prime})a^{\dagger}_{t\sigma}a^{\dagger}_{t^{\prime}\sigma^{\prime}}a_{t^{\prime\prime\prime}\sigma^{\prime}}a_{t^{\prime\prime}\sigma}\end{split},$
(12)

where the modified one-electron integrals $\tilde{h}_{tt^{\prime}}$ now
include core-active space interactions

$\tilde{h}_{tt^{\prime}}=h_{tt^{\prime}}+\sum_{ii^{\prime}}\left[(tt^{\prime}|ii^{\prime})-\frac{1}{2}(ti^{\prime}|it^{\prime})\right]\gamma_{ii^{\prime}}.$
(13)

The key approximation in the active space approach is that a (small) set of
“active” orbitals and electrons are defined a priori and the FCI expansion is
constrained to that subset of electrons and orbitals. The quality of the CAS-
CI and VQE results depends strongly on the selected active space. Veryazov
_et al._ (2011); Sun and Chan (2014); Sayfutyarova  _et al._ (2017).

###  A.3 Variational Quantum Eigensolver

This section briefly summarizes the VQE ansatz used in this work, it is
identical to the approach taken in our previous work Malone  _et al._ (2022)
and discussed in more detail. The SAPT post processing step only relies on the
reduced density matrices which can in principle be generated with an quantum
algorithm such as other VQE flavorsAnselmetti  _et al._ (2021) or Quantum
Krylov methodsParrish and McMahon (2019).

We generate the active space wave function of the strongly correlated monomer
using the VQE ansatz described below (it is also possible that both monomer
wave function are evaluated on the quantum computer but in practice usually
only one monomer will exhibit strong correlation and thus require a VQE
treatment):

$|\Psi_{\mathrm{VQE}}\rangle\equiv\hat{U}_{\mathrm{VQE}}|\Phi_{I}\rangle$ (14)

where $|\Phi_{I}\rangle$ guess state (typically the Hartree–Fock).

Throught this study we only use real active space wavefunction
$|\Psi_{\mathrm{VQE}}\rangle$ and they will be a definite eigenfunction of the
$\hat{N}_{\alpha}$, $\hat{N}_{\beta}$, and $\hat{S}^{2}$ operators.

The paper is using the following Jordan-Wigner representation:

$p^{\pm}=\bigotimes_{p^{\prime}=0}^{p^{\prime}=p-1}\hat{Z}_{p^{\prime}}(\hat{X}_{p}\mp
i\hat{Y}_{p})/2,$ (15)
$\bar{p}^{\pm}=\bigotimes_{p^{\prime}=0}^{p^{\prime}=N_{\alpha}-1}\hat{Z}_{p^{\prime}}\bigotimes_{\bar{p}^{\prime}=0}^{\bar{p}^{\prime}=\bar{p}-1}\hat{Z}_{\bar{p}^{\prime}}(\hat{X}_{\bar{p}}\mp
i\hat{Y}_{\bar{p}})/2,$ (16)

where $p^{+}=p^{\dagger}$ and $p^{-}=p$ and we order the Jordan-Wigner strings
in $\alpha$-then-$\beta$ order and $\hat{Z},\hat{Y}$ and $\hat{X}$ are the
usual Pauli operators.

In this work we use a modified version of the unitary cluster Jastrow
wavefunction Matsuzawa and Kurashige (2020) ($k$-uCJ) which takes the form

$|\Psi_{\mathrm{VQE}}\rangle=\prod_{k}\exp(-\hat{K}^{(k)})\exp(\hat{T}^{(k)})\exp(+\hat{K}^{(k)})|\Phi_{I}\rangle,$
(17)

where $\hat{K}^{(k)}$ and $\hat{T}^{(k)}$ are one- and two-body operators, and
$k$ is a parameter that controls the depth of the circuit and as a result its
variational freedom. The key difference in our $k$-uCJ ansatz from Ref. 61 is
that the two-body operator and the restriction to real anti-symmetric
matrices.

The one-body orbital transformations (spin restricted) are definad as

$\hat{K}^{(k)}\equiv\sum_{pp^{\prime}}\kappa_{pp^{\prime}}^{(k)}\left[(p^{\dagger}p^{\prime}-p^{\prime\dagger}p)+(\bar{p}^{\dagger}\bar{p}^{\prime}-\bar{p}^{\prime\dagger}\bar{p})\right]$
(18)

where $\kappa_{pp^{\prime}}^{(k)}=-\kappa_{p^{\prime}p}^{(k)}$ is a real,
antisymmetric $N_{a}\times N_{a}$ matrix of orbital rotation generators, which
equivalent to the one-particle spin-restricted orbital transformation:

$U_{pp^{\prime}}^{(k)}\equiv\left[\exp(\kappa^{(k)})\right]_{pp^{\prime}}$
(19)

This spin-restricted orbital rotation is expressed as quantum circuits using a
fabric of Givens rotations Kivlichan  _et al._ (2018).

The two-particle operator is defined as

$\begin{split}\hat{T}^{(k)}\equiv\sum_{p=0}^{M-1}\sum_{\begin{subarray}{c}p^{\prime}=p\
\mathrm{mod}\ 2\\\
p+=2\end{subarray}}^{M-2}\tau_{pp^{\prime}}^{(k)}[&(p^{\prime}+1)^{\dagger}\overline{(p^{\prime}+1)}^{\dagger}p^{\prime}\bar{p}^{\prime}\\\
&-p^{\prime\dagger}\bar{p}^{\prime\dagger}{(p^{\prime}+1)}\overline{(p^{\prime}+1)}]\end{split}$
(20)

The uCJ implementation is similar but not exact to Refs. Anselmetti  _et al._
(2021) and O’Gorman  _et al._ (2019). This it is denoted as $k$-muCJ for
clarity, with the ‘m’ standing for modified. It is important to point out that
the choice of VQE ansatz is largely irrelevant from a SAPT perspective and is
not a major point in this paper. An example of one layer of the muCJ circuit
ansatz is given in Fig. S1.

@C=1em @R=.7em & G [2] G [2] 3P_X G [2] G [2]

G [2] G [2] P_X G [2] G [2]

G G P_X G G

G G P_X G G

Figure S1: Quantum circuit of a single layer ($k=1$) $k$-muCJ VQE entangler
circuit for $M=2$ spatial orbitals or $N=4$ qubits. Even (odd) qubits
represent $\alpha$ ($\beta$) spin-orbitals. The quantum circuit starts with
two-quibit Givens rotation among $\alpha$ and $\beta$ orbitals. The next steps
are a double substitution operator (four qubit exchange gate) and another
layer of Givens rotations.

Using the ansatz discussed above, the VQE objective function is defined as

$\displaystyle E_{\mathrm{VQE}}(\kappa_{pq}^{k},\tau_{pq}^{k})$
$\displaystyle\equiv\langle\Psi_{\mathrm{VQE}}(\kappa_{pq}^{k},\tau_{pq}^{k})|\hat{H}|\Psi_{\mathrm{VQE}}(\kappa_{pq}^{k},\tau_{pq}^{k})\rangle$
(21)
$\displaystyle=\langle\Phi_{\mathrm{I}}|\hat{U}^{\dagger}(\kappa_{pq}^{k},\tau_{pq}^{k})\hat{H}\hat{U}(\kappa_{pq}^{k},\tau_{pq}^{k})|\Phi_{\mathrm{I}}\rangle.$
(22)

The number of measurements scales $\mathcal{O}(N_{a}^{4})$; however, we point
out that efficient estimation of the density matrices are currently
intensively investigated by several reaseach groups, e.g. Ref. 96 or Ref. 97.

##  Appendix B Derivation of the Second Order SAPT Terms

It is helpful to first recall standard Rayleigh-Schrödinger perturbation
theory for the intermolecular interaction energy. To begin we can write we
write the Hamiltonian of the combined system as

$\hat{H}=\hat{H}_{A}+\hat{H}_{B}+\hat{V},$ (23)

where we assume $\hat{H}_{X}\lvert\Psi_{X}\rangle=E_{X}\lvert\Psi_{X}\rangle$,
where $\lvert\Psi_{X}\rangle$ is the ground state wavefunction of monomer $X$
and $\hat{V}$ contains only the Coulombic interactions between monomer $A$ and
$B$:

$\hat{V}=\sum_{i=1}^{N_{A}}\sum_{j=1}^{N_{B}}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}-\sum_{i=1}^{N_{A}}\sum_{\beta}\frac{Z_{\beta}}{|\mathbf{R}_{\beta}-\mathbf{r}_{i}|}-\sum_{j}^{N_{B}}\sum_{\alpha}\frac{Z_{\alpha}}{|\mathbf{R}_{\alpha}-\mathbf{r}_{j}|}+\sum_{\alpha}\sum_{\beta}\frac{Z_{\alpha}Z_{\beta}}{|\mathbf{R}_{\alpha}-\mathbf{R}_{\beta}|},$
(24)

where it is understood that $i$ and $j$ are distinct indices and the sum over
$\alpha/\beta$ runs over the number of atoms in monomer $A/B$ respectively.
Under the assumption that the potential $\hat{V}$ is a small perturbation
(valid for weak intermolecular interactions) we can in principle now use
perturbation theory. Following Ref. 42 one typically starts SAPT by first
constructing a perturbation series for the so-called polarization energy

$|\Psi_{\mathrm{pol}}^{(n)}\rangle=-\hat{R}_{0}\hat{V}|\Psi^{(n)}_{\mathrm{pol}}\rangle+\sum_{k}^{n-1}E_{\mathrm{pol}}^{(k)}\hat{R}_{0}|\Psi_{\mathrm{pol}}^{(n-k)}\rangle,$
(25)

where

$\hat{R}_{0}=(\hat{H}_{0}-E_{0}+\hat{P}_{0})^{-1}\hat{Q}_{0},$ (26)

is the resolvent operator, $\hat{H}_{0}=\hat{H}_{A}+\hat{H}_{B}$,
$\hat{P}=|\Psi_{0}\rangle\langle\Psi_{0}|$, $\hat{Q}=1-\hat{P}$ and
$|\Psi_{\mathrm{pol}}^{0}\rangle=|\Psi_{0}\rangle=|\Psi_{A}^{0}\Psi_{B}^{0}\rangle$.
Here we have that $\hat{H}=\hat{H}_{A}+\hat{H}_{B}+\hat{V}$ and we assume
$\hat{H}_{0}|\Psi_{0}\rangle=(E^{0}_{A}+E^{0}_{B})|\Psi_{0}\rangle$. The first
order polarization energy (usually called the electrostatic energy
$E_{\mathrm{elst}}$) is then given by

$\displaystyle E_{\mathrm{pol}}^{(1)}$
$\displaystyle=\langle\Psi_{0}|\hat{V}|\Psi_{0}\rangle$ (27)
$\displaystyle=\bar{\gamma}_{pp^{\prime}}v_{pq}^{p^{\prime}q^{\prime}}\bar{\gamma}_{qq^{\prime}},$
(28)

where we used the definition
$E_{\mathrm{pol}}^{(n)}=\langle\Psi_{0}|\hat{V}|\Psi_{\mathrm{pol}}^{(n)}\rangle$

For the second order polarisation energy it is conventional in SAPT to split
into the induction and dispersion energy, i.e.,
$E_{\mathrm{pol}}^{(2)}=E_{\mathrm{ind}}+E_{\mathrm{disp}}$. The distinction
arises by considering that the induction term accounts for contributions from
terms where one monomer is in an excited state and the dispersion term
considers contributions from terms where both monomers are in excited states.
For the induction energy we have
$|\Psi_{A}^{(\mathrm{ind})}\rangle\equiv\hat{R}_{0}^{A}\hat{\Omega}^{B}|\Psi_{A}^{0}\rangle$,
where $\hat{R}_{0}^{A}=(\hat{H}_{A}-E_{0}^{A}+\hat{P}_{A})^{-1}\hat{Q}_{A}$
and

$\hat{\Omega}^{\mathrm{B}}=\sum_{i}^{N_{A}}\left(\hat{V}_{B}(\mathbf{r}_{i})+\int
d\mathbf{r}_{j}\frac{\hat{\rho}_{B}(\mathbf{r}_{j})}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}\right)$
(29)

is the effective electrostatic field of monomer $B$. Using the sum over states
formula for the resolvent we have then that

$\displaystyle|\Psi_{A}^{(\mathrm{ind})}\rangle$
$\displaystyle\equiv\hat{R}_{0}^{A}\hat{\Omega}^{B}|\Psi_{A}^{0}\rangle$ (30)
$\displaystyle=-\sum_{\mu}\frac{|\Psi^{\mu}_{A}\rangle\langle\Psi^{\mu}_{A}|\hat{\Omega}^{B}|\Psi_{A}^{0}\rangle}{E_{A}^{\mu}-E_{A}^{0}},$
(31)

from which it follows

$E_{\mathrm{ind}}(A\leftarrow
B)=-\sum_{\mu}\frac{|\langle\Psi^{\mu}_{A}|\hat{\Omega}^{B}|\Psi_{A}^{0}\rangle|^{2}}{E_{A}^{\mu}-E_{A}^{0}},$
(32)

with a similar expression for $E_{\mathrm{ind}}(B\leftarrow A)$. By
introducing the (spin-summed) transition one-particle reduced density matrix

$\bar{\gamma}_{pp^{\prime}}^{\mu}=\langle\Psi^{0}|p^{\dagger}p^{\prime}|\Psi^{\mu}\rangle+\langle\Psi^{0}|\bar{p}^{\dagger}\bar{p}^{\prime}|\Psi^{\mu}\rangle$
(33)

we can write Eq. 32 as

$E_{\mathrm{ind}}(A\leftarrow
B)=-\sum_{\mu}\frac{(\bar{\gamma}_{pp^{\prime}}^{\mu}\Omega_{pp^{\prime}})^{2}}{E_{A}^{\mu}-E_{A}^{0}}.$
(34)

Similar expressions exist for monomer $B$.

For the dispersion contribution we first define

$|\Psi^{\mathrm{disp}}_{AB}\rangle=\frac{|\Psi_{A}^{\mu}\Psi_{B}^{\nu}\rangle\langle\Psi_{A}^{\mu}\Psi_{B}^{\nu}|\hat{V}_{ee}|\Psi_{A}^{0}\Psi_{B}^{0}\rangle}{E_{A}^{\mu}+E_{B}^{\nu}},$
(35)

so that

$\displaystyle E_{\mathrm{disp}}^{(2)}$
$\displaystyle=\sum_{\mu\nu}\frac{|\langle\Psi_{A}^{0}\Psi^{0}_{B}|\hat{V}_{ee}|\Psi_{A}^{\mu}\Psi_{B}^{\nu}\rangle|^{2}}{E_{A}^{\mu}+E_{B}^{\nu}}$
(36)
$\displaystyle=\sum_{\mu\nu}\frac{\left(\bar{\gamma}^{\mu}_{pp^{\prime}}v_{pq}^{p^{\prime}q^{\prime}}\bar{\gamma}^{\nu}_{qq^{\prime}}\right)^{2}}{E_{A}^{\mu}+E_{B}^{\nu}}$
(37)

Equations 28, 34 and 37 provide all polarization energy contributions up to
second order in the intermolecular interaction energy. Unfortunately, without
further modification one finds the RS perturbation theory does not converge
for many-electron systems(Jeziorski  _et al._ , 1994). Symmetry adapted
perturbation theory (SAPT) attempts to fix this issue by explicitly accounting
for fermionic anti-symmetry in the wavefunction when electrons undergo
exchange processes between monomers.

The SAPT expression for the intermolecular interaction energy is given
by(Jeziorski  _et al._ , 1994)

$E_{\mathrm{SAPT}}^{(n)}=\frac{\langle\Psi_{0}|\hat{V}|\mathcal{A}\Psi_{\mathrm{pol}}^{(n-1)}\rangle-\sum_{k=1}^{n-1}E_{\mathrm{SAPT}}^{(n)}\langle\Psi_{0}|\mathcal{A}\Psi_{\mathrm{pol}}^{(n-k)}\rangle}{\langle\Psi_{0}|\mathcal{A}\Psi_{0}\rangle},$
(38)

where $|\Psi_{\mathrm{pol}}^{(n)}\rangle$ are the $n$th-order polarization
wavefunctions given previously and $\mathcal{A}$ is the antisymmetrizer
operator. With this definition the SAPT interaction takes the following form

$E_{\mathrm{SAPT}}^{(n)}=E_{\mathrm{pol}}^{(n)}+E_{\mathrm{exch}}^{(n)}.$ (39)

If one neglects all electron exchange processes other than those that exchange
a single electron between monomer $A$ and $B$ one arrives at the so-called
$S^{2}$ approximation to the exchange energies (through second order)(Korona,
2008c, 2009a; Hapka  _et al._ , 2019, 2021b):

$\displaystyle E_{\mathrm{exch}}^{(1)}$
$\displaystyle=\langle\Psi_{0}|(\hat{V}-\bar{V})(\hat{P}-\bar{P})|\Psi_{0}\rangle$
(40) $\displaystyle E_{\mathrm{exch-ind}}^{(2)}(B\rightarrow A)$
$\displaystyle=\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})(\hat{P}-\bar{P})|\Psi_{A}^{\mathrm{ind}}\Psi_{B}^{0}\rangle,$
(41) $\displaystyle E_{\mathrm{exch-ind}}^{(2)}$
$\displaystyle=\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})(\hat{P}-\bar{P})|\Psi^{\mathrm{disp}}_{AB}\rangle,$
(42)

where $|\Psi^{\mathrm{ind}}\rangle$ and $|\Psi^{\mathrm{disp}}_{AB}\rangle$
have been defined previously and $\hat{P}$ is following electron exchange
operator

$\hat{P}=\sum_{i}^{N_{A}}\sum_{j}^{N_{B}}\hat{P}_{ij}$ (43)

where $\hat{P}_{ij}$ exchanges an electron from monomer $A$ to monomer $B$.
Note a similar expression exists for $E_{\mathrm{exch-ind}}(A\rightarrow B)$.

It is helpful to note that Eqs. 40, 41 and 42 all share a similar structure.
Let us first example Eq. 40:

$E_{\mathrm{exch}}=\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})\hat{P}|\Psi^{0}_{A}\Psi_{B}^{0}\rangle$
(44)

which can be written in the density matrix formalism of SAPT as

$E_{\mathrm{exch}}^{(1)}(S^{2})=\int\gamma_{\mathrm{int}}(\mathbf{x}_{i},\mathbf{x}_{j})\left(\tilde{v}(\mathbf{r}_{i},\mathbf{r}_{j})-\frac{E_{\mathrm{pol}}^{(1)}}{N_{A}N_{B}}\right)d\mathbf{x}_{i}d\mathbf{x}_{j},$
(45)

where $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ denote both spin and spatial
coordinates of electrons in monomer $A$ and $B$ respectively and the
‘interaction’ density matrix is defined as

$\begin{split}&\gamma_{\mathrm{int}}(\mathbf{x}_{i},\mathbf{x}_{j})=\\\
&-\gamma_{A}(\mathbf{x}_{i},\mathbf{x}_{j})\gamma_{B}(\mathbf{x}_{j},\mathbf{x}_{i})\\\
&-\int\gamma_{A}(\mathbf{x}_{i},\mathbf{x}_{j^{\prime}})\Gamma_{B}(\mathbf{x}_{j},\mathbf{x}_{j^{\prime}},\mathbf{x}_{j},\mathbf{x}_{i})d\mathbf{x}_{j^{\prime}}\\\
&-\int\Gamma_{A}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}},\mathbf{x}_{i},\mathbf{x}_{j})\gamma_{B}(\mathbf{x}_{j},\mathbf{x}_{i^{\prime}})d\mathbf{x}_{i^{\prime}}\\\
&-\int\int\Gamma_{A}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}},\mathbf{x}_{i},\mathbf{x}_{j^{\prime}})\Gamma_{B}(\mathbf{x}_{j},\mathbf{x}_{j^{\prime}},\mathbf{x}_{j},\mathbf{x}_{i^{\prime}})d\mathbf{x}_{i^{\prime}}d\mathbf{x}_{j^{\prime}}.\end{split}$
(46)

By inspection Eq. 41 is _identical_ to Eq. 44 up to the replacement of
$|\Psi_{A}^{0}\rangle$ by the excited state wavefunction
$|\Psi_{A}^{\mu}\rangle$. Thus we have

$\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})\hat{P}|\Psi_{\mu}^{A}\Psi_{B}^{0}\rangle=\int\gamma^{\mu}_{\mathrm{int}}(\mathbf{x}_{i},\mathbf{x}_{j})\left(\tilde{v}(\mathbf{r}_{i},\mathbf{r}_{j})-\frac{E_{\mathrm{pol}}^{(1)}}{N_{A}N_{B}}\right)d\mathbf{x}_{i}d\mathbf{x}_{j},$
(47)

where

$\begin{split}&\gamma^{\mu}_{\mathrm{int}}(\mathbf{x}_{i},\mathbf{x}_{j})=\\\
&-\gamma^{\mu}_{A}(\mathbf{x}_{i},\mathbf{x}_{j})\gamma_{B}(\mathbf{x}_{j},\mathbf{x}_{i})\\\
&-\int\gamma^{\mu}_{A}(\mathbf{x}_{i},\mathbf{x}_{j^{\prime}})\Gamma_{B}(\mathbf{x}_{j},\mathbf{x}_{j^{\prime}},\mathbf{x}_{j},\mathbf{x}_{i})d\mathbf{x}_{j^{\prime}}\\\
&-\int\Gamma^{\mu}_{A}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}},\mathbf{x}_{i},\mathbf{x}_{j})\gamma_{B}(\mathbf{x}_{j},\mathbf{x}_{i^{\prime}})d\mathbf{x}_{i^{\prime}}\\\
&-\int\int\Gamma^{\mu}_{A}(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}},\mathbf{x}_{i},\mathbf{x}_{j^{\prime}})\Gamma_{B}(\mathbf{x}_{j},\mathbf{x}_{j^{\prime}},\mathbf{x}_{j},\mathbf{x}_{i^{\prime}})d\mathbf{x}_{i^{\prime}}d\mathbf{x}_{j^{\prime}}.\end{split}$
(48)

where $\gamma^{\mu}$ and $\Gamma^{\mu}$ are transition one- and two-particle
density matrices. Expanding the transition density matrix in orbital space and
performing spin summations we find

$\begin{split}\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})\hat{P}|\Psi^{\mu}_{A}\Psi_{B}^{0}\rangle&=-\frac{1}{2}(\bar{\gamma}^{\mu}_{pp^{\prime}}\bar{\gamma}_{qq^{\prime}}\tilde{v}_{pq}^{q^{\prime}p^{\prime}}\\\
&+\bar{\gamma}^{\mu}_{pp^{\prime}}\bar{\Gamma}_{q^{\prime\prime}q^{\prime\prime\prime}}^{qq^{\prime}}S_{p^{\prime}q^{\prime}}\tilde{v}_{pq}^{q^{\prime\prime\prime}q^{\prime\prime}}\\\
&+\bar{\gamma}_{qq^{\prime}}[\bar{\Gamma}^{\mu}]_{p^{\prime\prime}p^{\prime\prime\prime}}^{pp^{\prime}}S_{q^{\prime}p^{\prime}}\tilde{v}_{pq}^{p^{\prime\prime}p^{\prime\prime\prime}}\\\
&+[\bar{\Gamma}^{\mu}]_{p^{\prime\prime}p^{\prime\prime\prime}}^{pp^{\prime}}\bar{\Gamma}_{q^{\prime\prime}q^{\prime\prime\prime}}^{qq^{\prime}}S_{p^{\prime}q^{\prime\prime\prime}}S_{q^{\prime}p^{\prime\prime\prime}}\tilde{v}_{pq}^{p^{\prime\prime}q^{\prime\prime}}\\\
&-E_{\mathrm{pol}}^{(1)}S_{pq^{\prime}}\bar{\gamma}_{qq^{\prime}}S_{qp^{\prime}}\bar{\gamma}^{\mu}_{pp^{\prime}}).\end{split}$
(49)

It is helpful to define the exchange functional

$\begin{split}F^{\mu\nu}\equiv
F[\bar{\gamma}^{\mu},\bar{\Gamma}^{\mu},\bar{\gamma}^{\nu},\bar{\Gamma}^{\nu}]&=-\frac{1}{2}([\bar{\gamma}^{\mu}]_{pp^{\prime}}[\bar{\gamma}^{\nu}]_{qq^{\prime}}\tilde{v}_{pq}^{q^{\prime}p^{\prime}}\\\
&+[\bar{\gamma}^{\mu}]_{pp^{\prime}}[\bar{\Gamma}^{\nu}]_{q^{\prime\prime}q^{\prime\prime\prime}}^{qq^{\prime}}S_{p^{\prime}q^{\prime}}\tilde{v}_{pq}^{q^{\prime\prime\prime}q^{\prime\prime}}\\\
&+[\bar{\gamma}^{\nu}]_{qq^{\prime}}[\bar{\Gamma}^{\mu}]_{p^{\prime\prime}p^{\prime\prime\prime}}^{pp^{\prime}}S_{q^{\prime}p^{\prime}}\tilde{v}_{pq}^{p^{\prime\prime}p^{\prime\prime\prime}}\\\
&+[\bar{\Gamma}^{\mu}]_{p^{\prime\prime}p^{\prime\prime\prime}}^{pp^{\prime}}[\bar{\Gamma}^{\nu}]_{q^{\prime\prime}q^{\prime\prime\prime}}^{qq^{\prime}}S_{p^{\prime}q^{\prime\prime\prime}}S_{q^{\prime}p^{\prime\prime\prime}}\tilde{v}_{pq}^{p^{\prime\prime}q^{\prime\prime}}\\\
&-E_{\mathrm{pol}}^{(1)}S_{pq^{\prime}}[\bar{\gamma}^{\nu}]_{qq^{\prime}}S_{qp^{\prime}}[\bar{\gamma}^{\mu}]^{pp^{\prime}}),\end{split}$
(50)

with the understanding that $\gamma^{0}=\gamma$ etc., so that for the
exchange-dispersion contribution we have

$\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})\hat{P}|\Psi_{A}^{\mu}\Psi_{B}^{\nu}\rangle=F^{\mu\nu}=F[\bar{\gamma}_{A}^{\mu},\bar{\Gamma}_{A}^{\mu},\bar{\gamma}_{B}^{\nu},\bar{\Gamma}^{\nu}_{B}],$
(51)

which follows by inspection.

With these notational conventions we can evaluate exchange-induction term we
have (in the $S^{2}$ approximation)

$E_{\mathrm{exch-ind}}(B\rightarrow
A)=\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})(\hat{P}-\bar{P})|\Psi_{A}^{(\mathrm{ind})}\Psi_{B}^{0}\rangle.$
(52)

Inserting Eq. 31 into Eq. 52 and breaking the expectation value into two terms
we find (dropping $B\rightarrow A$ notation)

$E_{\mathrm{exch-
ind}}=-\sum_{\mu}\frac{\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})(\hat{P}-\bar{P})|\Psi^{\mu}_{A}\Psi_{B}^{0}\rangle\langle\Psi^{\mu}_{A}\Psi_{B}^{0}|\hat{\Omega}^{B}|\Psi_{A}^{0}\Psi_{B}^{0}\rangle}{E_{A}^{\mu}-E_{A}^{0}}.$
(53)

Let us first look at the term in the numerator,
$\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})(\hat{P}-\bar{P})|\Psi^{\mu}_{A}\Psi_{B}^{0}\rangle$
which we will break down into two pieces:

$\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})\hat{P}|\Psi^{\mu}_{A}\Psi_{B}^{0}\rangle$
(54)

and

$\bar{P}\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})|\Psi^{\mu}_{A}\Psi_{B}^{0}\rangle.$
(55)

Eq. 54 is given by $F^{\mu}$ as discussed previously while Eq. 55 is simply
(note $\langle\Psi_{A}^{0}|\Psi_{A}^{\mu}\rangle=0$ so the $\bar{P}\bar{V}$
contribution is zero)

$\displaystyle\bar{P}\langle\Psi_{A}^{0}\Psi_{B}^{0}|\hat{V}|\Psi^{\mu}_{A}\Psi_{B}^{0}\rangle$
$\displaystyle=\bar{P}\langle\Psi_{A}^{0}\Psi_{B}^{0}|\hat{\Omega}|\Psi^{\mu}_{A}\Psi_{B}^{0}\rangle$
(56)
$\displaystyle=\bar{P}\bar{\gamma}_{pp^{\prime}}^{\mu}\Omega_{pp^{\prime}}^{B}$
(57)

and by inspection with Eq. 45

$\displaystyle\bar{P}\equiv\langle\Psi_{A}^{0}\Psi_{B}^{0}|\hat{P}|\Psi_{A}^{0}\Psi_{B}^{0}\rangle$
$\displaystyle=\frac{1}{N_{A}N_{B}}\int d\mathrm{x}\ d\mathrm{x}^{\prime}\
\gamma_{\mathrm{int}}(\mathbf{x},\mathbf{x}^{\prime})$ (58)
$\displaystyle=-\frac{1}{2}S_{pq^{\prime}}\bar{\gamma}_{qq^{\prime}}S_{qp^{\prime}}\bar{\gamma}_{pp^{\prime}}.$
(59)

So, putting all of this together we have

$\displaystyle E_{\mathrm{exch-ind}}(B\rightarrow
A)=-\left(\frac{F^{\mu}t^{\mu}}{\omega_{A}^{\mu}}+E_{\mathrm{ind}}(B\rightarrow
A)\bar{P}\right),$ (60)

where

$t^{\mu}=\gamma_{pp^{\prime}}^{\mu}\Omega^{B}_{pp^{\prime}},$ (61)

and we have used the fact that

$\frac{(\bar{\gamma}^{\mu}_{pp^{\prime}}\Omega^{B}_{pp^{\prime}})^{2}}{\omega^{\mu}_{A}}\equiv-
E_{\mathrm{ind}}(B\rightarrow A).$ (62)

The full expression is then

$E_{\mathrm{exch-ind}}=E_{\mathrm{exch-ind}}(B\rightarrow A)+E_{\mathrm{exch-
ind}}(A\rightarrow
B)=-\left(\frac{F^{\mu}t^{\mu}}{\omega_{A}^{\mu}}+\frac{F^{\nu}t^{\nu}}{\omega_{B}^{\nu}}+E_{\mathrm{ind}}\bar{P}\right).$
(63)

The second-order exchange-dispersion energy can be found in a similar fashion.
We can start from the definition in the $S^{2}$ approximation

$E_{\mathrm{exch-
disp}}(S^{2})=\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})(\hat{P}-\bar{P})|\Psi^{\mathrm{disp}}_{AB}\rangle,$
(64)

In contrast to the induction term, the dispersion interaction includes terms
where both monomers are in an excited state. Much like the exchange-induction
interaction we can split Eq. 64 into two terms. Inserting Eq. Eq. 35 into Eq.
64 we first have to evaluate

$\langle\Psi_{A}^{0}\Psi_{B}^{0}|(\hat{V}-\bar{V})\hat{P}|\Psi_{A}^{\mu}\Psi_{B}^{\nu}\rangle=F^{\mu\nu}=F[\bar{\gamma}_{A}^{\mu},\bar{\Gamma}_{A}^{\mu},\bar{\gamma}_{B}^{\nu},\bar{\Gamma}^{\nu}_{B}],$
(65)

which follows by inspection. We also identify the usual induction numerator
term

$\langle\Psi_{A}^{0}\Psi_{B}^{0}|V_{ee}|\Psi_{A}^{\mu}\Psi_{B}^{\nu}\rangle=v_{pq}^{p^{\prime}q^{\prime}}\gamma_{pp^{\prime}}^{\mu}\gamma_{qq^{\prime}}^{\nu}=t^{\mu\nu}.$
(66)

Noting that

$E^{(2)}_{\mathrm{disp}}=-\sum_{\mu\nu}\frac{(t^{\mu\nu})^{2}}{\omega_{A}^{\mu}+\omega_{B}^{\nu}}$
(67)

and combining all of these with Eq. 59 we find

$E_{\mathrm{exch-
disp}}=-\left(\frac{F^{\mu\nu}t^{\mu\nu}}{\omega_{A}^{\mu}+\omega_{B}^{\nu}}+E_{\mathrm{disp}}\bar{P}\right).$
(68)

##  Appendix C Extended Random Phase Approximation

Evaluating the SAPT expressions given in Appendix B requires a knowledge of
the one- and two-particle transition density matrices at a given level of
theory. There are two challenges to achieving this given current NISQ
hardware. First is that we can only simulate a small fragment of the problem
in a quantum computer, which is achieved here using an active space approach.
Second is that in principle we would need to compute all eigenvalues and
excited states to compute the transition density matrices which is
significantly more challenging than just computing the ground state energy.

To overcome these issues, we instead include only the subset of particle-hole
excitations. In particular, we follow Hapka et al. (Hapka  _et al._ , 2018,
2019, 2021b) and use the extended random phase approximation (Chatterjee and
Pernal, 2012; van Aggelen  _et al._ , 2013; Pernal  _et al._ , 2014) (ERPA) to
approximately determined excited state properties from CASSCF quality
wavefunctions. For notational clarity in this section we will not distinguish
between monomers, and orbitals are general spatial orbitals.

Following Pernal (Chatterjee and Pernal, 2012; Pernal  _et al._ , 2014), in
the ERPA solves the following eigenvalue problem

$\begin{bmatrix}A&B\\\ B&A\end{bmatrix}\begin{bmatrix}X_{\nu}\\\
Y_{\nu}\end{bmatrix}=\omega_{\nu}\begin{bmatrix}-N&0\\\
0&N\end{bmatrix}\begin{bmatrix}X_{\nu}\\\ Y_{\nu}\end{bmatrix}.$ (69)

where

$\displaystyle A_{pq,rs}$
$\displaystyle=\langle\Psi_{0}|[p^{\dagger}q,[\hat{H},s^{\dagger}r]|\Psi_{0}\rangle$
(70) $\displaystyle B_{pq,rs}$
$\displaystyle=\langle\Psi_{0}|[p^{\dagger}q,[\hat{H},r^{\dagger}s]|\Psi_{0}\rangle$
(71) $\displaystyle N_{pq,rs}$
$\displaystyle=\langle\Psi_{0}|[p^{\dagger}q,s^{\dagger}r]|\Psi_{0}\rangle=\delta_{ps}\delta_{qr}(n_{q}-n_{p}),$
(72)

where the above equations are valid in the natural spin-orbital basis
corresponding to $|\Psi_{0}\rangle$ and we take $p>q$ and $r>s$. Computing the
Hessian matrices in Eq. 69 is somewhat tedious but expressions are available
in the literature (Chatterjee and Pernal, 2012; Pernal  _et al._ , 2014).

Note that using excited states from the ERPA introduces two approximations in
the quality of the resulting SAPT interaction energies. The first, is that
even with an exact FCI wavefunction the ERPA will not yield all the excited
states of the Hamiltonian as only particle-hole-like excitations are included.
The second issue is that the CASSCF wavefunction is not an eigenstate of the
zeroth order Hamiltonian in the SAPT perturbation theory. In principle, this
issue can be addressed within the context of the adiabatic-connection
formalism based upon the ERPA (Pernal, 2018; Pastorczak and Pernal, 2018;
Hapka  _et al._ , 2018, 2019, 2021b) which yields the so-called coupled
approximation used in this work. This amounts to using the full Hamiltonian in
Eq. 69 with the CASSCF (or VQE) wavefunction as $|\Psi_{0}\rangle$. The
adiabatic-connection approach assumes the one- and two-particle reduced
density matrices remain constant across the adiabatic-connection pathway. This
approximation appears to perform well in the sense that SAPT(CASSCF) compares
favourably to FCI interaction energies (Hapka  _et al._ , 2018, 2019, 2021b).

###  C.1 Spin-Summation

In this work we only target singlet excited states and thus
$X_{\beta\alpha}=X_{\alpha\beta}=0$, thus we explicitly form the spin-
restricted ERPA equations. In particular we form

$\displaystyle A_{pq,rs}$
$\displaystyle=\langle\Psi_{0}|\hat{E}_{pq},[\hat{H},\hat{E}_{sr}]|\Psi_{0}\rangle$
(73) $\displaystyle B_{pq,rs}$
$\displaystyle=\langle\Psi_{0}|[\hat{E}_{pq},[\hat{H},\hat{E}_{rs}]|\Psi_{0}\rangle$
(74) $\displaystyle N_{pq,rs}$
$\displaystyle=\langle\Psi_{0}|[\hat{E}_{pq},\hat{E}_{sr}]|\Psi_{0}\rangle=\delta_{ps}\delta_{qr}(n_{q}-n_{p}),$
(75)

where $\hat{E}_{rs}=(r^{\dagger}s+\bar{r}^{\dagger}\bar{s})$. It is helpful to
note the following identities:

$\displaystyle[A,BC]$ $\displaystyle=A\\{B,C\\}-\\{A,C\\}B$ (76)
$\displaystyle[AB,CD]$ $\displaystyle=[AB,C]D+C[AB,D]$ (77)

So, for example,

$\begin{split}[\hat{E}_{rs},\hat{E}_{sr}]&=\sum_{pq^{\prime}}[r^{\dagger}s+\bar{r}^{\dagger}\bar{s},q^{\dagger}p+\bar{q}^{\dagger}\bar{p}]\\\
&=\sum_{pq^{\prime}}[r^{\dagger}s,q^{\dagger}p]+[r^{\dagger}s,\bar{q}^{\dagger}\bar{p}]+[\bar{r}^{\dagger}\bar{s},q^{\dagger}p]+[\bar{r}^{\dagger}\bar{s},\bar{q}^{\dagger}\bar{p}],\end{split}$
(78)

and

$\displaystyle[r^{\dagger}s,q^{\dagger}p]$
$\displaystyle=[r^{\dagger}s,q^{\dagger}]p+q^{\dagger}[r^{\dagger}s,p]$ (79)
$\displaystyle=(r^{\dagger}\\{s,q^{\dagger}\\}-\\{r^{\dagger},q^{\dagger}\\}s)p+q^{\dagger}(r^{\dagger}\\{s,p\\}-\\{r^{\dagger},p\\}s)$
(80) $\displaystyle=r^{\dagger}p\delta_{sq}-q^{\dagger}s\delta_{rp}$ (81)
$\displaystyle[\bar{r}^{\dagger}\bar{s},\bar{q}^{\dagger}\bar{p}]$
$\displaystyle=\bar{r}^{\dagger}\bar{p}\delta_{\bar{s}\bar{q}}-\bar{q}^{\dagger}\bar{s}\delta_{\bar{r}\bar{p}}$
(82) $\displaystyle[\bar{r}^{\dagger}\bar{s},q^{\dagger}p]$
$\displaystyle=[r^{\dagger}s,\bar{q}^{\dagger}\bar{p}]=0$ (83)

Similarly(Helgaker  _et al._ , 2014),

$\displaystyle[\hat{E}_{mn},\hat{E}_{pq}]$
$\displaystyle=\hat{E}_{mq}\delta_{pn}-\hat{E}_{pn}\delta_{mq}$ (84)
$\displaystyle[\hat{E}_{rs},[\hat{E}_{ij},\hat{E}_{qp}]]$
$\displaystyle=[\hat{E}_{rs},\hat{E}_{ip}]\delta_{qj}-[\hat{E}_{rs},\hat{E}_{qj}]\delta_{ip}$
(85)
$\displaystyle=\hat{E}_{rp}\delta_{is}\delta_{qj}-\hat{E}_{is}\delta_{rp}\delta_{qj}-\hat{E}_{rj}\delta_{qs}\delta_{ip}+\hat{E}_{qs}\delta_{rj}\delta_{ip}$
(86) $\displaystyle[\hat{E}_{mn},\hat{e}_{pqrs}]$
$\displaystyle=\delta_{pn}\hat{e}_{mqrs}-\delta_{mq}\hat{e}_{pnrs}+\delta_{rn}\hat{e}_{pqms}-\delta_{ms}\hat{e}_{pqrn},$
(87)

where
$\hat{e}_{pqrs}=p^{\dagger}r^{\dagger}sq+p^{\dagger}\bar{r}^{\dagger}\bar{s}q+\bar{p}^{\dagger}r^{\dagger}s\bar{q}+\bar{p}^{\dagger}\bar{r}^{\dagger}\bar{s}\bar{q}$.
We can now evaluate the commutator on the left hand side of Eq. 73 where we
take

$\hat{H}=\sum_{pq}h_{pq}(p^{\dagger}q+\bar{p}^{\dagger}\bar{q})+\frac{1}{2}\sum_{pqrs}(pq|rs)\hat{e}_{pqrs}.$
(88)

Let us begin with the one-body term (dropping hats on operators for brevity,
and here all orbital labels refer to general spatial orbitals, no distinction
is made between active, core or inactive)

$\displaystyle\sum_{ij}h_{ij}[E_{rs},[E_{ij},E_{qp}]]$
$\displaystyle=\sum_{ij}h_{ij}(\hat{E}_{rp}\delta_{is}\delta_{qj}-\hat{E}_{is}\delta_{rp}\delta_{qj}-\hat{E}_{rj}\delta_{qs}\delta_{ip}+\hat{E}_{qs}\delta_{rj}\delta_{ip})$
(89)
$\displaystyle=h_{sq}\hat{E}_{rp}-\sum_{i}h_{iq}\hat{E}_{is}\delta_{rp}-\sum_{j}h_{pj}\hat{E}_{rj}\delta_{qs}+h_{pr}\hat{E}_{qs}$
(90)

If we work in the natural orbital basis
$\langle\Psi_{0}|\hat{E}_{pq}|\Psi_{0}\rangle=\bar{\gamma}_{pq}=n_{p}\delta_{pq}$,
with $0\leq n_{p}\leq 2$111Note we assume $n_{p}<n_{q}$ for $p>q$ (i.e. the
basis is ordered in descending order of natural orbital occupancy) which is
the opposite ordering that is usually taken in the ERPA literature (Chatterjee
and Pernal, 2012; Pernal  _et al._ , 2014), then we have

$\displaystyle\langle\sum_{ij}h_{ij}[E_{rs},[E_{ij},E_{qp}]]\rangle$
$\displaystyle=h_{sq}\delta_{rp}(n_{r}-n_{s})+h_{pr}\delta_{qs}(n_{q}-n_{r}).$
(91)

For the two-body part we have to evaluate expressions like (see Eq. 74)

$\displaystyle\sum_{ijkl}(ij|kl)[E_{rs},[e_{ijkl},E_{pq}]]$ (92)

Let us look at

$\displaystyle[E_{rs},[e_{ijkl},E_{pq}]]$
$\displaystyle=-[E_{rs},[E_{pq},e_{ijkl}]]=-[E_{rs},\delta_{iq}\hat{e}_{pjkl}-\delta_{pj}\hat{e}_{iqkl}+\delta_{kq}\hat{e}_{ijpl}-\delta_{pl}\hat{e}_{ijkq}]$
(93)
$\displaystyle=-[E_{rs},\delta_{iq}\hat{e}_{pjkl}-\delta_{pj}\hat{e}_{iqkl}+\delta_{kq}\hat{e}_{ijpl}-\delta_{pl}\hat{e}_{ijkq}]$
(94) $\displaystyle[\hat{E}_{rs},\delta_{iq}\hat{e}_{pjkl}]$
$\displaystyle=\delta_{iq}(\delta_{ps}\hat{e}_{rjkl}-\delta_{rj}\hat{e}_{pskl}+\delta_{ks}\hat{e}_{pjrl}-\delta_{rl}\hat{e}_{pjks})$
(95) $\displaystyle[\hat{E}_{rs},\delta_{pj}\hat{e}_{iqkl}]$
$\displaystyle=\delta_{pj}(\delta_{is}\hat{e}_{rqkl}-\delta_{rq}\hat{e}_{iskl}+\delta_{ks}\hat{e}_{iqrl}-\delta_{rl}\hat{e}_{iqks})$
(96) $\displaystyle[\hat{E}_{rs},\delta_{kq}\hat{e}_{ijpl}]$
$\displaystyle=\delta_{kq}(\delta_{is}\hat{e}_{rjpl}-\delta_{rj}\hat{e}_{ispl}+\delta_{ps}\hat{e}_{ijrl}-\delta_{rl}\hat{e}_{ijps})$
(97) $\displaystyle[\hat{E}_{rs},\delta_{pl}\hat{e}_{ijkq}]$
$\displaystyle=\delta_{pl}(\delta_{is}\hat{e}_{rjkq}-\delta_{rj}\hat{e}_{iskq}+\delta_{ks}\hat{e}_{ijrq}-\delta_{rq}\hat{e}_{ijks})$
(98)

so,

$\begin{split}\langle\sum_{ijkl}(ij|kl)[E_{rs},[e_{ijkl},E_{pq}]]\rangle=-\sum_{ijkl}(ij|kl)&\bigg{[}\delta_{iq}(\delta_{ps}\bar{\Gamma}_{rjkl}-\delta_{rj}\bar{\Gamma}_{pskl}+\delta_{ks}\bar{\Gamma}_{pjrl}-\delta_{rl}\bar{\Gamma}_{pjks})\\\
&-\delta_{pj}(\delta_{is}\bar{\Gamma}_{rqkl}-\delta_{rq}\bar{\Gamma}_{iskl}+\delta_{ks}\bar{\Gamma}_{iqrl}-\delta_{rl}\bar{\Gamma}_{iqks})\\\
&+\delta_{kq}(\delta_{is}\bar{\Gamma}_{rjpl}-\delta_{rj}\bar{\Gamma}_{ispl}+\delta_{ps}\bar{\Gamma}_{ijrl}-\delta_{rl}\bar{\Gamma}_{ijps})\\\
&-\delta_{pl}(\delta_{is}\bar{\Gamma}_{rjkq}-\delta_{rj}\bar{\Gamma}_{iskq}+\delta_{ks}\bar{\Gamma}_{ijrq}-\delta_{rq}\bar{\Gamma}_{ijks})\bigg{]},\end{split}$
(99)

where the spin-summed two-particle reduced density matrix is

$\bar{\Gamma}_{pqrs}=\langle\Psi_{0}|\hat{e}_{pqrs}|\Psi_{0}\rangle.$ (100)

Removing the Kroenecker deltas we find

$\begin{split}\frac{1}{2}\langle\sum_{ijkl}(ij|kl)&[E_{rs},[e_{ijkl},E_{pq}]]\rangle=\\\
&-\frac{1}{2}\bigg{[}\sum_{jkl}(qj|kl)\delta_{ps}\bar{\Gamma}_{rjkl}-\sum_{kl}(qr|kl)\bar{\Gamma}_{pskl}+\sum_{jl}(qj|sl)\bar{\Gamma}_{pjrl}-\sum_{jk}(qj|kr)\bar{\Gamma}_{pjks}\\\
&-\sum_{kl}(sp|kl)\bar{\Gamma}_{rqkl}+\sum_{ikl}(ip|kl)\delta_{rq}\bar{\Gamma}_{iskl}-\sum_{il}(ip|sl)\bar{\Gamma}_{iqrl}+\sum_{ik}(ip|kr)\bar{\Gamma}_{iqks}\\\
&+\sum_{jl}(sj|ql)\bar{\Gamma}_{rjpl}-\sum_{il}(ir|ql)\bar{\Gamma}_{ispl}+\sum_{ijl}(ij|ql)\delta_{ps}\bar{\Gamma}_{ijrl}-\sum_{ij}(ij|qr)\bar{\Gamma}_{ijps}\\\
&-\sum_{jk}(sj|kp)\bar{\Gamma}_{rjkq}+\sum_{ik}(ir|kp)\bar{\Gamma}_{iskq}-\sum_{ij}(ij|sp)\bar{\Gamma}_{ijrq}+\sum_{ijk}(ij|kp)\delta_{rq}\bar{\Gamma}_{ijks}\bigg{]}\end{split}$
(101)

which can be simplified a using symmetries in the ERIs and TPDMs. Note Eq. 101
reduces to the expression given in Ref.54 if working in the natural spin
orbital basis.

###  C.2 Numerical Solution

With the expressions for the Hessians at hand we can now solve Eq. 69. First,
it is helpful to bring this expression into a symmetric eigenvalue problem of
half the dimension. This has two advantages, namely the computational cost
will be reduced by a factor of eight and we can use a symmetric eigenvalue
solver.

To proceed it is important to note some features Eq. 69. First we choose the
normalization:

$\displaystyle Y_{\nu}^{\dagger}NY_{\nu}-X_{\nu}^{\dagger}NX_{\nu}$
$\displaystyle=1,$ (102) $\displaystyle(Y-X)^{\dagger}N(Y+X)$
$\displaystyle=1.$ (103)

Next, we note that

$\displaystyle(A+B)(X+Y)$ $\displaystyle=\omega_{n}N(Y-X)$ (104)
$\displaystyle(A-B)(Y-X)$ $\displaystyle=\omega_{n}N(X+Y),$ (105)

which implies

$\displaystyle(A+B)(X+Y)=\omega_{n}^{2}N(A-B)^{-1}N(X+Y)$ (106)

or

$\displaystyle(A+B)N^{-1/2}N^{1/2}(X+Y)$
$\displaystyle=\omega_{n}^{2}N^{1/2}N^{1/2}(A-B)^{-1}N^{1/2}N^{1/2}(X+Y)$
(107) $\displaystyle A_{+}T_{n}$ $\displaystyle=\omega_{n}^{2}A_{-}^{-1}T_{n}$
(108) $\displaystyle A_{-}A_{+}T_{n}$ $\displaystyle=\omega_{n}^{2}T_{n}$
(109)

if $(A-B)$ is positive definite.

If $(A+B)$ is also positive definite we can write

$\displaystyle A_{+}^{1/2}A_{+}^{1/2}A_{-}A_{+}^{1/2}A_{+}^{1/2}T_{n}$
$\displaystyle=\omega_{n}^{2}A_{+}T_{n}$ (110) $\displaystyle
A_{+}^{1/2}A_{-}A_{+}^{1/2}A_{+}^{1/2}T_{n}$
$\displaystyle=\omega_{n}^{2}A_{+}^{1/2}T_{n}$ (111) $\displaystyle
A_{+}^{1/2}A_{-}A_{+}^{1/2}Q_{n}$ $\displaystyle=\omega_{n}^{2}Q_{n}$ (112)

where

$Q_{n}=A_{+}^{1/2}T_{n}=A_{+}^{1/2}N^{1/2}(X+Y),$ (114)

$Q^{\dagger}Q=1$, and

$X+Y=N^{-1/2}A_{+}^{-1/2}Q_{n}.$ (115)

This allows us to write combinations of the $X$ and $Y$ matrices necessary for
building transition density matrices in terms of the solutions of the
Hermitian eigenvalue problem. The factors of $N^{-1/2}$ can be determined
using canonical orthogonalization with a threshold to discard small
eigenvalues.

###  C.3 Normalization

A useful sanity check at this point is that we have

$\displaystyle(A+B)(X+Y)$ $\displaystyle=\omega_{n}N(Y-X).$ (116)

Multiplying on the right by $(X+Y)^{\dagger}$ we have and using Eq. 103

$\displaystyle(X+Y)_{n}(A+B)(X+Y)_{n}$ $\displaystyle=\omega_{n}(X+Y)N(Y-X)$
(117) $\displaystyle(X+Y)_{n}(A+B)(X+Y)_{n}$ $\displaystyle=\omega_{n}$ (118)

Inserting Eq. 115 in Eq. 118 we find

$\displaystyle(X+Y)_{n}(A+B)(X+Y)_{n}$
$\displaystyle=Q_{n}A_{+}^{-1/2}N^{-1/2}(A+B)N^{-1/2}A_{+}^{-1/2}Q_{n}=Q_{n}A_{+}^{-1/2}A_{+}A_{+}^{-1/2}Q_{n}=1\neq\omega!$
(119)

A more consistent normalization of Eq. 115 can be found by

$X+Y\equiv T_{+}=\sqrt{\omega_{n}}N^{-1/2}A_{+}^{-1/2}Q_{n}$ (120)

and by Eq. 103 we can deduce

$Y-X\equiv T_{-}=\frac{1}{\sqrt{\omega_{n}}}N^{-1/2}A_{+}^{1/2}Q_{n}$ (121)

from which it follows that

$\displaystyle X$ $\displaystyle=\frac{1}{2}(T_{+}-T_{-})$ (122)
$\displaystyle Y$ $\displaystyle=\frac{1}{2}(T_{+}+T_{-}).$ (123)

To compute the symmetrized (spin summed) transition density matrix we follow
Pernal (Chatterjee and Pernal, 2012; Pernal, 2014) and use the identity

$\displaystyle\bar{\gamma}^{\nu}_{pq}$
$\displaystyle=\langle\Psi_{0}|\hat{E}_{pq}|\Psi_{\nu}\rangle=\langle
0|[\hat{E}_{pq},\hat{O}^{\dagger}_{\nu}]|0\rangle,$ (124)

then

$\displaystyle\gamma_{pq}^{\nu}+\gamma_{qp}^{\nu}$
$\displaystyle=[N(Y-X)]_{pq}\qquad\forall_{p>q}.$ (125)

For the transition two-particle reduced density matrix we have

$\displaystyle\bar{\Gamma}^{\nu}_{pqrs}=\langle\Psi_{0}|[\hat{e}_{pqrs},\hat{O}_{\nu}^{\dagger}]|\Psi_{0}\rangle.$
(126)

Using Eq. 87 and looking at the $X$ part of $\hat{O}_{\nu}^{\dagger}$ first we
have

$\displaystyle\bar{\Gamma}^{\nu(X)}_{pqrs}$
$\displaystyle=-\sum_{m>n}X_{mn}^{\nu}\langle[\hat{E}_{mn},\hat{e}_{pqrs}]\rangle$
(127)
$\displaystyle=-\sum_{m>n}X_{mn}^{\nu}\left(\delta_{pn}\bar{\Gamma}_{mqrs}-\delta_{mq}\bar{\Gamma}_{pnrs}+\delta_{rn}\bar{\Gamma}_{pqms}-\delta_{ms}\bar{\Gamma}_{pqrn}\right)$
(128)

and similarly we have

$\displaystyle\bar{\Gamma}^{\nu(Y)}_{pqrs}$
$\displaystyle=-\sum_{m>n}Y_{mn}^{\nu}\langle[\hat{E}_{nm},\hat{e}_{pqrs}]\rangle$
(129)
$\displaystyle=-\sum_{m>n}Y_{mn}^{\nu}\left(\delta_{pm}\bar{\Gamma}_{nqrs}-\delta_{nq}\bar{\Gamma}_{pmrs}+\delta_{rm}\bar{\Gamma}_{pqns}-\delta_{ns}\bar{\Gamma}_{pqrm}\right)$
(130)

and naturally
$\bar{\Gamma}^{\nu}_{pqrs}=\bar{\Gamma}^{\nu(X)}_{pqrs}+\bar{\Gamma}^{\nu(Y)}_{pqrs}$.

It is helpful remove the Kroenecker deltas so that

$\begin{split}\bar{\Gamma}_{pqrs}^{\nu}=&-\left(\sum_{m=p+1}^{M}X^{\nu}_{mp}\bar{\Gamma}_{mqrs}+\sum_{m=0}^{p}Y^{\nu}_{pm}\bar{\Gamma}_{mqrs}\right)\\\
&+\left(\sum_{m=0}^{q}X^{\nu}_{qm}\bar{\Gamma}_{pmrs}+\sum_{m=q+1}^{M}Y^{\nu}_{mq}\bar{\Gamma}_{pmrs}\right)\\\
&-\left(\sum_{m=r+1}^{M}X^{\nu}_{mr}\bar{\Gamma}_{pqms}+\sum_{m=0}^{r}Y^{\nu}_{rm}\bar{\Gamma}_{pqms}\right)\\\
&+\left(\sum_{m=0}^{s}X^{\nu}_{sm}\bar{\Gamma}_{pqrm}+\sum_{m=s+1}^{M}Y^{\nu}_{sm}\bar{\Gamma}_{pqrm}\right),\end{split}$
(131)

where $M$ is the number of (natural) spatial orbitals. If we define

$\displaystyle Q_{pq}^{\nu}$ $\displaystyle=X_{pq}^{\nu}\ \forall\ p>q$ (132)
$\displaystyle Q_{pq}^{\nu}$ $\displaystyle=Y_{qp}^{\nu}\ \ \forall\ q>p$
(133) $\displaystyle R_{pq}^{\nu}$ $\displaystyle=Y_{pq}^{\nu}\ \ \forall\
p>q$ (134) $\displaystyle R_{pq}^{\nu}$ $\displaystyle=X_{qp}^{\nu}\ \forall\
q>p$ (135)

then we have

$\bar{\Gamma}_{pqrs}^{\nu}=\sum_{m}\left(-Q^{\nu}_{mp}\bar{\Gamma}_{mqrs}+R^{\nu}_{mq}\bar{\Gamma}_{pmrs}-Q^{\nu}_{mr}\bar{\Gamma}_{pqms}+R^{\nu}_{ms}\bar{\Gamma}_{pqrm}\right).$
(136)

A helpful reference point to benchmark is the Hartree–Fock limit
$|\Psi_{0}\rangle\rightarrow|\mathrm{RHF}\rangle$ then

$\displaystyle\bar{\Gamma}_{pqrs}$
$\displaystyle=\bar{\gamma}_{pq}\bar{\gamma}_{rs}-\frac{1}{2}\bar{\gamma}_{ps}\bar{\gamma}_{rq}$
(137)
$\displaystyle=n_{q}\delta_{pq}n_{s}\delta_{rs}-\frac{1}{2}n_{s}\delta_{ps}n_{q}\delta_{rq},$
(138)

where $n_{i}=2$ and $n_{a}=0$ for $i$ and occupied MO and $a$ a virtual MO.
Let us first insert the Coulomb like terms from Eq. 138 into the first two
terms of Eq. 136. We have

$\displaystyle\sum_{m}-Q_{mp}\bar{\Gamma}_{mqrs}+R_{mq}\bar{\Gamma}_{pmrs}$
$\displaystyle\rightarrow\sum_{m}(R_{mq}n_{p}\delta_{mp}n_{s}-Q_{mp}n_{q}\delta_{mq})n_{s}\delta_{rs}$
(139) $\displaystyle=(R_{pq}n_{p}-Q_{qp}n_{q})n_{s}\delta_{rs}.$ (140)

Now for $p>q$ we have

$\displaystyle\sum_{m}-Q_{mp}\bar{\Gamma}_{mqrs}+R_{mq}\bar{\Gamma}_{pmrs}$
$\displaystyle\rightarrow(Y_{pq}^{\nu}n_{p}-Y_{pq}n_{q})n_{s}\delta_{rs}$
(141) $\displaystyle=(n_{p}-n_{q})Y_{pq}n_{s}\delta_{rs}$ (142)

and for $p<q$ we have

$\displaystyle\sum_{m}-Q_{mp}\bar{\Gamma}_{mqrs}+R_{mq}\bar{\Gamma}_{pmrs}$
$\displaystyle\rightarrow(X_{qp}^{\nu}n_{p}-X_{qp}n_{q})n_{s}\delta_{rs}$
(143) $\displaystyle=(n_{p}-n_{q})X_{qp}n_{s}\delta_{rs}$ (144)
$\displaystyle=-(n_{p}-n_{q})X_{pq}n_{s}\delta_{rs}\ \forall\ p>q$ (145)

where we relabelled $p$ and $q$ in the last line. Thus we have

$\sum_{m}-Q_{mp}\bar{\Gamma}_{mqrs}+R_{mq}\bar{\Gamma}_{pmrs}\xrightarrow[]{\text{Coulomb-
like}}\gamma_{pq}^{\nu}n_{s}\delta_{rs},$ (146)

and by inspection

$\sum_{m}-Q_{mr}\bar{\Gamma}_{pqms}+R_{ms}\bar{\Gamma}_{pqrm}\xrightarrow[]{\text{Coulomb-
like}}\gamma_{rs}^{\nu}n_{q}\delta_{pq}.$ (147)

For the exchange terms we just need to set $q\rightarrow s$ so putting it
together we find

$\bar{\Gamma}^{\nu}_{pqrs}=\gamma_{pq}^{\nu}n_{s}\delta_{rs}+\gamma_{rs}^{\nu}n_{p}\delta_{pq}-\frac{1}{2}\gamma_{ps}^{\nu}n_{q}\delta_{rq}-\frac{1}{2}\gamma_{rq}^{\nu}n_{s}\delta_{ps}.$
(148)

This expression is useful for checking that exchange-induction and exchange-
dispersion terms reduce to their appropriate SAPT0 expressions.

###  C.4 Numerical Implementation of SAPT Expressions Within the ERPA

In principle one could proceed to evaluate the SAPT expressions in Appendix B
by building the transition one- and two-particle density matrices. Indeed,
this is possible for the first-order SAPT terms that do not contain any
response term, and the second-order polarization energies which are functions
of the transition one-particle reduced density matrix only. However, forming
the transition two-particle density matrix, which is required for exchange-
induction and exchange-induction is not practical in general so some work
needs to be done. We follow Refs. 56; 57 and try to optimize the contraction
order of the intermediate tensors. Let us first look at the exchange-induction
term. The first two terms only contain the one-particle (transition) density
matrix; we start with the derivation of the third term:

$\displaystyle F^{\mu}_{3a}$
$\displaystyle=-\sum_{m>p}\tilde{Q}^{\mu}_{mp}\bar{\Gamma}^{mp^{\prime}}_{p^{\prime\prime}p^{\prime\prime\prime}}\bar{\gamma}_{qq^{\prime}}S_{q^{\prime}p^{\prime}}\tilde{v}^{p^{\prime\prime}p^{\prime\prime\prime}}_{pq}.$
(149)

Let

$T^{mp^{\prime}}_{pq}=\bar{\Gamma}^{mp^{\prime}}_{p^{\prime\prime}p^{\prime\prime\prime}}\tilde{v}_{pq}^{p^{\prime\prime}p^{\prime\prime\prime}}$
(150)

and

$N_{p^{\prime}q}=\bar{\gamma}_{qq^{\prime}}S_{q^{\prime}p^{\prime}}$ (151)

then

$W_{pm}=T_{pq}^{mp^{\prime}}N_{p^{\prime}q}$ (152)

so

$F_{3a}^{\mu}=-\sum_{m>p}\tilde{Q}^{\mu}_{mp}W_{pm}.$ (153) $\displaystyle
F^{\mu}_{3b}$
$\displaystyle=\sum_{m>p^{\prime\prime}}\tilde{R}^{\mu}_{mp^{\prime\prime}}\bar{\Gamma}^{pp^{\prime}}_{mp^{\prime\prime\prime}}\bar{\gamma}_{qq^{\prime}}S_{q^{\prime}p^{\prime}}\tilde{v}^{p^{\prime\prime}p^{\prime\prime\prime}}_{pq}$
(154)
$\displaystyle=\sum_{m>p^{\prime\prime}}\tilde{R}^{\mu}_{mp^{\prime\prime}}\left(\bar{\Gamma}_{mp^{\prime\prime\prime}}^{pp^{\prime}}\tilde{v}_{pq}^{p^{\prime\prime}p^{\prime\prime\prime}}\right)\left(\bar{\gamma}_{qq^{\prime}}S_{q^{\prime}p^{\prime}}\right)$
(155)
$\displaystyle=\sum_{m>p^{\prime\prime}}\tilde{R}^{\mu}_{mp^{\prime\prime}}T_{mq}^{p^{\prime\prime}p^{\prime}}N_{p^{\prime}q}$
(156)
$\displaystyle=\sum_{m>p^{\prime\prime}}\tilde{R}^{\mu}_{mp^{\prime\prime}}W_{p^{\prime\prime}m}.$
(157)

Next we have

$\displaystyle F^{\mu}_{3c}$
$\displaystyle=-\sum_{m>p^{\prime}}\tilde{Q}^{\mu}_{mp^{\prime}}\bar{\Gamma}^{pm}_{p^{\prime\prime}p^{\prime\prime\prime}}\bar{\gamma}_{qq^{\prime}}S_{q^{\prime}p^{\prime}}\tilde{v}^{p^{\prime\prime}p^{\prime\prime\prime}}_{pq}$
(158)
$\displaystyle=-\sum_{m>p^{\prime}}\tilde{Q}^{\mu}_{mp^{\prime}}\left(\bar{\Gamma}_{p^{\prime\prime}p^{\prime\prime\prime}}^{pm}\tilde{v}_{pq}^{p^{\prime\prime}p^{\prime\prime\prime}}\right)\left(\bar{\gamma}_{qq^{\prime}}S_{q^{\prime}p^{\prime}}\right)$
(159)
$\displaystyle=-\sum_{m>p^{\prime}}\tilde{Q}^{\mu}_{mp^{\prime}}T_{mq}N_{p^{\prime}q}$
(160)
$\displaystyle=-\sum_{m>p^{\prime}}\tilde{Q}^{\mu}_{mp^{\prime}}W_{p^{\prime}m}$
(161)

and finally

$\displaystyle F^{\mu}_{3d}$
$\displaystyle=\sum_{m>p^{\prime\prime\prime}}\tilde{R}^{\mu}_{mp^{\prime\prime\prime}}\bar{\Gamma}^{pp^{\prime}}_{p^{\prime}m}\bar{\gamma}_{qq^{\prime}}S_{q^{\prime}p^{\prime}}\tilde{v}^{p^{\prime\prime}p^{\prime\prime\prime}}_{pq}$
(162)
$\displaystyle=\sum_{m>p^{\prime\prime\prime}}\tilde{R}^{\mu}_{mp^{\prime\prime\prime}}\left(\bar{\Gamma}_{p^{\prime\prime}m}^{pp^{\prime}}\tilde{v}_{pq}^{p^{\prime\prime}p^{\prime\prime\prime}}\right)\left(\bar{\gamma}_{qq^{\prime}}S_{q^{\prime}p^{\prime}}\right)$
(163)
$\displaystyle=\sum_{m>p^{\prime\prime\prime}}\tilde{R}^{\mu}_{mp^{\prime\prime\prime}}T_{mq}^{p^{\prime}p^{\prime\prime\prime}}N_{p^{\prime}q}$
(164)
$\displaystyle=\sum_{m>p^{\prime\prime\prime}}\tilde{R}^{\mu}_{mp^{\prime\prime\prime}}W_{p^{\prime\prime}m}.$
(165)

For the fourth term we have to evaluate something like

$F^{\mu}_{4a}=\left(\sum_{m>p}X_{mp}^{\nu}+\sum_{p>m}Y_{pm}^{\nu}\right)\bar{\Gamma}^{mp^{\prime}}_{p^{\prime\prime}p^{\prime\prime\prime}}\bar{\Gamma}_{q^{\prime\prime}q^{\prime\prime\prime}}^{qq^{\prime}}S_{p^{\prime}q^{\prime\prime\prime}}S_{q^{\prime}p^{\prime\prime\prime}}\tilde{v}^{p^{\prime\prime}q^{\prime\prime}}_{pq}$
(166)

we can define

$U_{pp^{\prime}}^{p^{\prime\prime}q^{\prime}}=\bar{\Gamma}_{q^{\prime\prime}q^{\prime\prime\prime}}^{qq^{\prime}}\tilde{v}_{pq}^{p^{\prime\prime}q^{\prime\prime}}S_{p^{\prime}q^{\prime\prime\prime}}$
(167)

and

$M_{p^{\prime\prime}q^{\prime}}^{mp^{\prime}}=\bar{\Gamma}_{p^{\prime\prime}p^{\prime\prime\prime}}^{mp^{\prime}}S_{q^{\prime}p^{\prime\prime\prime}}$
(168)

then

$Z_{pm}=U_{pp^{\prime}}^{p^{\prime\prime}q^{\prime}}M_{p^{\prime\prime}q^{\prime}}^{mp^{\prime}}$
(169)
$F^{\mu}_{4a}=\left(\sum_{m>p}X_{mp}^{\nu}+\sum_{p>m}Y_{pm}^{\nu}\right)Z_{pm}$
(170)

The exchange-dispersion part is much more verbose and let us look at the most
complex part which contains contributions like

$F^{\mu\nu}_{4}=\left(\sum_{m>p}X_{mp}^{\mu}+\sum_{p>m}Y_{pm}^{\mu}\right)\left(\sum_{n>q}X_{nq}^{\nu}+\sum_{p>n}Y_{pn}^{\nu}\right)\bar{\Gamma}^{mp^{\prime}}_{p^{\prime\prime}p^{\prime\prime\prime}}\bar{\Gamma}_{q^{\prime\prime}q^{\prime\prime\prime}}^{nq^{\prime}}S_{p^{\prime}q^{\prime\prime\prime}}S_{q^{\prime}p^{\prime\prime\prime}}\tilde{v}^{p^{\prime\prime}q^{\prime\prime}}_{pq}$
(171)

of which there are sixteen different terms arising from all combinations of
the $m$ and $n$ indices in the TPDMs. Let us go about this systematically, and
build some common intermediates (inspired by Hapka(Hapka  _et al._ , 2019)).
For convenience let us use upper case indices for monomer B and lower case
indices for monomer A and lower all indices. Given the finite number of
letters in the alphabet, intermediates will often have reuse symbols. For
safety we will have to assume they are only defined locally (i.e. only for one
of the 16 terms), although by inspection many of them are identical (up to
permutation of indices). Let us first at least define the protected
intermediates (similar to Hapka(Hapka  _et al._ , 2019)) which can be used in
practice:

$\displaystyle N_{mqrR}^{A}$ $\displaystyle=\Gamma_{mqrs}S_{sR}$ (172)
$\displaystyle N_{MQRr}^{B}$ $\displaystyle=\Gamma_{MQRS}S_{rS}$ (173)
$\displaystyle Z_{mpMP}^{AB}$
$\displaystyle=N_{mqrR}N_{MQRr}\tilde{v}_{pP}^{qQ}$ (174) $\displaystyle
U_{rsPQ}$ $\displaystyle=\Gamma_{pqrs}\tilde{v}_{pP}^{qQ}$ (175)
$\displaystyle W_{rsRS}^{AB}$
$\displaystyle=\Gamma_{pqrs}\Gamma_{PQRS}\tilde{v}_{pP}^{qQ}$ (176)
$\displaystyle=U_{rsPQ}\Gamma_{PQRS},$ (177)

with obvious analogues if monomer $B$’s indices are contracted over.

$\displaystyle\Gamma_{mqrs}\Gamma_{MQRS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{mqrR}N^{B}_{MQRr}\tilde{v}_{pqPQ}$ (178)
$\displaystyle=T^{AB}_{mqMQ}\tilde{v}_{pqPQ}$ (179)
$\displaystyle=V^{AB}_{mpMP}$ (180)

then

$\displaystyle F^{\mu\nu}_{4(1)}$
$\displaystyle=Q^{\mu}_{mp}Q^{\nu}_{MP}V^{AB}_{mpMP}$ (181)
$\displaystyle=G^{\mu}_{MP}Q^{\nu}_{MP}.$ (182)

Now there is some pattern in this, namely we will have contributions where the
contracted $Q/R$ matrix index ($m,M$) will be in either half of the 2PDM. Let
us proceed systematically along $\mu$ first then $\nu$. For the second term we
have:

$\displaystyle\Gamma_{pmrs}\Gamma_{MQRS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{pmrR}N^{B}_{MQRr}\tilde{v}_{pqPQ}$ (183)
$\displaystyle=T^{AB}_{pmMQ}\tilde{v}_{pqPQ}$ (184)
$\displaystyle=V^{AB}_{mqMP}$ (185) $\displaystyle F^{\mu\nu}_{4(2)}$
$\displaystyle=R^{\mu}_{mq}Q^{\nu}_{MP}V^{AB}_{mqMP}$ (186)
$\displaystyle=G^{\mu}_{MP}Q^{\nu}_{MP}.$ (187)

For the third term we have:

$\displaystyle\Gamma_{pqms}\Gamma_{MQRS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{pqmR}N^{B}_{MQRr}\tilde{v}_{pqPQ}$ (188)
$\displaystyle=T^{AB}_{mRPQ}N^{B}_{MQRr}$ (189) $\displaystyle=V^{AB}_{mrMP}\
\ \text{or}$ (190) $\displaystyle=U_{rsPQ}(S_{sR}N^{B}_{MQRr})$ (191)
$\displaystyle=V^{AB}_{mrMP}$ (192) $\displaystyle F^{\mu\nu}_{4(3)}$
$\displaystyle=Q^{\mu}_{mr}Q^{\nu}_{MP}V_{mrMP}$ (193)
$\displaystyle=G^{\mu}_{MP}Q^{\nu}_{MP}.$ (194)

For the fourth term we have:

$\displaystyle\Gamma_{pqrm}\Gamma_{MQRS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{pqSm}N^{B}_{MQsS}\tilde{v}_{pqPQ}$ (195)
$\displaystyle=T^{AB}_{mSPQ}N^{B}_{MQsS}$ (196) $\displaystyle=V^{AB}_{msMP}$
(197) $\displaystyle F^{\mu\nu}_{4(4)}$
$\displaystyle=R^{\mu}_{ms}Q^{\nu}_{MP}V^{AB}_{msMP}$ (198)
$\displaystyle=G^{\mu}_{MP}Q^{\nu}_{MP}.$ (199)

For the fifth term we have:

$\displaystyle\Gamma_{mqrs}\Gamma_{PMRS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{mqrR}N^{B}_{PMRr}\tilde{v}_{pqPQ}$ (200)
$\displaystyle=T^{AB}_{mqMP}\tilde{v}_{pqPQ}$ (201)
$\displaystyle=V^{AB}_{mpMQ}$ (202) $\displaystyle F^{\mu\nu}_{4(5)}$
$\displaystyle=Q^{\mu}_{mp}R^{\nu}_{MQ}V^{AB}_{mpMQ}$ (203)
$\displaystyle=G^{\mu}_{MQ}R^{\nu}_{MQ}.$ (204)

For the sixth term we have:

$\displaystyle\Gamma_{pmrs}\Gamma_{PMRS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{pmrR}N^{B}_{PMRr}\tilde{v}_{pqPQ}$ (205)
$\displaystyle=T^{AB}_{pmPM}\tilde{v}_{pqPQ}$ (206)
$\displaystyle=V^{AB}_{mqMQ}$ (207) $\displaystyle F^{\mu\nu}_{4(6)}$
$\displaystyle=R^{\mu}_{mq}R^{\nu}_{MQ}V^{AB}_{mqMQ}$ (208)
$\displaystyle=G^{\mu}_{MQ}R^{\nu}_{MQ}.$ (209)

For the seventh term we have:

$\displaystyle\Gamma_{pqms}\Gamma_{PMRS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{pqmR}N^{B}_{PMRr}\tilde{v}_{pqPQ}$ (210)
$\displaystyle=T^{AB}_{mRPQ}N^{B}_{PMRr}$ (211) $\displaystyle=V^{AB}_{mrMQ}$
(212) $\displaystyle F^{\mu\nu}_{4(7)}$
$\displaystyle=Q^{\mu}_{mr}R^{\nu}_{MQ}V_{mrMQ}$ (213)
$\displaystyle=G^{\mu}_{MQ}R^{\nu}_{MQ}.$ (214)

For the eight term we have:

$\displaystyle\Gamma_{pqrm}\Gamma_{PMRS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{pqSm}N^{B}_{PMsS}\tilde{v}_{pqPQ}$ (215)
$\displaystyle=T^{AB}_{mSPQ}N^{B}_{PMsS}$ (216) $\displaystyle=V^{AB}_{msMQ}$
(217) $\displaystyle F^{\mu\nu}_{4(8)}$
$\displaystyle=R^{\mu}_{ms}R^{\nu}_{MQ}V^{AB}_{msMQ}$ (218)
$\displaystyle=G^{\mu}_{MQ}Q^{\nu}_{MQ}.$ (219)

For the ninth term we have:

$\displaystyle\Gamma_{mqrs}\Gamma_{PQMS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{mqrR}N^{B}_{PQMr}\tilde{v}_{pqPQ}$ (220)
$\displaystyle=N^{A}_{mqrR}T^{AB}_{Mrpq}$ (221) $\displaystyle=V^{AB}_{mpMR}$
(222) $\displaystyle F^{\mu\nu}_{4(9)}$
$\displaystyle=Q^{\mu}_{mp}Q^{\nu}_{MR}V^{AB}_{mpMR}$ (223)
$\displaystyle=G^{\mu}_{MR}Q^{\nu}_{MR}.$ (224)

For the tenth term we have:

$\displaystyle\Gamma_{pmrs}\Gamma_{PQMS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{pmrR}N^{B}_{PQMr}\tilde{v}_{pqPQ}$ (225)
$\displaystyle=N^{A}_{pmrR}T^{AB}_{Mrpq}$ (226) $\displaystyle=V^{AB}_{mqMR}$
(227) $\displaystyle F^{\mu\nu}_{4(10)}$
$\displaystyle=R^{\mu}_{mq}Q^{\nu}_{MR}V^{AB}_{mqMR}$ (228)
$\displaystyle=G^{\mu}_{MR}Q^{\nu}_{MR}.$ (229)

For the eleventh term we have:

$\displaystyle\Gamma_{pqms}\Gamma_{PQMS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=T^{AB}_{msPQ}N^{B}_{PQMr}S_{sR}$ (230)
$\displaystyle=V^{AB}_{mrMR}$ (231) $\displaystyle F^{\mu\nu}_{4(11)}$
$\displaystyle=Q^{\mu}_{mr}Q^{\nu}_{MR}V^{AB}_{mrMR}$ (232)
$\displaystyle=G^{\mu}_{MR}Q^{\nu}_{MR}.$ (233)

For the twelfth term we have:

$\displaystyle\Gamma_{pqrm}\Gamma_{PQMS}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=T^{AB}_{mrPQ}N^{B}_{PQMr}S_{sR}$ (234)
$\displaystyle=V^{AB}_{msMR}$ (235) $\displaystyle F^{\mu\nu}_{4(12)}$
$\displaystyle=R^{\mu}_{ms}Q^{\nu}_{MR}V^{AB}_{msMR}$ (236)
$\displaystyle=G^{\mu}_{MR}Q^{\nu}_{MR}.$ (237)

For the thirteenth term we have:

$\displaystyle\Gamma_{mqrs}\Gamma_{PQRM}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{mqSR}T^{AB}_{pqRM}$ (238) $\displaystyle=V^{AB}_{mpMS}$
(239) $\displaystyle F^{\mu\nu}_{4(13)}$
$\displaystyle=Q^{\mu}_{mp}R^{\nu}_{MS}V^{AB}_{mpMS}$ (240)
$\displaystyle=G^{\mu}_{MS}Q^{\nu}_{MS}.$ (241)

For the fourteenth term we have:

$\displaystyle\Gamma_{pmrs}\Gamma_{PQRM}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=N^{A}_{mqSR}T^{AB}_{pqRM}$ (242) $\displaystyle=V^{AB}_{mqMS}$
(243) $\displaystyle F^{\mu\nu}_{4(14)}$
$\displaystyle=R^{\mu}_{mq}R^{\nu}_{MS}V^{AB}_{mqMS}$ (244)
$\displaystyle=G^{\mu}_{MS}Q^{\nu}_{MS}.$ (245)

For the fifteenth term we have:

$\displaystyle\Gamma_{pqms}\Gamma_{PQRM}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=\Gamma_{pqms}S_{rS}S_{sR}T^{AB}_{pqRM}$ (246)
$\displaystyle=U^{AB}_{msRM}S_{rS}S_{sR}$ (247) $\displaystyle=V^{AB}_{mrMS}$
(248) $\displaystyle F^{\mu\nu}_{4(15)}$
$\displaystyle=Q^{\mu}_{mr}R^{\nu}_{MS}V^{AB}_{mrMS}$ (249)
$\displaystyle=G^{\mu}_{MS}Q^{\nu}_{MS}.$ (250)

For the sixteenth term we have:

$\displaystyle\Gamma_{pqrm}\Gamma_{PQRM}S_{rS}S_{sR}\tilde{v}_{pqPQ}$
$\displaystyle=\Gamma_{pqrm}S_{rS}S_{sR}T^{AB}_{pqRM}$ (251)
$\displaystyle=U^{AB}_{rmRM}S_{rS}S_{sR}$ (252) $\displaystyle=V^{AB}_{msMS}$
(253) $\displaystyle F^{\mu\nu}_{4(16)}$
$\displaystyle=R^{\mu}_{ms}R^{\nu}_{MS}V^{AB}_{msMS}$ (254)
$\displaystyle=G^{\mu}_{MS}Q^{\nu}_{MS}.$ (255)

Naturally,

$F_{4}^{\mu\nu}=\sum_{i}F_{4(i)}^{\mu\nu}.$ (256)

##  Appendix D Methods (Computational Details)

All density functional theory calculations were performed with the pyscf
software packageSun  _et al._ (2020) and used a Lebedev grid used for the
numerical integration of the exchange correlation functionals with 70, 105,
140 radial grid points and 590 770, 770 angular grid points, for atoms of the
first, second and fourth period respectively (grid level 5). The geometry
optimizations were performed with pyscf and the geomeTRICWang and Song (2016)
package. All CASSCF calculation were performed with pyscf software packageSun
_et al._ (2020) and selected CI (semistochastic heat-bath configuration
interaction (SHCI)) calculations with the DICEHolmes  _et al._ (2016); Sharma
_et al._ (2017) program package via the pyscf plugin. We used a loose
screening parameter of 0.001 (sweep $\epsilon$) for computational efficiency
and in order to afford large active spaces. We used the 6-31GHehre  _et al._
(1972); Rassolov  _et al._ (1998) basis for the “stretched water” dimer and a
mixed basis (C: 6-31G; H: STO-3GHehre  _et al._ (1969); Pietro and Hehre
(1983)) for the benzene p-benzyne dimer. We also employed a mixed basis for
the [Mn(NH3)_3(CN)_2NO]^0$\cdots$X (X = HF,H2O,NH3 and CH4), where the crucial
moieties were described with the 6-31G basis (Mn, NO and X) and the ligand
framework (CN and NH3) STO-3G was used. The mixed basis was chosen due to the
memory constraints of the GPU for the subsequent ERPA calculations. This is
because our current implementation is limited in system size due to a non-
optimal treatment of core/active/virtual simpliciation and a lack of density
fitting of the response functions. Note that all VQE calculations were
performed on ideal statevector simulators, which were carried out on an GPU
accelerated in-house QC Ware package (quasar/vulcan). In addition, double
factorization was used for evaluating the total energy.Motta  _et al._ (2021);
Huggins  _et al._ (2021) All SAPT and VQE calculations were performed on the
newest NVIDIA A100 GPUs provided by Amazon Web Services.

The geometries of the water dimer and benzene were taken from the SI of our
previous study.Malone  _et al._ (2022). The [Mn(NH3)_3(CN)_2NO]^0 monomer was
optimized without constraints and the geometry was frozen for subsequent
geometry optimization of the hydrogen bonded dimer complexes. This ensures
that the active space orbitals remain as similar as possible in each complex.

The active spaces for the “stretched” water dimer and benzene p-benzyne dimer,
were selected based on MP2 natural orbitals followed by a CASSCF calculation
with (8e, 8o) and (6e, 6o), respectively. The active space of
[Mn(NH3)_3(CN)_2NO]^0 was selected via the automated construction of molecular
active spaces from atomic valence orbitals (AVAS)Sayfutyarova  _et al._
(2017), where we included all Mn 3-d and the p-orbitals of the NO ligand.
Additionally, we systematically added orbitals to the active space based on
SHCI resulting in a (16e, 22o) active space. The natural occupation numbers
are plotted in Fig. S2 and show several orbitals with strong deviation from
integer values confirming the multi-reference character of the system. The
subsequent CASSCF calculation used the SHCI natural orbitals as a starting
guess and included (6e, 6o) in the active space. The corresponding natural
orbitals are depicted in Fig. S5. The strongest deviation are seen in Fig. S5
(c)–(f); the orbitals show both bonding and antibonding $\pi$ type orbitals of
Mn-d_xz and Mn-d_yz with two NO $\pi^{*}$ orbitals. These orbitals correspond
to metal-to-ligand $\pi$ backbonding in the Dewar–Chatt–Duncanson pictureChatt
and Duncanson (1953); Chatt  _et al._ (1955). The strong static correlation of
these four electrons indicates that the electronic structure of this complex
is a superposition of the two configurations [Mn(II)-NO∙] and [Mn(II)-NO+].
The corresponding NOONs are depicted in Fig. S2 and the values only slightly
change in comparison to the (approximate) larger active space SHCI results.

In order to converge the $k$-mUCJ ansatz with larger repetition factors a
read-in algorithm inspired by Ref. 24 was implemented. For a given circuit
repetition number $N$, the optimal VQE parameters from the $N-1$^th step are
used as a starting guess and the remaining extra parameters are populated with
random parameters from a normal Gaussian distribution ($\mathcal{N}$) with
mean $\mu=0.0$ and variance $\sigma^{2}=0.001$ ($\mathcal{N}$(0.0, 0.001)). We
increased the repetition factor in the following order: $N=1,2,4,\dots,m$. All
calculations were performed in the dimer-centered basis (ghost atoms on the
other monomer) and a gradient threshold of $1\times 10^{-6}$ for the L-BFGS-B
solver in scipy or with a maximum iteration of 1500 whichever occurred first.
In cases when the maximum iteration was reached first we found the norm of the
gradient to be $\approx 1\times 10^{-5}$.

For the subsequent SAPT calculations we constructed the one- and two-particle
reduced density matrices in the full set of natural orbitals. In the case of
RHF these correspond to the canonical Hartree–Fock orbitals. For SAPT(CASSCF)
and SAPT(VQE) we used the CASSCF natural orbitals. These density matrices were
then used to construct the ERPA and SAPT equations given in Eq. 69 and Section
C.4. We used a threshold of 1 $\times 10^{-5}$ for the canonical
orthogonalization step when solving the ERPA generalized (symmetric)
eigenvalue problem. The ERPA and SAPT expressions given in in Eq. 69 and
Section C.4 were implemented in an in-house python code. NumPy(Van Der Walt
_et al._ , 2011; Harris  _et al._ , 2020) was used for tensor manipulation and
linear algebra operations with certain steps accelerated using CuPy(Nishino
and Loomis, 2017). This prototyping code does not exploit core-active-virtual
partitioning and is thus limited to systems containing roughly 130 orbitals.

##  Appendix E Additional Information for the Results Section

Table SI: Detailed SAPT energy contributions for all system studied in this manuscript from Fig. 2, 3 and 5 from the main manuscript (energies in kcal/mol, intermolecular distance $r$ in Å) System | Method | $E^{(1)}_{\mathrm{elst}}$ | $E^{(1)}_{\mathrm{exch}}$ | $E^{(2)}_{\mathrm{disp}}$ | $E^{(2)}_{\mathrm{{exch-disp}}}$ | $E^{(2)}_{\mathrm{ind,u}}$ | $E^{(2)}_{\mathrm{exch-ind,u}}$ | $E_{\mathrm{int}}$ | $r$
---|---|---|---|---|---|---|---|---|---
H2O…H2O | SAPT(CAS-CI) | -11.41 | 7.89 | -0.88 | 0.20 | -3.57 | 2.15 | -5.63 | 2.0
H2O…H2O | SAPT(VQE) ($k$=1) | -10.89 | 7.75 | -0.90 | 0.19 | -3.33 | 1.93 | -5.25 | 2.0
H2O…H2O | SAPT(VQE) ($k$=1) | -11.40 | 7.93 | -0.88 | 0.20 | -3.53 | 2.11 | -5.58 | 2.0
|  |  |  |  |  |  |  |  |
p-Bz…Bz | SAPT(CAS-CI) | 0.13 | 0.26 | -0.61 | 0.03 | -0.10 | 0.05 | -0.23 | 3.9
p-Bz…Bz | SAPT(VQE) ($k$=1) | 0.08 | 0.26 | -0.61 | 0.03 | -0.10 | 0.05 | -0.29 | 3.9
p-Bz…Bz | SAPT(VQE) ($k$=1) | 0.12 | 0.26 | -0.61 | 0.03 | -0.10 | 0.05 | -0.25 | 3.9
|  |  |  |  |  |  |  |  |
MnNO…HF | SAPT(CAS-CI) | -11.88 | 5.6 | -1.09 | 0.13 | -3.21 | 1.12 | -9.33 | 1.8
MnNO…HF | SAPT(VQE) ($k$=1) | -12.19 | 5.58 | -1.08 | 0.13 | -3.23 | 1.12 | -9.67 | 1.8
MnNO…HF | SAPT(VQE) ($k$=1) | -11.87 | 5.60 | -1.09 | 0.13 | -3.21 | 1.12 | -9.31 | 1.8
|  |  |  |  |  |  |  |  |
MnNO…HOH | SAPT(CAS-CI) | -7.59 | 3.62 | -0.86 | 0.14 | -1.60 | 0.74 | -5.55 | 2.1
MnNO…HOH | SAPT(VQE) ($k$=1) | -7.78 | 3.60 | -0.85 | 0.14 | -1.61 | 0.74 | -5.76 | 2.1
MnNO…HOH | SAPT(VQE) ($k$=1) | -7.6 | 3.62 | -0.85 | 0.14 | -1.6 | 0.74 | -5.56 | 2.1
|  |  |  |  |  |  |  |  |
MnNO…HNH2 | SAPT(CAS-CI) | -2.77 | 1.31 | -0.62 | 0.06 | -0.56 | 0.20 | -2.37 | 2.3
MnNO…HNH2 | SAPT(VQE) ($k$=1) | -2.83 | 1.31 | -0.61 | 0.06 | -0.57 | 0.19 | -2.45 | 2.3
MnNO…HNH2 | SAPT(VQE) ($k$=1) | -2.77 | 1.31 | -0.62 | 0.06 | -0.56 | 0.20 | -2.37 | 2.3
|  |  |  |  |  |  |  |  |
MnNO…HCH3 | SAPT(CAS-CI) | -0.10 | 0.02 | -0.13 | 0.00 | -0.07 | 0.00 | -0.27 | 3.5
MnNO…HCH3 | SAPT(VQE) ($k$=1) | -0.10 | 0.02 | -0.13 | 0.00 | -0.07 | 0.00 | -0.28 | 3.5
MnNO…HCH3 | SAPT(VQE) ($k$=1) | -0.10 | 0.02 | -0.13 | 0.00 | -0.07 | 0.00 | -0.27 | 3.5

###  E.1 Detailed SAPT Analysis of the MnNO Complexes

The detailed analysis of the energy term decomposition helps to unravel the
origin of the interaction. Fig. 5 (b) plots each component of the SAPT(VQE) (k
= 4) calculation of the series of hydrogen bonded complexes plus the
“stretched” water dimer as a reference of a typical hydrogen bond. We see that
the electrostatic term dominates the interaction as expected for hydrogen
bonds. The MnNO…HF complex has the largest electrostatic contribution ($-$11.9
kcal/mol) which is stronger than a water hydrogen bond ($-$11.4 kcal/mol). The
electrostatic interaction of MnNO…HOH is significantly weaker despite the fact
that the dipole moments of both H-F and H2O are around 1.8 D. When comparing
the exchange terms, we see that the terms rapidly decrease from MnNO…HF to
MnNO…HCH3 which can be rationalized by the change in bond distances and by the
exponential decay of this term. Interestingly, the water dimer and MnNO…HF
have similar bond distances but MnNO…HF exhibits a significantly less
repulsive term. This difference is the main driving force for the difference
in interaction energies and can be rationalized by the diffuseness of the lone
pairs. The bound NO becomes (partly) NO^+, which makes the lone pair more
compact in space than the lone pair in the water dimer; thus, resulting in
less exchange repulsion.

###  E.2 Additional Comments on the Comparison Between DFT and SAPT(VQE)

We compare the SAPT(VQE) interaction energies to DFT based supermolecular
(BSSE correctedBoys and Bernardi (1970b)) interaction energies. We note that
exact comparisons are difficult for two reasons: first, the static correlation
makes it difficult to generate reliable reference energies as the “gold
standard” CCSD(T) for non-covalent interactions is not reliable anymoreBulik
_et al._ (2015); second, SAPT0 produces most accurate energies with the jun-
cc-pVDZ basis setParker  _et al._ (2014) which cannot be employed due to
technical limits in the current implementation. However, nitrosyl complexes
are an example of non-universality problem of approximate density functionals
as the hydrogen bonding moiety and the nitrosyl moiety prefer different
approximate density functionalsRadon and Pierloot (2008); Radon  _et al._
(2010); Mardirossian and Head-Gordon (2014) and thus reliable prediction are
only possible with careful system specific benchmarking when experimental data
is available.Lehnert  _et al._ (2013) In contrast, SAPT is expected to give
accurate results for hydrogen bonds given that proper monomer wavefunctions
are used. The calculation of accurate interaction energies poses a challenge
for DFT as different types of approximate density functionals are recommended
for the different molecular moieties: non-hybrid functionals are recommended
for an accurate description of the metal-NO moietyRadon and Pierloot (2008);
Radon  _et al._ (2010), while for hydrogen bonds range-separated hybrids are
recommendedMardirossian and Head-Gordon (2017). These two competing
requirements make these systems very sensitive to the choice of specific
approximate exchange correlation functional illustrating the non-universality
problem of approximate density functionals. When using DFT, careful benchmark
against experimental data is necessary.Lehnert  _et al._ (2013) Our findings
indicate that the $k$-uCJ VQE ansatz is able to accurately describe the
difficult electronic structure of the nitrosyl monomer in the hydrogen bonding
complex (see Fig. S6) To illustrate this point, Fig. 6 (main text) plots the
interaction energies of SAPT(VQE), SAPT(CAS-CI) and several popular DFT
functionals (supramolecular). We included many popular functionals as well as
several top performing functionals for non-covalent interactionsMardirossian
and Head-Gordon (2017). We see in Fig. 6 (main text) that the SAPT(VQE) (k =
4) is almost identical to the SAPT(CAS-CI) in all four cases. We also note
that the DFT functionals exhibit a significant spread for each complex. The
B97-D functional predicted the smallest binding energy in all four cases, but
the highest interaction energy is predicted by a different functional.
Furthermore, we see the relative ordering of the functionals change for each
system (color sequence in each plot). This illustrates the non-universality
problem for approximate exchange correlation functionals even for very similar
nitrosyl complexes (this also holds true for larger basis set as illustrated
in Fig. S8). We note that the SAPT(CAS-CI) results are the reference for the
SAPT(VQE) calculations and do _not_ represent the true interaction energy,
thus, only the SAPT(VQE) not the DFT interaction energies should be compared
against this reference. B97-DGrimme (2006), BP86Perdew and Yue (1986); Becke
(1988), SCANSun  _et al._ (2015), PBEPerdew  _et al._ (1996), M06-LZhao and
Truhlar (2006), TPSShStaroverov  _et al._ (2003), B3LYPBecke (1988); Lee  _et
al._ (1988); Becke (1993), PBE0Adamo and Barone (1999), MN15Haoyu  _et al._
(2016), PWB6KZhao and Truhlar (2005), CAM-B3LYPYanai  _et al._ (2004),
$\omega$B97X-DChai and Head-Gordon (2008). This diverse list includes local
DFT functionals, hybrid functionals with a wide range of exact exchange
(10%–46%) and range separated hybrids.

Figure S2: Natural occupation number of the SHCI/CASSCF natural orbitals for
the strongly correlated monomer of each test case: upper panel “stretched” H2O
and p-benzyne; lower panel: [Mn(CN)2(NH3)3NO]^0.

(a)

(b)

Figure S3: (a) potential energy scan along $r$ using SAPT(CAS-CI), SAPT(RHF)
and SAPT(4-uCJ); (b) distance dependence of each SAPT(4-uCJ) energy term with
respect to the intermolecular distance r for the “stretched” water dimer.
Figure S4: Potential energy scan along $r$ using SAPT(CAS-CI), SAPT(RHF) and
SAPT(VQE) ($k$ = 4) for the p-benzyne-benzene dimer.

(a)

(b)

(c)

(d)

(e)

(f)

Figure S5: Natural orbitals of the (6e, 6o) CASSCF calculation of [Mn(CN)2(NH3)3NO]^0. Figure S6: Error in the one and two particle density matrix as a function of the circuit repetition factor $k$ for the [Mn(CN)2(NH3)3NO]^0 monomer. The errors are defined as $|\Delta\gamma|_{\infty}=|\gamma_{\mathrm{CASCI}}-\gamma_{k\mathrm{-muCJ}}|_{\infty}$ and $|\Delta\Gamma|_{\infty}=|\Gamma_{\mathrm{CASCI}}-\Gamma_{k\mathrm{-muCJ}}|_{\infty}$.  Table SII: Quantum hardware resource requirements with respect to the repetition factor $k$ (Number of qubits, number of two qubit gates (2-Q-G), number of parameter, quantum circuit depth) for the VQE simulations using the $k$-muCJ ansatz for the simulation of the [Mn(CN)2(NH3)3NO]^0 with an (6e, 6o) active space. k | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
# Qubits | 12 | 12 | 12 | 12 | 12 | 12 | 12 | 12
# 2-Q-G | 750 | 1065 | 1380 | 1695 | 2010 | 2325 | 2640 | 2955
# Param. | 45 | 75 | 105 | 135 | 165 | 195 | 225 | 255
Depth | 216 | 378 | 540 | 702 | 864 | 1026 | 1188 | 1350
Figure S7: Term-by-term decomposition of the SAPT(4-uCJ) binding energies of
each heme-nitrosyls hydrogen complex at a fixed bond distance of r(O-H) = 2.08
Å (equilibrium distance of the H2O complex). Figure S8: Supramolecuar DFT
interaction energies in the def2-TZVPDWeigend and Ahlrichs (2005) basis for
the HF, H2O and NH3 nitrosyl hydrogen complexes.

(a)

(b)

(c)

Figure S9: (a) Comparison between supramolecular DFT interaction energies for
the nitrosyl and nitric acid hydrogen complexes (small basis); (b) Comparison
between supramolecular DFT and SAPT in their most common basis sets using the
nitric acid hydrogen complexes. High level CCSD(T) reference interaction
energies are also provided for reference (DFT: def2-TZVPD & BSSE corrected;
SAPT: SAPT in the jun-cc-pVDZ basis (skips the VQE calculation and uses the
RHF wave function instead) and CCSD(T): aug-cc-pVQZ & BSSE corrected; for the
DFT functional color coding see (c)).

Table SIII: Detailed interaction energies from Fig. 6 (main text) using the
supramolecular approach for DFT and CCSD(T) calculation (BSSE corrected).
CCSD(T)/aug-cc-pVQZ can be considered very close to the true E_int using both
a very accurate wave function ansatz (for non strongly correlated systems) and
large basis set; DFT calculations using the def2-TZVPD basis set as
recommended by practical best practices from Ref. 133 and SAPT calculation use
the jun-cc-pVDZ basis set as recommended by Ref. 49.

System | Basis | B97-D | BP86 | SCAN | PBE | M06-L | TPSSh | B3LYP | PBE0 | MN15 | PWB6K | CAM-B3LYP |  $\omega$B97X-D
---|---|---|---|---|---|---|---|---|---|---|---|---|---
MnNO…HF | mixed | -6.58 | -7.06 | -8.37 | -7.88 | -7.99 | -7.66 | -7.91 | -8.36 | -8.89 | -8.51 | -8.94 | -7.84
MnNO…HOH | mixed | -3.44 | -3.88 | -5.23 | -4.63 | -4.83 | -4.33 | -4.45 | -4.87 | -5.35 | -5.11 | -5.26 | -4.45
HNO…HF | 6-31G | -4.82 | -5.83 | -7.5 | -6.81 | -5.7 | -6.44 | -6.38 | -7.07 | -6.22 | -6.75 | -6.9 | -5.76
H2O…H2O | 6-31G | -1.95 | -2.43 | -3.44 | -3.28 | -2.81 | -2.9 | -2.98 | -3.34 | -3.31 | -3.62 | -3.68 | -2.81
|  |  |  |  |  |  |  |  |  |  |  |  |
HNO…HF | def2-TZVPD | -4.23 | -5.06 | -6.32 | -5.88 | -5.57 | -5.39 | -5.45 | -5.79 | -5.59 | -5.92 | -6.22 | -5.22
HNO…H2O | def2-TZVPD | -1.24 | -1.75 | -2.84 | -2.57 | -2.36 | -2.09 | -2.17 | -2.49 | -2.49 | -2.66 | -2.75 | -2.02

System | CCSD(T) (Basis) | SAPT(RHF) (Basis)
---|---|---
HNO…HF | -5.5 (aug-cc-pVQZ) | -5.1 (jun-cc-pVDZ)
HNO…H2O | -2.6 (aug-cc-pVQZ) | -2.4 (jun-cc-pVDZ)

## References

  * Cao  _et al._ (2018) Y. Cao, J. Romero, and A. Aspuru-Guzik, IBM J. Res. Dev. 62, 6 (2018).
  * Bauer  _et al._ (2020) B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, Chem. Rev. 120, 12685 (2020).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Reiher  _et al._ (2017) M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, Proc. Natl. Acad. Sci. U.S.A. 114, 7555 (2017).
  * Li  _et al._ (2019) Z. Li, S. Guo, Q. Sun, and G. K. Chan, Nat. Chem. 11, 1026 (2019).
  * White and Martin (1999) S. R. White and R. L. Martin, J. Chem. Phys. 110, 4127 (1999).
  * Chan and Sharma (2011) G. K.-L. Chan and S. Sharma, Annu. Rev. Phys. Chem. 62, 465 (2011).
  * Lüchow (2011) A. Lüchow, Wiley Interdiscip. Rev. Comput. Mol. Sci. 1, 388 (2011).
  * Garniron  _et al._ (2018) Y. Garniron, A. Scemama, E. Giner, M. Caffarel, and P.-F. Loos, J. Chem. Phys. 149, 064103 (2018).
  * Levine  _et al._ (2020) D. S. Levine, D. Hait, N. M. Tubman, S. Lehtola, K. B. Whaley, and M. Head-Gordon, J. Chem. Theory Comput. 16, 2340 (2020).
  * Holmes  _et al._ (2016) A. A. Holmes, N. M. Tubman, and C. Umrigar, J. Chem. Theory Comput. 12, 3674 (2016).
  * Sharma  _et al._ (2017) S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi, and C. J. Umrigar, J. Chem. Theory Comput. 13, 1595 (2017).
  * Aspuru-Guzik  _et al._ (2005) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Science 309, 1704 (2005).
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020).
  * Cao  _et al._ (2019) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik, Chem. Rev. 119, 10856 (2019).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Cohn  _et al._ (2021) J. Cohn, M. Motta, and R. M. Parrish, PRX Quantum 2, 040352 (2021).
  * Bharti  _et al._ (2022) K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Rev. Mod. Phys. 94, 015004 (2022).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nat. Commun. 2014 51 5, 1 (2014).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. 18, 023023 (2016).
  * Parrish and McMahon (2019) R. M. Parrish and P. L. McMahon, arXiv preprint arXiv:1909.08925 (2019).
  * Stair  _et al._ (2020) N. H. Stair, R. Huang, and F. A. Evangelista, J. Chem. Theory Comput. 16, 2236 (2020).
  * Huggins  _et al._ (2020) W. J. Huggins, J. Lee, U. Baek, B. O’Gorman, and K. B. Whaley, New J. Phys. 22, 073009 (2020).
  * Klymko  _et al._ (2022) K. Klymko, C. Mejuto-Zaera, S. J. Cotton, F. Wudarski, M. Urbanek, D. Hait, M. Head-Gordon, K. B. Whaley, J. Moussa, N. Wiebe, _et al._ , PRX Quantum 3, 020323 (2022).
  * Huggins  _et al._ (2022) W. J. Huggins, B. A. O’Gorman, N. C. Rubin, D. R. Reichman, R. Babbush, and J. Lee, Nature , 1 (2022).
  * Elfving  _et al._ (2020) V. E. Elfving, B. W. Broer, M. Webber, J. Gavartin, M. D. Halls, K. P. Lorton, and A. Bochevarov, arXiv preprint arXiv:2009.12472 (2020).
  * Liu  _et al._ (2022) H. Liu, G. H. Low, D. S. Steiger, T. Häner, M. Reiher, and M. Troyer, Mater. theory 6, 11 (2022).
  * Kirsopp  _et al._ (2021) J. J. M. Kirsopp, C. Di Paola, D. Z. Manrique, M. Krompiec, G. Greene-Diniz, W. Guba, A. Meyder, D. Wolf, M. Strahm, and D. M. Ramo, arXiv preprint arXiv:2110.08163 (2021).
  * Rice  _et al._ (2021) J. E. Rice, T. P. Gujarati, M. Motta, T. Y. Takeshita, E. Lee, J. A. Latone, and J. M. Garcia, J. Chem. Phys. 154, 134115 (2021).
  * Kim  _et al._ (2022) I. H. Kim, Y.-H. Liu, S. Pallister, W. Pol, S. Roberts, and E. Lee, Phys. Rev. Res. 4, 023019 (2022).
  * Greene-Diniz  _et al._ (2022) G. Greene-Diniz, D. Z. Manrique, W. Sennane, Y. Magnin, E. Shishenina, P. Cordier, P. Llewellyn, M. Krompiec, M. J. Rančić, and D. M. Ramo, arXiv preprint arXiv:2203.15546 (2022).
  * Heifetz (2020) A. Heifetz,  _Quantum mechanics in drug discovery_ (Springer, New York, 2020).
  * O’Brien  _et al._ (2021) T. E. O’Brien, M. Streif, N. C. Rubin, R. Santagati, Y. Su, W. J. Huggins, J. J. Goings, N. Moll, E. Kyoseva, M. Degroote, _et al._ , arXiv preprint arXiv:2111.12437 (2021).
  * Goings  _et al._ (2022) J. J. Goings, A. White, J. Lee, C. S. Tautermann, M. Degroote, C. Gidney, T. Shiozaki, R. Babbush, and N. C. Rubin, arXiv preprint arXiv:2202.01244 (2022).
  * Malone  _et al._ (2022) F. D. Malone, R. M. Parrish, A. R. Welden, T. Fox, M. Degroote, E. Kyoseva, N. Moll, R. Santagati, and M. Streif, Chem. Sci. (2022).
  * Rezac and Hobza (2016) J. Rezac and P. Hobza, Chemical reviews 116, 5038 (2016).
  * Kodrycka and Patkowski (2019) M. Kodrycka and K. Patkowski, J. Chem. Phys. 151, 070901 (2019).
  * Boys and Bernardi (1970a) S. Boys and F. Bernardi, Mol. Phys. 19, 553 (1970a).
  * McClean  _et al._ (2018) J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nat. Commun. 9, 4812 (2018).
  * Jeziorski  _et al._ (1976) B. Jeziorski, M. Bulski, and L. Piela, Int. J. Quantum Chem. 10, 281 (1976).
  * Jeziorski  _et al._ (1994) B. Jeziorski, R. Moszynski, and K. Szalewicz, Chem. Rev. 94, 1887 (1994).
  * Patkowski (2020) K. Patkowski, Wiley Interdiscip. Rev. Comput. Mol. Sci. 10, e1452 (2020).
  * Korona  _et al._ (1997) T. Korona, R. Moszynski, and B. Jeziorski, in  _Advances in quantum chemistry_ , Vol. 28 (Elsevier, United States, 1997) pp. 171–188.
  * Hapka  _et al._ (2021a) M. Hapka, M. Przybytek, and K. Pernal, J. Chem. Theory Comput. 17, 5538 (2021a).
  * Casida  _et al._ (1998) M. E. Casida, C. Jamorski, K. C. Casida, and D. R. Salahub, J. Chem. Phys. 108, 4439 (1998).
  * McClean  _et al._ (2017) J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, and W. A. De Jong, Phys. Rev. A 95, 042308 (2017).
  * Korona (2013) T. Korona, Mol. Phys. 111, 3705 (2013).
  * Parker  _et al._ (2014) T. M. Parker, L. A. Burns, R. M. Parrish, A. G. Ryno, and C. D. Sherrill, J. Chem. Phys. 140, 094106 (2014).
  * Korona (2008a) T. Korona, J. Chem. Phys. 128, 224104 (2008a).
  * Korona (2008b) T. Korona, Phys. Chem. Chem. Phys. 10, 6509 (2008b).
  * Korona and Jeziorski (2008) T. Korona and B. Jeziorski, J. Chem. Phys. 128, 144107 (2008).
  * Korona (2009a) T. Korona, J. Chem. Theory Comput. 5, 2663 (2009a).
  * Chatterjee and Pernal (2012) K. Chatterjee and K. Pernal, J. Chem. Phys. 137, 204109 (2012).
  * Hapka  _et al._ (2018) M. Hapka, M. Przybytek, and K. Pernal, J. Chem. Theory Comput. 15, 1016 (2018).
  * Hapka  _et al._ (2019) M. Hapka, M. Przybytek, and K. Pernal, J. Chem. Theory Comput. 15, 6712 (2019).
  * Hapka  _et al._ (2021b) M. Hapka, M. Przybytek, and K. Pernal, J. Chem. Theory Comput. 17, 5538 (2021b).
  * Parrish  _et al._ (2019) R. M. Parrish, E. G. Hohenstein, P. L. McMahon, and T. J. Martínez, Phys. Rev. Lett. 122, 230401 (2019).
  * Cai  _et al._ (2020) X. Cai, W.-H. Fang, H. Fan, and Z. Li, Phys. Rev. Res. 2, 033324 (2020).
  * Ollitrault  _et al._ (2020) P. J. Ollitrault, A. Kandala, C.-F. Chen, P. K. Barkoutsos, A. Mezzacapo, M. Pistoia, S. Sheldon, S. Woerner, J. M. Gambetta, and I. Tavernelli, Phys. Rev. Res. 2, 043140 (2020).
  * Matsuzawa and Kurashige (2020) Y. Matsuzawa and Y. Kurashige, J. Chem. Theory Comput. 16, 944 (2020).
  * Maia and Moura (2014) L. B. Maia and J. J. Moura, Chem. Rev. 114, 97 (2014).
  * Moncada (1991) S. Moncada, Pharmacol rev 43, 109 (1991).
  * Crawford  _et al._ (2001) T. D. Crawford, E. Kraka, J. F. Stanton, and D. Cremer, J. Chem. Phys. 114, 10638 (2001).
  * Bruckdorfer (2005) R. Bruckdorfer, Mol. Aspects Med. 26, 3 (2005).
  * Lehnert  _et al._ (2018) N. Lehnert, H. T. Dong, J. B. Harland, A. P. Hunt, and C. J. White, Nat. Chem. Rev. 2, 278 (2018).
  * Kadish  _et al._ (2000) K. Kadish, K. M. Smith, and R. Guilard,  _The Porphyrin Handbook, Volume 3_ , Vol. 3 (Elsevier, San Diego, 2000).
  * Lehnert  _et al._ (2013) N. Lehnert, W. R. Scheidt, and M. W. Wolf, in  _Nitrosyl Complexes in Inorganic Chemistry, Biochemistry and Medicine II_ (Springer, United States, 2013) pp. 155–223.
  * Hunt and Lehnert (2015) A. P. Hunt and N. Lehnert, Acc. Chem. Res. 48, 2117 (2015).
  * Santolini (2011) J. Santolini, J. Inorg. Biochem. 105, 127 (2011).
  * Walker (2005) F. A. Walker, J. Inorg. Biochem. 99, 216 (2005).
  * Hayton  _et al._ (2002) T. W. Hayton, P. Legzdins, and W. B. Sharp, Chem. Rev. 102, 935 (2002).
  * Radon and Pierloot (2008) M. Radon and K. Pierloot, J. Phys. Chem. A 112, 11824 (2008).
  * Lehnert  _et al._ (2021) N. Lehnert, E. Kim, H. T. Dong, J. B. Harland, A. P. Hunt, E. C. Manickas, K. M. Oakley, J. Pham, G. C. Reed, and V. S. Alfaro, Chem. Rev. 121, 14682 (2021).
  * Tfouni  _et al._ (2010) E. Tfouni, F. G Doro, L. E Figueiredo, J. CM Pereira, G. Metzke, and D. W Franco, Curr. Med. Chem. 17, 3643 (2010).
  * Hickok and Thomas (2010) J. R. Hickok and D. D. Thomas, Curr. Pharm. Des. 16, 381 (2010).
  * Serafim  _et al._ (2012) R. Serafim, M. Primi, G. Trossini, and E. Ferreira, Curr. Med. Chem. 19, 386 (2012).
  * Ampßler  _et al._ (2020) T. Ampßler, G. Monsch, J. Popp, T. Riggenmann, P. Salvador, D. Schröder, and P. Klüfers, Angew. Chem. Int. Ed. 59, 12381 (2020).
  * Boguslawski  _et al._ (2011) K. Boguslawski, C. R. Jacob, and M. Reiher, J. Chem. Theory Comput. 7, 2740 (2011).
  * Siegbahn  _et al._ (2010) P. E. Siegbahn, M. R. Blomberg, and S.-L. Chen, J. Chem. Theory Comput. 6, 2040 (2010).
  * Goodrich and Lehnert (2013) L. E. Goodrich and N. Lehnert, J. Inorg. Biochem. 118, 179 (2013).
  * Radon  _et al._ (2010) M. Radon, E. Broclawik, and K. Pierloot, J. Phys. Chem. B 114, 1518 (2010).
  * Boys and Bernardi (1970b) S. F. Boys and F. Bernardi, Mol. Phys. 19, 553 (1970b).
  * Mardirossian and Head-Gordon (2014) N. Mardirossian and M. Head-Gordon, Phys. Chem. Chem. Phys. 16, 9904 (2014).
  * Mardirossian and Head-Gordon (2017) N. Mardirossian and M. Head-Gordon, Mol. Phys. 115, 2315 (2017).
  * Motta  _et al._ (2021) M. Motta, E. Ye, J. R. McClean, Z. Li, A. J. Minnich, R. Babbush, and G. K. Chan, npj Quantum Information 7, 1 (2021).
  * Huggins  _et al._ (2021) W. J. Huggins, J. R. McClean, N. C. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley, and R. Babbush, npj Quantum Information 7, 1 (2021).
  * Moszynski  _et al._ (1994) R. Moszynski, B. Jeziorski, S. Rybak, K. Szalewicz, and H. L. Williams, J. Chem. Phys. 100, 5080 (1994).
  * Korona (2009b) T. Korona, J. Chem. Theory Comput. 5, 2663 (2009b).
  * Veryazov  _et al._ (2011) V. Veryazov, P. A. Malmqvist, and B. O. Roos, Int. J. Quantum Chem. 111, 3329 (2011).