Estimation of Electrostatic Interaction Energies on a Trapped-Ion Quantum Computer

We present the first hardware implementation of electrostatic interaction energies by using a trapped-ion quantum computer. As test system for our computation, we focus on the reduction of NO to N2O catalyzed by a nitric oxide reductase (NOR). The quantum computer is used to generate an approximate ground state within the NOR active space. To efficiently measure the necessary one-particle density matrices, we incorporate fermionic basis rotations into the quantum circuit without extending the circuit length, laying the groundwork for further efficient measurement routines using factorizations. Measurements in the computational basis are then used as inputs for computing the electrostatic interaction energies on a classical computer. Our experimental results strongly agree with classical noise-less simulations of the same circuits, finding electrostatic interaction energies within chemical accuracy despite hardware noise. This work shows that algorithms tailored to specific observables of interest, such as interaction energies, may require significantly fewer quantum resources than individual ground state energies would require in the straightforward supermolecular approach.


I. INTRODUCTION
The potential of quantum computers to simulate molecules is promising for many applications in industry and research [1][2][3].While quantum computers which could simulate classically intractable molecules yet have to be built, the task of identifying relevant industrial applications of these future machines is equally crucial.This goal involves not only identifying classically challenging molecules that could benefit from quantum computations [4][5][6][7], but also developing quantum algorithms which perform the computations required for specific applications in industry.
To date, most quantum algorithms have been developed to capture more accurate ground state total energies [8,9].However, accurate computations of properties beyond ground state energy are essential for many industrial applications [10][11][12].A prime example is the estimation of interaction energies between two molecules, a critical first step in computational drug design for ranking the efficacy of ligands against a target such as a protein [13].On a quantum computer, interaction energies could be obtained through three distinct ground state energy computations: one for each molecule and one for the entire system.Given that interaction energies typically constitute only a fraction of the total ground state energies, this method -known as the supermolecular approach-requires a high precision in each ground state energy calculation [14,15].Symmetry-adapted perturbation theory (SAPT) offers an alternative approach by expressing the interaction energy from a perturbative treatment of the intermolecular potential [16][17][18].Recent theoretical advancements have laid the foundation for harnessing emerging quantum computing capabilities to compute more accurate interaction energies via SAPT [19,20].Specifically, it was demonstrated that a near-term quantum algorithm, the variational quantum eigensolver (VQE) [8], could enhance the accuracy of SAPT calculations and provide an alternative to the supermolecular approach.However, this was only tested via the emulation of the quantum computer on classical computers.
In this work, we present the first experimental demonstration toward the computation of SAPT energies on a quantum computer, focusing specifically on the quantum computation of the electrostatic energy.We note that the computation of the electrostatic term by itself does not require symmetry adaption.However, we still refer to SAPT interaction energies as this work presents the first step toward the goal of implementing the full second order SAPT energies on quantum hardware [19,20].
An important part of this quantum algorithm is its embedding into a classical framework.Therefore, it is desirable to apply it to molecules whose sizes and complexities are representatives of standard problems in industrial use cases.This allows to resolve any classical bottleneck, ensuring a smooth transition to an era with improved quantum hardware.Hence, as a test system for our experiments, we focus on the active site of a nitric oxide reductase (NOR), a key player in the chemical reaction that reduces nitric oxide (NO) to nitrous oxide (N 2 O), i.e. an important step of the nitrogen cycle.Several intermediates in the catalytic cycle have electronic structures that are strongly correlated and thus challenging to simulate with classical methods.As a result, the NOR active site serves as a valuable benchmark system for quantum algorithms.Notably, NOR is a member of the cytochrome P450 superfamily, renowned for its monoxygenases that play a central role in drug metabolism.The fact that NOR belongs to the P450 superfamily renders the exploration of the catalytic cycle of NO reduction not only intrinsically valuable but also broadly significant due to its structural parallels with other P450 enzymes.
In order to make our computations possible on today's quantum computers with limited qubit counts and gate fidelities, we rely on classical quantum chemistry methods to preprocess and heavily simplify the system.Specifically, while we treated the studied molecules with more than a thousand orbitals classically, we only utilize four orbitals (mapped to eight qubits) to account for the molecule's strong correlation on the quantum computer.Despite the quantum algorithm only modeling a very small part of the system, a task that could be readily handled by any classical computer, our work serves as the first demonstration of such a computation on quantum hardware.An interesting feature emerging from this experiment is the measurement in a single basis allowed by the implementation of a fermionic basis change in the quantum circuit.The same circuit primitives could also be used to measure two-particle density matrix (2-PDM) with the help of double factorization [9,21,22], greatly reducing the number of measurement circuits compared to the "naïve" approach.This feature highlights our work also as a demonstration of experimental capabilities relevant beyond the SAPT framework.
In Sec.II we review the SAPT method and explain how to efficiently compute the first term of its perturbation series, the electrostatic contribution.Next, in Sec.III, we highlight the crucial classical steps that account for the preparation of the molecular systems and the optimization of the quantum circuits.Finally, in Sec.IV, we review the implementation of the quantum circuits on quantum hardware and discuss the experimental results.

