Accurate non-covalent interaction energies on noisy intermediate-scale quantum computers via second-order symmetry-adapted perturbation theory

Abstract

The calculation of non-covalent interaction energies on noisy intermediate-scale quantum (NISQ) computers appears to be challenging with straightforward application of existing quantum algorithms. For example, the 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 (Chem. Sci., 2022, 13, 3094), this provides a recipe for complete SAPT(VQE) interaction energies up to second order, which is a well established truncation. 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 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 a quantum computer, simulated through ideal statevectors. The errors of the total interaction energy are orders of magnitude lower than the corresponding VQE total energy errors of the monomer wavefunctions. In addition, we present heme-nitrosyl model complexes as a system class for near term quantum computing simulations. They are strongly correlated, biologically relevant and difficult to simulate with classical quantum chemical methods. This is illustrated with density functional theory (DFT) as the predicted interaction energies exhibit a strong sensitivity with respect to the choice of functional. Thus, this work paves the way to obtain accurate interaction energies on a NISQ-era quantum computer with few quantum resources. It is the first step in alleviating 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.

---

1. Introduction

Quantum computing has emerged as a promising platform to approach classically challenging problems in chemistry.^1,2 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^3,4 (KS-DFT) may fail to describe the electronic structure correctly; popular examples are Fe–S clusters or the FeMo-cofactor.^5,6

Classically, a proper treatment of these strongly correlated systems is achieved with multi-reference methods where the naive 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.^7–13 Alternatively, quantum algorithms^2,14–16 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¹⁷ suffers from relatively poor gate fidelity and a low qubit count,¹⁸ which poses two key challenges. First, it is important for NISQ-tailored quantum algorithms¹⁹ to minimize quantum resources. The most prominent NISQ methods are hybrid quantum-classical algorithms like the variational quantum eigensolver (VQE),^20,21 quantum Krylov methods,^18,22–26 or the fermionic quantum Monte Carlo method.²⁷ The second challenge is to find specific applications that could harness quantum computing.^28,29 Many application studies in chemistry use either reduced model systems or molecules with a simple electronic structure.^30–33 There are several promising application areas for quantum chemistry in computer aided drug design³⁴ namely, exploring potential energy surfaces,³⁵ simulating metalloenzymes³⁶ and computing protein-ligand interaction energies,³⁷ the last of which we consider in this work.

The computation of non-covalent interaction energies is a routine task in classical quantum chemistry^38,39 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.⁴⁰ 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 set of parameters, where it becomes increasingly difficult to reach the global minimum on the parameter surface.⁴¹ Second, accounting for the basis set superposition error (BSSE) in the supermolecular approach⁴⁰ 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 the VQE.^20,21

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)^42–44 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,³⁷ 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.⁴³ 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 authors³⁰ 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. The combination of first order energies (from our previous work³⁷) and second order energies (presented in this work) is an accurate truncation of the SAPT perturbation series.⁴³ Those methods usually provide interaction energies with an accuracy below 1 kcal mol⁻¹ for a large spectrum of non-covalent binding motifs, even with modest basis sets.^45,46 Therefore, this work provides a pathway to obtain accurate intermolecular interaction energies using wavefunctions from quantum computations (in this work using VQE wavefunctions). The improved SAPT(VQE) workflow is depicted in Fig. 1 (using a water dimer as an example).

Fig. 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 the VQE. The converged wavefunction yields 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,³⁷ and the latter two are presented in this study and employ 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).

Our approach follows the SAPT(FCI) approach of Korona⁴⁷ 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.,⁴⁸ with the 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⁴⁹ in TD-DFT or the quantum subspace expansion (QSE) for VQE excited states.⁵⁰ 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, 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² 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 the ESI† 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 overheads (though optimization techniques such as hole/active/particle separations and density fitting might significantly reduce these overheads). Complexity notwithstanding, the approach presented provides a recipe for a SAPT(VQE) doppelganger of SAPT(FCI) complete through all second order terms.

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 ESI (Sections 1–3).† We then demonstrate the numerical performance of all four SAPT(VQE) terms and total interaction energies for several small multi-reference dimers and a model heme-nitrosyl hydrogen bonding complex, using classical ideal statevector simulators to emulate the VQE.

2. 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_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

Ĥ = Ĥ_A + Ĥ_B + V̂ (2)