II. QUANTUM COMPUTATIONS OF THE ELECTROSTATIC ENERGY
The interaction energy between two monomers A and B is given in the supermolecular approach by where E AB , E A , and E B denote the ground state energies of the dimer, monomer A, and monomer B, respectively.Consequently, using this approach, the computation of the interaction energy requires three separate groundstate energy calculations, one of which is performed on the larger dimer system.An alternative approach to computing interaction energies is SAPT [16,17].In SAPT, the interaction between the two molecules is treated as a perturbation, and the solution of the symmetrized Rayleigh-Schrödinger equation yields the interaction energy as a perturbation series, written in polarization terms, E pol and exchange terms, E exch , as where n is the perturbation order.The polarization and exchange terms of Eq. ( 2) are commonly relabelled to their interaction energy kind.The four main interaction types are the electrostatic energy, the exchange energy, the induction energy and the dispersion energy.The exact expression of these main terms depends on the level of truncation in the chosen SAPT method.We direct the reader to Refs.[23,24] for a more detailed introduction to SAPT.It is important to remark here that on top of providing insight to the underlying forces of the intermolecular interaction, this approach also eliminates the need for a ground state energy calculation of the dimer system and the propagation of large errors from energy calculations into typically small energy differences.
In this study, we focus on the quantum computation of a single interaction type, E elst , which describes the electrostatic energy between the two monomers.We employ the full configuration interaction (FCI) SAPT formalism of Ref. [25] wherein a complete (or VQE-type near-complete) treatment of electron correlation is assumed at each level of intermolecular perturbation.In this case, the electrostatic energy can be calculated via where γ A pp ′ is the spin-summed one-particle density matrix (1-PDM) of monomer A in a (non-orthogonal) spinrestricted spatial "atomic orbital (AO)" basis {ϕ p (⃗ r)} and is the generalized electrostatic potential matrix.Here, N A and N B are the numbers of electrons of monomer A and monomer B respectively.The first term comprises the two-body integrals, (pp ′ |qq ′ ) and the spin-summed 1-PDM of monomer B, γ B qq ′ .In the second and third terms, V X pp ′ is the nuclear potential from monomer X and S pp ′ the overlap matrix in the AO basis.In the last term, V AB is the inter-monomer nuclear repulsion term.All formal definitions can be found in Appendix A. In this work, the orbital basis is always defined as the dimer basis.
The expression of Eq. ( 3) is particularly advantageous when monomer B is well described by classical methods such as Hartree-Fock (HF), density functional theory (DFT) or other low polynomial-scaling approaches as, in this case, J B pp ′ can easily be obtained in classical preprocessing.The remaining task is to find the (crucial) interaction of the electronic density of A with the full electrostatic potential of B. To do so γ A pp ′ needs to be computed.When the ground state wavefunction of A, |Ψ A ⟩, is represented in whole, or in part, on a quantum computer, the use of an orthonormal "molecular orbital (MO)" spatial basis is commonly required.The MO basis is defined by with the orthonormality condition The 1-PDMs in the AO and MO basis are then related by In the MO basis Eq. ( 3) becomes The rewriting of E elst in the MO basis, in Eq. (8), is interesting because it makes apparent that a different set of MOs, Cpv can be chosen such that J B vv ′ is diagonal, i.e., This "electrostatic potential natural orbital basis" is simply obtained by rotating the canonical MOs as where U tv are the eigenvectors of J B tt ′ and wv the corresponding eigenvalues.Finally, Eq. ( 8) becomes with γA vv = tt ′ U tv γ A tt ′ U t ′ v .Importantly, when the wavefunction is represented on a quantum computer with a spin-restricted Jordan-Wigner formalism, the diagonal elements of the MO 1-PDM correspond to commuting, diagonal single-particle Pauli measurements where the Pauli operator Ẑv,α acts on the qubit representing MO v with spin α and | ΨA ⟩ is the ground state of A in the electrostatic potential natural orbital basis.Note that the determination of the ground state |Ψ A ⟩ can still be performed as usual with the canonical MOs and the transformation to the electrostatic potential natural orbitals can simply be implemented by decomposing the fermionic basis rotation U tv into Givens rotations [26,27] and appending them to the quantum circuit right before the measurements.
Ultimately, we see that this procedures allows to obtain E elst from a single set of measurements since all the operators which expectation values are required to reconstruct γ A vv , as defined in Eq. ( 12), commute with each other.The employed method contrasts with the "naïve" approach where the number of distinct measurement circuits to reconstruct the 1-PDM grows with quadratically with the number of orbitals N .
In comparison, the supermolecular approach, as per Eq. ( 1), requires the calculation of the ground state energies of the both monomer and dimer systems, which confers a clear advantage to the perturbative alternative.Naively, this method requires up to O(N 4 ) measurement circuits.However, an interesting parallel with the present work is that the same change of fermionic basis approach can be used to measure the 2-PDM with the help of double factorization [9,21,22], reducing the number of measurement circuits.Furthermore randomizing over these basis rotations would allow for a classical shadow approach only scaling with the number of occupied fermionic modes [28].

III. CLASSICAL PROCESSING
As a test system, we exploit the reduction of NO to N 2 O catalysed by NOR.Details on this reaction, such as the proposed catalytic cycle are discussed in Appendix B. Ref. [29] discusses the importance of intermolecular interactions in the reaction mechanism.For instance, water molecules surrounding the NOR active site form a hydrogen bonding network and act as a proton delivery pathway for the formation of one key intermediate.One intermediate step involves a hydride transfer from a surrounding NADH molecule to the NO-bound heme (see Appendix B).The NO-bound heme is expected to exhibit strong correlations and hence can be a good target for quantum computations.The strong correlation arises as π/π * orbitals of both NO and porphyrin moieties are known to mix wit the 3-d orbitals FIG. 1. Visualization of the quantum computation of the electrostatic interaction energy.(left) The output of the classical Hartree-Fock calculations is fed into the quantum computer, where a previously generated approximate ground state of monomer A is rotated in the basis given in Eq. (11).Measuring the diagonal elements of the 1-RDM allows for an efficient calculation of the electrostatic energy.(right) Visualization of the quantum circuit, including preparing the approximate ground state and the subsequent fermionic basis rotation, see also Fig. 3.
of the Fe center [5,30,31].In this context, we inquire whether electrostatic interactions favour the stabilisation of the NADH molecule in the heme pocket.In practice, this involves calculating the difference in electrostatic energies between both intermediates, denoted as The preparation of the model systems for the two intermediates of interest, shown in Fig. 2, is discussed in great detail in Appendix B. To calculate the electrostatic interactions, we break down both intermediates into two separate monomers.For intermediate A, the first monomer consists of the NO-bound heme, while the second monomer encompasses the amino acids and water molecules.In the case of intermediate B, the first monomer includes both the NO-bound heme and the NADH molecule, while the amino acids and the remaining water molecule form the second monomer.In both cases the ground state wavefunction of the first monomer is partially treated with a quantum computer.To match with the hardware constrains, we select a small four electrons in four MOs (4e, 4o) active space.Capturing the important electronic correlations in our systems would obviously require a much larger active space [5] and the choice we make is strictly driven by the experimental setup.However, the orbitals selected by looking at the natural orbital occupation numbers (see Appendix B and Fig. 2) include character from the iron 3d, NO (anti)bonding π, porphyrin (anti)bonding π orbitals which govern the strong corrleation in this system [5,30,31].
To find the ground state energy in the active space, the VQE algorithm [8] is used.The fermionic Hamiltonian in the (4e, 4o) active space is mapped through the Jordan-Wigner mapping onto a 8-qubit system, in α-then-β ordering, meaning the first (last) 4 qubits are assigned to the α/spin-up orbitals (β/spin-down orbitals) with a qubit in state |1⟩/|0⟩ representing an occupied/unoccupied orbital.We adapt the VQE ansatz from Ref. [27]  to generate an approximate ground state.This ansatz is composed out of an alternating application of Givens rotations, G(θ) (i.e.fermionic basis rotations represented after Jordan-Wigner transformation) and fermionic double excitations, P X (θ), where we restrict our ansatz to a single layer of each, see Fig. 3 and Appendix D. The set FIG. 3. The quantum circuit used in this work.The two-qubit gates labeled with G denote Given rotations, which are equivalent to local fermionic basis rotations after a Jordan-Wigner mapping.The four-qubit PX gate represents a PairExchange gate [27].We refer to Appendix D for more details.The Givens rotations of the VQE circuit were merged with the basis rotation circuit (hatched area) to reduce the circuit depth.
of parameters yielding an approximation of the ground state is found by optimizing the expectation value of the second quantized Hamiltonian with the L-BFGS-B algorithm using a classical simulator implemented in Qiskit [32].Independent from the number of layers, the last action of the ansatz corresponds to a Givens rotation network.As the electrostatics require only the diagonal part of the 1-PDM γvv in the electrostatic potential natural orbital basis, we implement the fermionic basis rotation U tv that diagonalizes J B tt ′ by another network of Givens rotation gates identical in structure to this last part of the ansatz.We leverage the property that the multiplication of two Givens rotation networks is itself another Givens rotation network, see Appendix D for details.This allows us to merge the two Givens rotation networks at the end of our circuit and reduce the circuit depth, see Fig. 3.We finally measure the qubits in their computational basis to obtain the expectation value of the electrostatics operator in the active space and combine it with the contribution from the core as explained in Appendix C.