where we assume Ĥ_X|Ψ_X〉 = E_X|Ψ_X〉, where |Ψ_X〉 is the ground state wavefunction of monomer X and V̂ contains only the coulombic interactions between monomers 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 havewhere E_pol⁽ⁿ⁾ and E_exch⁽ⁿ⁾ are nth-order polarization and exchange energies respectively.

In our previous work,³⁷ we derived and implemented the first order polarisation energy, usually denoted as the electrostatic interaction energy E⁽¹⁾_elst and the first order exchange energy E⁽¹⁾_exch.

In the following, we provide a brief overview over the second order terms (the detailed equations are presented in ESI Section 1).† For the second order polarisation energy, it is conventional in SAPT to split it into induction and dispersion energy, i.e., E⁽²⁾_pol = E⁽²⁾_ind + E⁽²⁾_disp.

The second order induction energy is E⁽²⁾_ind defined as:

E⁽²⁾_ind = E_ind(A ← B) + E_ind(B ← A) (4)

where E_ind(A ← B) is the induction of monomer A in the presence of monomer B and defined as:

E⁽²⁾_ind(A ← B) = 〈Ψ⁰_AΨ⁰_B|V̂|Ψ^ind_AΨ⁰_B〉 (5)

and the first order induction wave function for monomer A|Ψ^ind_A〉 is defined as:where |Ψ_A^μ〉 corresponds to the μth excited state wavefunction and E_A^μ to the μth excited state energy of monomer A and [capital Omega, Greek, circumflex]^B is the effective electrostatic potential of monomer B (see ESI 2† for an exact definition). It follows thatand E_ind(B ← A) follows analogously.

The dispersion energy is defined as:

E⁽²⁾_disp = 〈Ψ⁰_AΨ⁰_B|V̂|Ψ^disp_AB〉 (8)

where the first order dispersion wave function is defined as:yielding:where V̂_ee is the usual electron–electron Coulomb operator.

If one neglects all electron exchange processes other than those that exchange a single electron between monomers A and B one arrives at the so-called S² approximation to the exchange energies,^51–54 and the second order exchange terms are defined as:

E⁽²⁾_exch-ind(A → B) = 〈Ψ⁰_AΨ⁰_B|(V̂ − V̄)(P̂ − P̄)|Ψ^ind_AΨ⁰_B〉 (11)

E⁽²⁾_exch-disp = 〈Ψ⁰_AΨ⁰_B|(V̂ − V̄)(P̂ − P̄)|Ψ^disp_AB〉 (12)

where P̂ is defined as the single electron exchange operator, P̄ is defined as 〈Ψ⁰_AΨ⁰_B|P̂|Ψ⁰_AΨ⁰_B〉 and V̄ is defined as E⁽¹⁾_pol.

The combination of both first and second order terms yields the following expression for the interaction energy:

E_int ≈ E_SAPT = E⁽¹⁾_elst + E⁽¹⁾_exch + E⁽²⁾_ind + E⁽²⁾_disp + E⁽²⁾_exch-disp + E⁽²⁾_exch-ind (13)

Both the dispersion and induction energies are often combined with their exchange counterparts and the order superscripts are dropped to yield a simplified energy expression:

E_SAPT = E_elst + E_exch + E_ind + E_disp (14)

(for further details on the derivations see ESI Section 1).†