IV. HARDWARE IMPLEMENTATION & RESULTS
As quantum hardware, we use the quantum processor aqt marmot hosted by AQT [33].The aqt marmot system is based on trapped 40 Ca + ions and supports a universal set of gates.The native gate set comprises singlequbit gates with arbitrary rotation angles and axis.The entangling operation is a two-qubit Mølmer-Sørensen (MS) gate [34].The MS gate allows for entangling operations with arbitrary rotation angles that can be implemented between any qubit pair (see Appendix D).In this work, we utilize an 8-qubit register featuring allto-all connectivity.We choose R X (θ) = exp[−iθ/2 Xi ], R Z (θ) = exp[−iθ/2 Ẑi ] and R XX (θ) = exp[−iθ Xi Xj ] with Xi and Ẑi being the Pauli-X and Pauli-Z matrices acting on qubit i respectively, as the basis gate set since it closely resembles the native gate set of the aqt marmot system.The gate errors can be approximated by depolarizing noise acting on the addressed qubits.The error rates for the single-qubit gates (R X error) are approximately 3 • 10 −4 on average, whereas the error rate for the two-qubit gate (R XX error) is around 1.5 • 10 −2 on average.The typical gate times are 15 µs for single-qubit gates and 200 µs for two-qubit gates.The T 1 and T 2 times are 1.14 ± 0.06 seconds and 0.452 ± 0.068 seconds respectively, resulting in a coherence/gate time ratio of about 10 3 [35].
We use the quantum circuit described in Section III to prepare an approximate ground state on the quantum computer.We compile the quantum circuits' gates into the hardware's native gate set to implement the circuit on the trapped-ion quantum computer.Subsequently, we use Qiskit's transpiler [32] to further optimize the circuit and minimize the total number of gates.After optimization, our circuit consists of 279 single-qubit gates and 63 two-qubit gates.The transpiled circuit for intermediate A is displayed in Appendix D.
In Figure 4(a) we present the output statistics obtained from N meas.= 4 × 10 4 measurements for intermediate A and B. We compare the experimental results to results from an error-free VQE simulation on a classical computer.We also put these distributions in contrast with the uniform one (random sampling).Despite a nonnegligible level of noise, we observe strong agreement in the dominant computational basis states.To obtain the electrostatic interaction energies for both intermediates, we first exclude all computational states with incorrect alpha and beta particle numbers from the measured data, technique commonly referred to as postselection.In both experiments, we discard roughly 50% of the samples.We then construct the diagonal parts of the rotated 1-PDM using Eq. ( 12) and compute the electrostatics according to Eq. (11).We summarize our findings in Fig. 4(b).We find that the electrostatic energies of both systems are within chemical accuracy (1 kcal mol −1 ) to the error-free VQE simulations, with deviations of ∆E = −0.204(45)kcal mol −1 and ∆E = −0.354(30)kcal mol −1 for intermediate A and intermediate B, respectively.To estimate the statistical errors, we calculate the electrostatics energy for each measurement with correct particle number separately.From the list of electrostatics, we determine the standard error of the mean, given by σ E = σ/N cor., where σ is the standard deviation of the list and N cor. the number of measurements with correct particle number.In Fig. 4(c), we also report the convergence of the electrostatics when increasing the number of measurements to construct the required 1-PDM.Here, for clarity, we report estimates obtained from using more than 1000 measurements.In Appendix E we show the same convergence results using 2 to 1000 measurements.In this case, the calculated electrostatic energies lie in a much broader range (from −26.98 kcal mol −1 to 6.80 kcal mol −1 ).This emphasizes that the spectrum of our estimate is non-trivial but that its expectation value converges rapidly when increasing the number of measurements.This also indicates that the remaining ∆E is a consequence of device noise biasing our estimate.In Appendix E, we also explore results from random sampling, which yields fairly accurate electrostatic energies.By comparing the output statistics between experimental results, noise-free simulations, and random sampling, we argue that the observed phenomenon is an attribute of the system under investigation and does not undermine the experimental outcome.
From the experimental results, we calculate a difference in electrostatic energies of ∆E exp.elst = −31.32(5)kcal mol −1 between intermediate A and B, allowing us to answer the question raised in Sec.B and conclude that the electrostatics favors the stabilization of NADH in the heme pocket.Although, for this model system, this could have been qualitatively described using Hartree-Fock theory (∆E RHF elst = −28.57kcal mol −1 ), our results quantitatively align with FCI calculations (∆E FCI elst = −31.14kcal mol −1 ).
Multiple efforts to reduce the noise bias ∆E were explored: To investigate and mitigate potential biases stemming from keeping the ion-qubit assignment fixed during the quantum computation, we performed additional ex-periments, where we introduced a random alteration of the qubit-ion assignment every 100 experimental runs.Additionally, to mitigate the impact of qubit decay from |1⟩ → |0⟩, we treated the qubits' zero state |0⟩ to represent an occupied orbital, while |1⟩ represent an unoccupied orbital.This can be easily implemented by adding a full layer of Pauli-X gates in the beginning of the quantum circuit while concurrently substituting all circuit parameters with their negated values.As demonstrated in Appendix E, incorporating these error mitigation techniques does not significantly improve the accuracy of the electrostatic energy.We also classically simulated the effect of Zero-Noise-Extrapolation (ZNE), an error mitigation method well-suited for depolarizing errors as they predominate in the used quantum processor (see Appendix G).While we found an improvement of the electrostatic energy, the error of the estimate increased significantly.To reduce the error, a substantial increase in the number of experimental measurements would have been necessary.
In order to draw a comparison between our approach and the supermolecular approach given in Eq. ( 1), we would require to estimate the ground state energies of both monomers and the larger dimer system.However, since we only measure the diagonal elements of the rotated 1-PDM of monomer A, as per Eq. ( 12), we cannot use these measurements to estimate the ground state energy of monomer A, as this would necessitate the full 1-PDM and the 2-PDM.Instead, we use these measurements to compute the expectation value of the diagonal elements of the one-body part of the Hamiltonian, with further details provided in App.F. An deviation of 1.16 kcal mol −1 was observed when comparing the experimental results with the noise-less simulation.Considering the error in the diagonal contribution to the one body energy, it is unlikely that the accuracy of the ground state energies of the monomer systems would be sufficient to estimate interaction energies within chemical accuracy, let alone the ground state energy of the larger dimer system.

V. CONCLUSION
In our study, we performed the first direct quantum computations of interaction energies; specifically we computed the electrostatic energies of two important intermediates in the catalytic conversion of NO to N 2 O by P450nor.Using a VQE ansatz, we generated an approximate ground state for one of both monomers on a trapped-ion quantum computer, while the other was treated with classical HF theory.A fermionic basis rotation at the end of the quantum computation allowed the efficient measurement of the electrostatics energy with a single quantum circuit and without extending the circuit length.We found that the quantum-computed electrostatics strongly agree with our simulation results, within a chemical accuracy of 1 kcal mol −1 .As measuring under these fermionic basis rotations is also used when e.g.measuring ground state energies efficiently with double factorization [9,21,22], this experiment implies that these techniques are starting to be in reach of current hardware.
In order to execute the demonstrated quantum experiment, we had to considerably simplify the computational problem to conform to limited number of qubits, coherence times and gate fidelities available on the quantum computer.Furthermore, we employed classical simulations of the quantum circuit to optimize the variational parameters.This allowed us to leave aside the task of finding the ground state in a noisy environment and focus our study solely on the computation of the electrostatic energy.While this work does not offer new insights into the systems' chemistry, it does provide the first experimental demonstration of the potential usefulness of quantum algorithms which are tailored to the computation of quantities beyond ground state energies.In particular, we showed how the measurement count can be reduced in practice, resulting, for the small models considered, in a shot noise resilient estimate.
Scaling our approach to industry-relevant problems requires further research.First, the electrostatic energy is only a single term of SAPT, which by itself does not require symmetry adaption.Obtaining other terms such as the exchange, induction and dispersion is vital to accurately describe the interaction between two molecules.To this end, efficient measurement strategies similar to this work must also be developed.Another crucial aspect for future research resides in finding efficient measurement schemes when both monomers exhibit strong correlations and require accurate descriptions on quantum computers.
While tailored quantum computations, such as the one presented, allow to push the boundaries of existing quantum hardware, it is unlikely that a quantum advantage will become possible without exploiting any form of error mitigation or limited error correction.Faulttolerant quantum computations of SAPT energies [36] could potentially offer another approach to the quantum computation of interaction energies.and the NO ligand.In the model of intermediate B, additional steps are required to model the active site.A water (WAT565) from the X-ray structure, hydrogen-bonded to Thr243 is included.The carboxylic acid group is modified to a carboxamide group (as in NADH) so that the carbonyl group still forms a hydrogen bond to the backbone NH from Ala239.The pyridine ring is reduced by adding a hydrogen atom (to have the reduced form: NADH), and the ligand is capped after the first ribose-carbon.The rest of the NADH is assumed not to contribute to the hydride-transfer reaction.The NO-ligand is manually put to the heme to finalize the model system with all relevant and correct species.Then, a force field optimization of the hydrogen placements in MOE (AMBER10:EHT) is done, followed by a tethered minimization (tether on all heavy atoms).An overlay of the two resulting model structures is shown in Fig. 6.
Finally, both structures are optimized with DFT.The position of the heme is held fixed and the position of the C α atoms of the amino acids is also constrained as a surrogate for the protein backbone environment that keeps amino acid side chains and co-factors in place.We use the GPU accelerated DFT implementation of Promethium [39] on a single A100 GPU, allowing us to perform the full geometry optimization at B3LYP-D3/def2-SV(P) level of theory in 3h58 for intermediate A (1042 basis functions) and 6h36 for intermediate B (1130 basis functions).The resulting geometries are shown in Fig. 7.
Because we are limited quantum resource available on current hardware, we devise an active space of four MOs.The size of this active space is solely chosen to fit the hardware constraints and does not reflect the important chemistry in our systems.Therefore, we choose to employ an automated active space selection approach rather than picking the orbitals fully based on their symmetries.To do so, we first run a semistochastic heat bath configuration interaction (SHCI) (as implemented in Dice [40,41]) for the two systems in the singlet state with ϵ = 0.01, the def2-SVP basis set, in a 20 electrons in 20 orbitals (taken around the HOMO-LUMO gap).This starting active space is sizeable enough to include the iron 3d, iron 4s, iron-nitrogen anti-bonding and axial ligand (anti-)bonding orbitals [5].From the resulting 1-PDM, we obtain the natural orbitals and order them according to their occupation numbers.We show the four orbitals around the gap in Fig. 8(left) along with their occupation numbers.These orbitals are a mixture of the iron 3d, the iron-nitrosil (anti)bonding π, the sulfur p as well as heme π orbitals which are all important to capture the essential mixing.Hence, as a final step we run a complete active space self-consistent field (CASSCF) calculation with PySCF [42] using the SHCI natural orbital as a starting guess and a 4 electrons in 4 orbitals active space (4e, 4o).The resulting orbitals define our final (4e, 4o) active space.They are displayed in Fig. 8  In the following, we describe how to recover the electrostatic energy in the case where the quantum computer is used to generate an approximate ground state in an active space.
Since the 1-PDM is block diagonal, we can rewrite Eq. 8 as where in the first (second) term the indices run over the MOs in the core (active) space.In the core space, by definition γ A tt ′ is diagonal (with 2s on its diagonal).Hence, we only rotate the orbital basis in the active space and, as per Eq. ( 11), we find that two Givens rotation network can be easily merged.As every Givens rotation network represents a fermionic basis rotation U tv , two Givens rotation networks can be easily merged by multiplying the underlying fermionic basis rotations.From the found new basis rotation, we generate a new Givens rotation network.The circuit architecture together with the merging of the Givens rotation networks is more clearly depicted on Fig. 9.
As explained in the main text, we choose the gates R X (θ), R Z (θ), R XX (θ) as target gate set for the Qiskit transpiler, where the R XX (θ) entangling operation is a two-qubit Mølmer-Sørensen (MS) gate that can also be described as, with c = cos (θ/2) and s=sin (θ/2).In Figure 10 and Figure 11 we show the full transpiled quantum circuit used to compute the electrostatic energy in intermediate A. The circuit comprises a total number of gates of (R Z , 168), (R X , 111), (R XX , 63).The circuit for intermediate B is similar to the circuit for intermediate A with a total number of gates of (R Z , 168), (R X , 113), (R XX , 63).
Note that for both systems, the optimized angles of the last two P X (θ) gates are small (≈ 10 −3 ).Their contribution to the energy is therefore negligible here.However, the standard Qiskit transpiler does not do circuit approximation, which is why we implement every circuit unitary as specified.In principle, the circuit size could be further minimized using approximation methods, where the level of approximation would depend on the noise of the specific hardware.While it seems trivial for a 8 qubit circuit with few gates, this emphasizes the need for robust transpilation codes or for guided circuit construction techniques [43] when scaling up to relevant size systems.Here, in the pursuit of reporting the performance of state of the art quantum hardware and software kit, we keep the resulting transpiled circuit while bearing in mind that the noise introduced by applying the previously mentioned P X (θ) gates will probably be greater than their contribution to the energy.