If we use high quality monomer wavefunctions approaching full configuration interaction (FCI), we conceptually approach 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 static (induction) and frequency dependent dynamic susceptibility tensors (dispersion), followed by contractions of these property tensors to form polarization and exchange SAPT contributions.^{47,52,55–57} In this work, we instead employ extended random phase approximation (ERPA)⁵⁸ as pioneered within SAPT(CASSCF) by Hapka^48,53,54,59 to avoid computing excited state properties on a quantum computer as they often require a significant measurement overhead on NISQ-era quantum computers.^22,60–62 We note that the use of ERPA as a proxy for the explicit response properties of the VQE resembles the quantum subspace expansion (QSE) method,⁵⁰ 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 ESI† along with explicit expressions for each term in terms of ERPA transition density matrices (ESI Sections 2 and 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.⁴⁶ Third- and higher-order extensions are likely possible along similar response properties 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 calculated on a quantum computer (see ESI 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:

|Ψ_VQE〉 ≡ Û_VQE|Φ_I〉 (15)

where |Φ_I〉 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⁶³ (VQE) which takes the formwhere K̂^(k) and 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. 63, which we denote as k-muCJ for clarity, with the ‘m’ standing for modified (see ESI 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.

3. Results

As a first step, we test SAPT(VQE) with two classic intermolecular interaction motifs (water dimer and t-shaped benzene dimer). As a second step, we apply SAPT(VQE) to a heme-nitrosyl model complex; these systems are highly relevant in both biological⁶⁴ and pharmaceutical⁶⁵ chemistry. The SAPT(VQE) results presented in this section are the result of ideal statevector VQE simulations (see ESI Section 4† 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 Fig. 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.

Fig. 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, and 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.

3.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,³⁷ 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. 2a. 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 and 8o) in the active space (for a detailed procedure on how the active orbitals were selected for the CASSCF calculations see ESI Section 4).† 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 from the large deviation of the absolute energy of the monomer (see Fig. 2b 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⁻¹. This can lead to large resource requirements for the supermolecular ansatz: in a worst case scenario for the supermolecular ansatz, both monomers are treated with an active space which then increases the resource requirement for the dimer computation even further (assuming that the active spaces combine). This is illustrated for the water dimer, we chose the same (8e, 8o) active space for the “stretched” water and add a (4e, 4o) active space for the second water yielding a (12e, 12o) and 24 qubits for the dimer. The VQE calculation for the dimer system did not converge sufficiently even with k > 20 (ε ∼30 kcal mol⁻¹). This illustrates the difficulties of using the supermolecular ansatz in the NISQ era.

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 the next step, we probed the bond dissociation of the water complex (along the intermolecular O–H bond labeled r in Fig. 2a). 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. 3a. The p-benzyne monomer has a biradical ground state, which is difficult to describe with classical single reference methods.⁶⁶ The k-muCJ VQE ansatz requires k = 8 for a moderate (6e, 6o) active space to reach sub 1 kcal mol⁻¹ accuracy as illustrated in Fig. 3b. The key findings are identical to the previous test case and also hold for different intermolecular distances (see Fig. 3b and c). Thus, these findings hold for very different types of intermolecular interaction motifs, intermolecular distances, different active spaces and different type difficult electronic structures.

Fig. 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, and the dotted gray line represents the chemical accuracy 1 kcal mol⁻¹ threshold); (c) absolute error 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.

3.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^65,67 or as a key intermediate in the nitrogen cycle.^64,68 At the center of these processes are metalloporphyrins,^69–71 where NO binds to the metal center as a nitrosyl ligand.^72,73

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⁷⁴ as illustrated by theoretical,⁷⁵ experimental⁷⁶ and medicinal^77–79 work. Unfortunately, the metal–NO bond in nitrosyl complexes poses a challenge for many quantum chemistry methods due to the redox active nature of the NO ligand.⁸⁰ There are three possible oxidation states for the NO moiety: NO⁻, NO˙ and NO⁺, which is illustrated for a generic M–NO complex in Fig. 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.⁸¹ This results in a large variety of recommended functionals depending on the specific nitrosyl complex studied.^70,75,82,83

Fig. 4 Upper left panel: different oxidation states of a transition metal bound NO; bottom left panel: structure of the [Mn(CN)₂(NH₃)₃NO]⁰ complex and the natural orbitals (NOON) with the largest deviation from integer value occupation (π type metal to ligand backbonding); right panel: simplification of the heme ligand framework and schematic representation of the four manganese nitrosyl hydrogen bond complexes.

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 the ligand field, total charge and π-acidity (see Fig. 4). As hydrogen bond donors we chose HF, H₂O, NH₃ and CH₄, 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⋯HNH₂ and MnNO⋯HCH₃.

As a first step, we analyzed the electronic structure of [Mn(CN)₂(NH₃)₃NO]⁰ and found the system to be strongly correlated (see ESI Section 5† for a more detailed discussion). We used a (6e, 6o) active space for the 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⁸⁴ (see Fig. S2 and Section 5†). The k-uCJ VQE 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. 5a). 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. However, we reiterate that this work uses an ideal statevector simulator; thus, the effect of noise is yet to be seen.

Fig. 5 (a) Top figure: absolute errors of the VQE total energy for the [Mn(CN)₂(NH₃)₃NO]⁰ 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, and 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 Table S1† for details).

In addition to the interaction energy, SAPT provides a decomposition into physical meaningful terms helping to unravel the origin of the interaction. Fig. 6b 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 ESI 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 the 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.

Fig. 6 Interaction energies using several approximate DFT exchange correlation functionals (supermolecular 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).

Finally, we compare the SAPT(VQE) interaction energy to DFT based supermolecular (BSSE corrected⁸⁵) 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^75,84,86 and thus reliable predictions are only possible with careful system specific benchmarking when experimental data are available.⁷⁰ In contrast, SAPT is expected to robustly give accurate results for hydrogen bonds given appropriate 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.⁸⁷ 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 functional 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 a 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, and 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 to the system studied here, we calculate azanone hydrogen complexes with HF and H₂O (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 functional 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 ESI Section 5.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.

4. Conclusion

With the development of the present manuscript, we have what we believe represents a minimally complete path for the accurate determination of intermolecular interaction energies on a NISQ-type computer. Our previous study³⁷ 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. This results in a level of SAPT that is expected to produce accurate interaction energies with or near chemical accuracy.⁴⁶ In this hybrid quantum-classical procedure, we obtain monomer wavefunctions on a 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 a 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 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 the first step in alleviating 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 simplification and a lack of density fitting of the response functions.⁵⁷ The implementation of the latter would alleviate the limitation to small basis sets, as medium size basis sets are required for accurate interaction energies with second order SAPT methods. 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^88,89 could significantly reduce the number of required measurements when the time comes to deploy this method in the presence of shot noise. As the present work used ideal statevector simulators, the effect of noisy RDMs remains unclear. However, efforts are already underway to gauge this effect. 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.

Data availability

Data for all application examples availe in the ESI.† In addition, XYZ structure and NPZ files for the MnNO⋯HF example. We also include python code to obtain and save the active space Hamiltonian via PySCF.

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. M. L. wrote the initial draft of the manuscript. All authors assisted with the production and analysis of the application results and with the writing of the manuscript.

Conflicts of interest

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

Acknowledgements

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 graphic design of Fig. 1.

References

1. Y. Cao, J. Romero and A. Aspuru-Guzik, IBM J. Res. Dev., 2018, 62, 6.
2. B. Bauer, S. Bravyi, M. Motta and G. K. L. Chan, Chem. Rev., 2020, 120, 12685.
3. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864.
4. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133.
5. M. Reiher, N. Wiebe, K. M. Svore, D. Wecker and M. Troyer, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 7555.
6. Z. Li, S. Guo, Q. Sun and G. K. Chan, Nat. Chem., 2019, 11, 1026.
7. S. R. White and R. L. Martin, J. Chem. Phys., 1999, 110, 4127.
8. G. K. L. Chan and S. Sharma, Annu. Rev. Phys. Chem., 2011, 62, 465.
9. A. Lüchow, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 388.
10. Y. Garniron, A. Scemama, E. Giner, M. Caffarel and P. F. Loos, J. Chem. Phys., 2018, 149, 064103.
11. D. S. Levine, D. Hait, N. M. Tubman, S. Lehtola, K. B. Whaley and M. Head-Gordon, J. Chem. Theory Comput., 2020, 16, 2340.
12. A. A. Holmes, N. M. Tubman and C. Umrigar, J. Chem. Theory Comput., 2016, 12, 3674.
13. S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi and C. J. Umrigar, J. Chem. Theory Comput., 2017, 13, 1595.
14. A. Aspuru-Guzik, A. D. Dutoi, P. J. Love and M. Head-Gordon, Science, 2005, 309, 1704.
15. S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin and X. Yuan, Rev. Mod. Phys., 2020, 92, 015003.
16. 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., 2019, 119, 10856.
17. J. Preskill, Quantum, 2018, 2, 79.
18. J. Cohn, M. Motta and R. M. Parrish, PRX Quantum, 2021, 2, 040352.
19. 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., 2022, 94, 015004.
20. 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, 515, 1.
21. J. R. McClean, J. Romero, R. Babbush and A. Aspuru-Guzik, New J. Phys., 2016, 18, 023023.
22. R. M. Parrish and P. L. McMahon, 2019, preprint, arXiv:1909.08925,  DOI:10.48550/arXiv.1909.08925Focus.
23. N. H. Stair, R. Huang and F. A. Evangelista, J. Chem. Theory Comput., 2020, 16, 2236.
24. W. J. Huggins, J. Lee, U. Baek, B. O'Gorman and K. B. Whaley, New J. Phys., 2020, 22, 073009.
25. K. Klymko, C. Mejuto-Zaera, S. J. Cotton, F. Wudarski, M. Urbanek, D. Hait, M. Head-Gordon, K. B. Whaley, J. Moussa and N. Wiebe, et al. , PRX Quantum, 2022, 3, 020323.
26. A. Asthana, A. Kumar, V. Abraham, H. Grimsley, Y. Zhang, L. Cincio, S. Tretiak, P. A. Dub, S. E. Economou, E. Barnes, et al., arXiv, 2022, preprint, arXiv:2206.10502,  DOI:10.48550/arXiv.2206.10502Focus.
27. W. J. Huggins, B. A. O’Gorman, N. C. Rubin, D. R. Reichman, R. Babbush and J. Lee, Nature, 2022, 603(7901), 416–420.
28. V. E. Elfving, B. W. Broer, M. Webber, J. Gavartin, M. D. Halls, K. P. Lorton, and A. Bochevarov, arXiv, 2020, preprint, arXiv:2009.12472,  DOI:10.48550/arXiv.2009.12472.
29. H. Liu, G. H. Low, D. S. Steiger, T. Häner, M. Reiher and M. Troyer, Mater. theory, 2022, 6, 11.
30. 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, Int. J. Quantum Chem., 2022, 122(22), e26975.
31. J. E. Rice, T. P. Gujarati, M. Motta, T. Y. Takeshita, E. Lee, J. A. Latone and J. M. Garcia, J. Chem. Phys., 2021, 154, 134115.
32. I. H. Kim, Y. H. Liu, S. Pallister, W. Pol, S. Roberts and E. Lee, Phys. Rev. Res., 2022, 4, 023019.
33. 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, 2022, preprint, arXiv:2203.15546,  DOI:10.48550/arXiv.2203.15546.
34. A. Heifetz, Quantum Mechanics in Drug Discovery, Springer, New York, 2020.
35. 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, 2021, preprint, arXiv:2111.12437,  DOI:10.48550/arXiv.2111.12437.
36. J. J. Goings, A. White, J. Lee, C. S. Tautermann, M. Degroote, C. Gidney, T. Shiozaki, R. Babbush and N. C. Rubin, Proc. Natl. Acad. Sci., 2022, 119(38), e2203533119.
37. 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, 13(11), 3094–3108.
38. J. Rezac and P. Hobza, Chem. Rev., 2016, 116, 5038.
39. M. Kodrycka and K. Patkowski, J. Chem. Phys., 2019, 151, 070901.
40. S. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553.
41. J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush and H. Neven, Nat. Commun., 2018, 9, 4812.
42. B. Jeziorski, M. Bulski and L. Piela, Int. J. Quantum Chem., 1976, 10, 281.
43. B. Jeziorski, R. Moszynski and K. Szalewicz, Chem. Rev., 1994, 94, 1887.
44. K. Patkowski, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2020, 10, e1452.
45. T. Korona, Mol. Phys., 2013, 111, 3705.
46. T. M. Parker, L. A. Burns, R. M. Parrish, A. G. Ryno and C. D. Sherrill, J. Chem. Phys., 2014, 140, 094106.
47. T. Korona, R. Moszynski, and B. Jeziorski, in Advances in quantum chemistry, Elsevier, United States, 1997, vol. 28, pp. 171–188.
48. M. Hapka, M. Przybytek and K. Pernal, J. Chem. Theory Comput., 2021, 17, 5538.
49. M. E. Casida, C. Jamorski, K. C. Casida and D. R. Salahub, J. Chem. Phys., 1998, 108, 4439.
50. J. R. McClean, M. E. Kimchi-Schwartz, J. Carter and W. A. De Jong, Phys. Rev. A, 2017, 95, 042308.
51. T. Korona, Phys. Chem. Chem. Phys., 2008, 10, 6509.
52. T. Korona, J. Chem. Theory Comput., 2009, 5, 2663.
53. M. Hapka, M. Przybytek and K. Pernal, J. Chem. Theory Comput., 2019, 15, 6712.
54. M. Hapka, M. Przybytek and K. Pernal, J. Chem. Theory Comput., 2021, 17, 5538.
55. T. Korona, J. Chem. Phys., 2008, 128, 224104.
56. T. Korona, Phys. Chem. Chem. Phys., 2008, 10, 6509.
57. T. Korona and B. Jeziorski, J. Chem. Phys., 2008, 128, 144107.
58. K. Chatterjee and K. Pernal, J. Chem. Phys., 2012, 137, 204109.
59. M. Hapka, M. Przybytek and K. Pernal, J. Chem. Theory Comput., 2018, 15, 1016.
60. R. M. Parrish, E. G. Hohenstein, P. L. McMahon and T. J. Martínez, Phys. Rev. Lett., 2019, 122, 230401.
61. X. Cai, W. H. Fang, H. Fan and Z. Li, Phys. Rev. Res., 2020, 2, 033324.
62. 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., 2020, 2, 043140.
63. Y. Matsuzawa and Y. Kurashige, J. Chem. Theory Comput., 2020, 16, 944.
64. L. B. Maia and J. J. Moura, Chem. Rev., 2014, 114, 97.
65. S. Moncada, Pharmacol. Rev., 1991, 43, 109.
66. T. D. Crawford, E. Kraka, J. F. Stanton and D. Cremer, J. Chem. Phys., 2001, 114, 10638.
67. R. Bruckdorfer, Mol. Aspects Med., 2005, 26, 3.
68. N. Lehnert, H. T. Dong, J. B. Harland, A. P. Hunt and C. J. White, Nat. Rev. Chem., 2018, 2, 278.
69. K. Kadish, K. M. Smith, and R. Guilard, The Porphyrin Handbook, Elsevier, San Diego, 2000, vol. 3.
70. 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.
71. A. P. Hunt and N. Lehnert, Acc. Chem. Res., 2015, 48, 2117.
72. J. Santolini, J. Inorg. Biochem., 2011, 105, 127.
73. F. A. Walker, J. Inorg. Biochem., 2005, 99, 216.
74. T. W. Hayton, P. Legzdins and W. B. Sharp, Chem. Rev., 2002, 102, 935.
75. M. Radon and K. Pierloot, J. Phys. Chem. A, 2008, 112, 11824.
76. 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., 2021, 121, 14682.
77. E. Tfouni, F. G. Doro, L. E Figueiredo, J. CM Pereira, G. Metzke and D. W Franco, Curr. Med. Chem., 2010, 17, 3643.
78. J. R. Hickok and D. D. Thomas, Curr. Pharm. Des., 2010, 16, 381.
79. R. Serafim, M. Primi, G. Trossini and E. Ferreira, Curr. Med. Chem., 2012, 19, 386.
80. T. Ampßler, G. Monsch, J. Popp, T. Riggenmann, P. Salvador, D. Schröder and P. Klüfers, Angew. Chem., Int. Ed., 2020, 59, 12381.
81. K. Boguslawski, C. R. Jacob and M. Reiher, J. Chem. Theory Comput., 2011, 7, 2740.
82. P. E. Siegbahn, M. R. Blomberg and S. L. Chen, J. Chem. Theory Comput., 2010, 6, 2040.
83. L. E. Goodrich and N. Lehnert, J. Inorg. Biochem., 2013, 118, 179.
84. M. Radon, E. Broclawik and K. Pierloot, J. Phys. Chem. B, 2010, 114, 1518.
85. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553.
86. N. Mardirossian and M. Head-Gordon, Phys. Chem. Chem. Phys., 2014, 16, 9904.
87. N. Mardirossian and M. Head-Gordon, Mol. Phys., 2017, 115, 2315.
88. M. Motta, E. Ye, J. R. McClean, Z. Li, A. J. Minnich, R. Babbush and G. K. Chan, NPJ Quantum Inf., 2021, 7, 1.
89. W. J. Huggins, J. R. McClean, N. C. Rubin, Z. Jiang, N. Wiebe, K. B. Whaley and R. Babbush, NPJ Quantum Inf., 2021, 7, 1.
90. O. Ryosuke, U. Yuya, N. Daisuke, H. Shohei and L. Crissman, CuPy: A NumPy-Compatible Library for NVIDIA GPU Calculations, Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS), 2017, http://learningsys.org/nips17/assets/papers/paper_16.pdf.

---

Footnotes

† Electronic supplementary information (ESI) available: Computational details, extended derivation of the ERPA-based SAPT terms, more background on the VQE algorithm used, additional data for the application examples, XYZ structures and NPZ files for the MnNO⋯HF example. We also include the python code to obtain and save the active space Hamiltonian via PySCF. See DOI: https://doi.org/10.1039/d2sc05896k

‡ These two authors contributed equally to this work.

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

---