Appendix E: Additional experimental results
The electrostatic energy of our model system, constrained to a small active space to align with experimental constraints, exhibits minimal variation in response to state errors.To validate the reliability of our results as a consequence of successful quantum computation rather than mere chance, we present additional experimental data below.In Table I, we first compare the electrostatic energies obtained in classical simulation of an exact superposition state (random sampling) and of the VQE state.The difference between the two intermediate, ∆E elst = E elst (B) − E elst (A), is -31.751kcal mol −1 for random sampling, -31.314 kcal mol −1 for the exact VQE.These values are to be compared to the -31.140 kcal mol −1 predicted by FCI.While the VQE results bring us closer to FCI, the random sampling outcomes also fall within chemical accuracy, making it challenging to definitively conclude the true success of the experiment from looking exclusively at the electrostatic energy.For this reason, in Table II, we illustrate the correspondence between the statistics derived from the noise-less simulation of the quantum circuit and both the experimental results and the exact superposition distribution associated with random sampling.Here, the consistency of the results becomes more evident, with the overlaps of the experimental distributions significantly surpassing those obtained from random sampling.
For each intermediate, we also present 4 different experiments: two solely derived from the VQE quantum circuit (plain), and the remaining two incorporating error mitigation techniques (mitigated).In particular, we aim to mitigate potential biases coming from the difference in qubit performance and from qubit relaxation (|1⟩ → |0⟩).Therefore, in these additional experiments, we introduced a random alteration of the qubit-ion assignment every 100 experimental runs along with an inversion of the qubits' |0⟩ and |1⟩ states to represent an occupied orbital and an unoccupied orbital respectively.This last bit can be easily implemented by adding a full layer of Pauli-X gates in the beginning of the quantum circuit while concurrently substituting all circuit parameters with their negated values.The maximum absolute deviation in the resulting four electrostatic energies from the expected exact VQE results is very low, 0.293 kcal mol −1 for intermediate A and 0.341 kcal mol −1 for intermediate B. This demonstrates the concurrence of our results.Nevertheless, it is noteworthy that no discernible improvement is observed with the implementation of error mitigation protocols.
Finally, in Fig. 12, we report the convergence of the electrostatics when increasing the number of measurements to construct the required 1-PDM using from 2 to 1000 measurements.This is similar to Fig. 4(c) of the main text but with a focus on the first thousand measurements.As clearly visible, at low number of measurements, the spread of values is much higher.We note that the lowest/highest possible value of electrostatics energy within the space of computational states with correct particle sector for intermediate A is −26.98 kcal mol −1 and 6.80 kcal mol −1 respectively.This indicates that the spectrum of E elst in the computational basis states is wide and that, although the average of all of them (cf.Random) takes us close to the right value, in general, a state with the wrong structure would lead to an inaccurate E elst . Intermediate

FIG. 2 .
FIG. 2. (left) 2D representations of intermediate A and B. They include the NO-bound heme, three surrounding amino acids, water molecules and NADH (if present).In both cases, monomer A and B are shown in black and blue, respectively.(right) The four orbitals in the active space of monomer A for each intermediate.

FIG. 4 .
FIG. 4. (a) The output statistics of the VQE circuit, including the final basis rotation for intermediate A (top) and intermediate B (bottom).The experimental outcome (blue) is compared against exact classical simulation results (orange) and random sampling (black line).For each experiment, Nmeas.= 4 × 10 4 measurements were used.(b) A comparison between the electrostatic energies obtained from various classical and quantum methods.(c) The convergence of the electrostatics of intermediate A as a function of the number of measurements used to construct the 1-PDM in Eq. (12).

FIG. 5 .
FIG.5.Reaction scheme of the reduction of NO to N2O catalysed by P450nor.In our study we focus on NO-bound intermediates A and B and the question whether the electrostatic interaction favour the stabilisation of the NADH molecule in the heme pocket.

FIG. 6 .
FIG. 6. Overlay of the model complexes before DFT optimization.Intermediate A and B are shown in orange and green, respectively.Dashed lines denote hydrogen bonds and the water labels correspond to pdb entry 1CL6.
(right).In both intermediates, two orbitals show strong deviation from integer values indicating the multi-reference character of the systems and exhibit the expected metal ligand backbonding of Fe and NO and porphyrin π and π * orbitals.Note that the inclusion of the first and last orbitals in the active space (the one with occupation numbers close to 2 and 0, respectively) lowers the total energy by 18 mHa for intermediate A and 10 mHa for intermediate B, justifying their importance in the ground state.

FIG. 8 .
FIG. 8. (left) The four SHCI natural orbitals around the gap and their occupation numbers, (right) The four CASSCF molecular orbitals around the gap and their occupation numbers

FIG. 9 .
FIG. 9. (top)The quantum circuit used in this work.The two-qubit gates labeled with G denote Given rotations, which are equivalent to local fermionic basis rotations after a Jordan-Wigner mapping.The four-qubit PX gate represents a PairExchange gate, see below.The green circuits represent the VQE ansatz, while the blue circuit represents the basis rotation to efficiently measure electrostatic interaction.On the right hand side, the Givens rotations network of the VQE circuit was merged with the basis rotation circuit (hatched area).(bottom) Gate definitions of the Givens rotation gate G(θ) and the PairExchange Gate PX (θ) as introduced in[27].Matrix entries are c := cos(θ/2), s := sin(θ/2).

FIG. 10 .
FIG. 10.Part 1 of transpiled quantum circuit for intermediate A. Rx, Rz and Rxx gates are shown in green, blue and orange, respectively.

FIG. 11 .
FIG. 11.Part 2 of transpiled quantum circuit for intermediate A. Rx, Rz and Rxx gates are shown in green, blue and orange, respectively.

TABLE I .
This table provides an overview over the experiments ran for this study.The labeling plain corresponds to the basic setup described in the main text of this work.The experiments with the label Relabeling exploited a change of the qubit-ion assignement after each 100 shots.Moreover in these experiments, the labeling of 1 and 0 was exchanged.All energies are given in kcal mol −1 .