Molecular Quantum Circuit Design: A Graph-Based Approach

Science is rich in abstract concepts that capture complex processes in astonishingly simple ways. A prominent example is the reduction of molecules to simple graphs. This work introduces a design principle for parametrized quantum circuits based on chemical graphs, providing a way forward in three major obstacles in quantum circuit design for molecular systems: Operator ordering, parameter initialization and initial state preparation. It allows physical interpretation of each individual component and provides an heuristic to qualitatively estimate the difficulty of preparing ground states for individual instances of molecules.

: Illustration of the basic concepts developed in this article connecting chemical graphs with orbitals and quantum circuits. The circuit can be employed to the π-system of C 4 H 4 (sketched in the Figure) as well as the linear H 4 and BeH 2 (see the explicit examples in the main text). simulation methods as in Ref. [12] where a problem-aware but hardware-efficient circuit design was developed. Pioneering physical motivated circuit designs in electronic structure [3,[13][14][15] relied on trotterized time evolutions leading to comparable deep circuits and the need for heuristics on the optimal ordering [16][17][18] of the primitive operations. Adaptive approaches, most prominent in the form of qubit coupled-cluster [19,20], adapt-vqe [21][22][23] and related selective approaches [24] tackle the ordering problem through a greedy approach. Together with [12,25,26], these approaches often provide a middle-ground between hardware-efficient and physically-motivated design principles and are promising candidates for a path towards fully automatized and reliable circuit construction. A rule of thumb is, that the success of the (adaptive) process is more likely when it starts from a good initial state. [12,[27][28][29][30] In a similar fashion, algorithms like the quantum phase estimation [31] or quantum imaginary time evolution [32,33] will profit from improved initial states through increased success probabilities. In this work we will focus on the construction of shallow quantum circuits with reasonably local qubit connectivity using the concept of chemical graphs. The most important graph is directly represented by a classically tractable circuit, and for each additional chemical graph the quantum circuit is extended with an U R U C U † R motif illustrated in Fig. 1. Here U R rotates the molecular orbitals such that they represent the given chemical graph and U C is built from (fermionic) correlators. The developed design principles provide a way forward in all three major obstacles in quantum circuit design for molecular systems: • Operator ordering: A given chemical graph G gives rise to a specific ordering of fermionic unitary single and double excitations grouped in U G = U R U C U † R motifs. The ordering of the individual structural motifs is given by the relative importance of the graphs which can be estimated through physically-motivated rules.
• Parameter Initialization: Angles of double excitations are initialized to zero in order to guarantee an overall energy improvement. Single excitations are initialized based on suitable initial guesses for orbitals representing the chemical graph. This ensures improved convergence in the optimization of the angles.
• Initial state: The correspondence between a single chemical graph and the separable pair approximation [12] is exploited to allow for a classically simulable initial state.  (12) for the H 2 /6-31G(2,8) system using standard Hartree-Fock (HF) or optimized orbitals (opt) with FCI/6-31G (2,8) as reference. Orbital optimization is necessary to reach the exact ground state energy. Center: Performance of the related U 12 2 circuit for the LiH/STO-3G (2,12) with frozen core orbital and FCI/STO-3G (4,14) as reference, showing that orbital optimization is necessary and the frozen-core error is negligible. Right: U (4,8) SPA circuit (the same as in Eq. (18)) for LiH/MRA-PNO (4,8) with an adaptive basis that ensures significant core correlation and reference energy FCI/MRA-PNO (4,8) showing that the 4 electrons of LiH can be can be separated in a core and valence pair. More details on the LiH results are provided in the appendix. The millihartree accuracy threshold is marked with blue.

Molecular Quantum Circuit Design
In the following we will describe how molecular quantum circuits can be constructed through chemical model systems that we will organize in three layers of abstraction: the chemical graph that represents the molecule of interest through a set of vertices (nuclei) and edges (bonds), the orbitals (one-body wavefunctions in spatial space that are mapped to qubits) which can be rotated to resemble such a graph, and the quantum circuits build from those rotations combined with two-electron correlators.
We will start with two electron systems and demonstrate that an orbital-optimized form of the classically tractable separable pair approximation of Ref. [12] solves those systems exactly. The corresponding circuits can be transferred to effective two-electron systems like the LiH molecule -also a prominent proof-of-principle example in the literature. The ground states of all those systems can be approximated almost exactly (with energetic errors below the millihartree threshold) with the constructed circuits. We will use this furthermore to introduce and illustrate important concepts like orbital rotations, separable pair approximations and the frozen core approximation in an illustrative way. In the sections thereafter we will then discuss how these concepts can be extended beyond (effective) two-electron systems.

Exact and Near-Exact Circuits: (Effective) Two Electron Systems
The hydrogen molecule was used prominently in the pioneering work on variational quantum algorithms [34] and has since then served as a toy model for a myriad of approaches in variational quantum computing. Here we will use it to introduce important concepts for the following parts of this article and to clarify why it is not an ideal benchmark system for novel algorithmic approaches to quantum chemistry as its eigenstates can be prepared efficiently by a low-depth and local circuit with low variable count that belongs to a class of classically tractable quantum circuits. We will begin with the smallest representation of H 2 that still admits correlation beyond a mean-field treatment. Here the electrons are allowed to move within two spatial orbitals, usually represented by linear combinations of atomic orbitals placed on the center of the two hydrogen atoms (see A and B for more details) . We will represent this graphically as where gray dots represent the atomic centers and the dark spheres the spherical symmetric atomic orbitals. A prominent basis set that provides those atomic orbitals is the STO-3G set (S later T ype Orbitals made from 3 contracted Gaussian functions). For more clarity we will denote the number of electrons and spin-orbitals (qubits) in parenthesis, so in this case we have H 2 /STO-3G (2,4). The standard approach to represent molecules graphically is given by molecular graphs (or Lewis formulas) with atomic nuclei as vertices and chemical bonds 1 as edges. For H 2 this looks like where the edges are interpreted as electron pairs that connect two atoms. We can prepare a 4-qubit wavefunction that represents this single bond situation with the circuit given as (3) The free parameter θ enters via the parametrized e −i θ 2 σy single qubit rotation gate in (3), with the four qubits representing four spin orbitals ϕ i ↑ , ϕ i ↓ that are formed by linear combinations of the atomic spin orbitals with ψ k↑ ∈ {ψ R↑ , ψ L↑ } (analog for the spin-down orbitals). In most applications, including this work, the spin-up and spin-down orbitals are assembled through the same coefficients For our 4-qubit (2 spatial orbitals) H 2 system, we can write Eq. (4) explicitly as matrix equation (omitting the spin label) In order to describe the resulting orbitals ϕ k , we will resort to the graphical depiction of Eq. (1) and represent the coefficients c ik as small circles placed on the corresponding atomic centers. The size of the circles will correspond to the relative size between the coefficients and the color will indicate sign changes, with red being used for negative coefficients. Lets consider molecular orbitals for H 2 where all coefficients are equal in size forming a bonding and anti-bonding molecular orbital from the atomic basis functions Combined with the right choice of orbitals -in this case, the canonical orbitals from Eq. (6) -the wavefunction prepared by circuit (3) describes the 2 orbital (4 qubit) H 2 system exactly for all possible bond distances.
The orbital rotations in Eq. (4) can be represented as a unitary operation in the second quantized formulation c. a so-called fermionic basis rotation (see appendix B for more details). Translated to a quantum circuit this orbital rotating unitary looks like where the two gates represent single-electronic excitations within R and L orbitals (qubits 0 and 2 for spin-up and qubits 1 and 3 for spin-down electrons respectively), each representable with two three-qubit Pauli rotations (XZY and Y ZX -see appendix B for an explicit decomposition). Alternatively the orbitals can be optimized in a more direct way by using a first order expansion of the transformed Hamiltonian H ′ = U † R HU R and optimizing the angles φ in the resulting functional. After the optimized angles are found, a new molecular Hamiltonian can be constructed from the optimized orbitals. In the case of H 2 /STO-3G(2,4), the optimized orbitals are the symmetric and anti-symmetric combination of the atomic orbitals from Eq. (6) meaning, adding the circuit (7) to circuit (3) and optimizing the angles has the same effect as using only circuit (3) with optimized orbitals. In the following part of this work, we will frequently resort to this graphical representation of the action of from fermionic rotations U R . In this case, the optimal orbitals with respect to the wavefunction prepared by U 4 2 are identical to the canonical Hartree-Fock orbitals and their optimal linear combination is the same for all bond distances. This is however not always the case and is due to the high degree of symmetry in the H 2 /STO-3G(2,4) system.
On first glance, it might seem counter-intuitive that the orbitals in (8) are optimal even for large bond distances, where the molecule is essentially dissociated. In the dissociated limit, the optimal angle in circuit (3) is π/2 making the generated qubit wavefunction with the corresponding fermionic expression where the first line denotes the wavefunction with the optimal molecular orbitals (8) (indexed with 0 and 1) and the second line uses atomic orbitals L and R through decomposition (4) with optimal coefficients (8). The representation with atomic orbitals clearly shows, that the electrons in the ground state of the dissociated molecule are spatially separated with one being located at the left and the other at the right atomic orbital. In other words, the two atomic orbitals are never doubly occupied. On the other hand, the representation with molecular orbitals has only double occupied orbitals. This illustrates how two spatially separated electrons can be represented by using only doubly occupied orbitals in the wavefunction, given that the right choice of orbitals is used.
In a larger basis, e.g. H 2 /6-31G (2,8), the circuit in Eq. (12) is extended to In general, the N qubit form of this circuit prepares with 2 k+1 + 2 k denoting the computational basis states with two 1s at qubits 2k and 2k + 1 in binary (e.g. for k = 2 we have |48⟩ ≡ |. . . 00110000⟩). As before, this wavefunction is restricted to doubly occupied orbitals and can describe the molecule exactly given that the orbitals are optimized accordingly. In this case however, the optimized orbitals are not identical to the Hartree-Fock orbitals witnessed by the energy difference in Fig. 2.
Other than H 2 the Helium atom is another obvious two-electron system, here the ground states of the 4-qubit He/6-31G(2,4) or ten qubit He/cc-pVDZ(2,10) can be represented exactly with the same type of circuits as used for H 2 . Meaning the circuit (3) (adapted to the corresponding qubit number) preparing wavefunction (13). We can essentially view an Helium atom as H 2 with zero bond distance (omitting of course the constant nuclearnuclear potential).
The beryllium atom and the lithium hydride (LiH) molecule are formally four electron systems but they behave like effective two electron systems where one electron pair remains close to the beryllium, respectively lithium, atom. These so called core electrons do usually show only negligible correlation within most basis sets -the frozen core error within the STO-3G basis is for example below the millihartree threshold for all bond distances (see Fig. 2 and Fig. 4 in the appendix). There are specialized basis sets (e.g. the cc-pCVXZ family where the capital C denotes core-polarization) that allow to represent these correlation effects. In the case of LiH there is however still no noteworthy correlation between the pairs, so a separable pair approximation [12] U is sufficient for an accurate description (see Fig. 2 and additional details in the appendix C). In a similar manner we can describe other single bond effects, like the dissociation of the C-C bond in C 2 H 6 -see Ref. [12] where this is explicitly discussed.

Chemical Graphs and Separated Electron Pairs
In the previous section it was demonstrated how (near-) exact wavefunctions of (effective) two-electron systems can be prepared with shallow parametrized circuits combined with an optimized orbital representation. In (2) the chemical graph (or Lewis formula) of the Hydrogen molecule -where the vertices represent the two hydrogen atoms and the edge represents a chemical bond -was shortly introduced. Graph based representations were applied successfully in quantum optimization approaches in the past [35][36][37] with the graph representations of optical systems specifically developed for this purpose. For chemical systems a graph representation comes more naturally and can be constructed by determining and connecting the valence electrons N A e for each atom A. In the following a simple scheme that achieves this for neutral molecules with even number of electrons built from atoms of the the first two main-row elements of the periodic table (atomic number n A < 18) is briefly described without going into too much detail. This is sufficient for the techniques described in this work but also for the majority of current day research in variational quantum algorithms.
Chemical graphs intrinsically use the frozen-core approximation described in the last section and connect only the valence electrons (all electrons that are not core electrons). Using the periodic system of elements, the number of core electrons for atom A is N A core = n B with B being the next noble gas with atomic number n B ≤ n A . For the first three rows we have for example with the noble gases Helium (n He = 2) and Neon (n N e = 10). The number of valence electrons is then Now each atom in the molecule is assigned N A valence connectors and all that is left is to connect them to obtain the edges of the graph. Note, that the connectors can also form loops (so called lone pairs) which is often a preferred choice once N A valence > 4. 2 Usually the most important chemical graph is the one where all atoms are maximally connected to their closest neighbors (see the linear H 4 example introduced further below). Note however, that the most important graph is not always unique (a famous example being conjugated ring systems [38] like H 4 or Benzene).
Once a chemical graph is determined we can construct a circuit from it by treating all edges as separable electron pairs whose wavefuncions can be prepared by the same circuits as in (12) and the optimal orbitals are determined as the linear combination of basis orbitals that yield the lowest energy for this circuit. This is the so called separable pair approximation (SPA) [12] that was originally introduced in combination with system-adapted orbitals determined through a surrogate model [39,40] and is related to classical (generalized) valence-bond approaches (see for example [41][42][43] or [44] for a modern approach in the context of tensor networks). In Heuristic 1 an SPA is constructed from a chemical graph of a molecule (note that the first cr y can be implemented as just r y -see 12). Here the pair ranks R e determine how many spatial orbitals are assigned to each edge -the H2/STO-3G(2,4) system with circuit 3 for example has a pair rank of 2 and the H2/6-31G (2,8) with circuit 12 has a pair rank of 4. A minimally correlated approach (employed within most examples in this work) is to set R e = 2 for all edges. Like most variational quantum approaches the separable pair approximation described in Heuristic 1 is not an exact algorithm as successful convergence to the best solution for a given graph depends on the initial guess for the optimal orbitals. In this work we explicitly exploit insight from the molecular graph in order to generate good initial guesses as demonstrated within the examples below. In particular the instruction assign orbitals to edge in Heuristic 1 can be realized with different strategies. In this work, we are assigning atomic orbitals to the corresponding vertices v. In the case of hydrogenic systems, this is straightforward, as there is just a single vertex for each atom, meaning that for every edge e = (v k , v l ) we assign the atomic orbitals centered at the atoms corresponding to the two vertices that form the edge. For other systems, such as BeH 2 we need to decide which atomic orbital to map to which vertex, as for example the Be atom will result in two vertices. This will be illustrated explicitly in the next section.
As it was shown in [12] the SPA wavefunctions are classically tractable with linear memory requirement with respect to the number of qubits. This can be seen from Eq. (13) that describes a single electron pair in R e spatial orbitals with R e coefficients and Heuristic 1 that constructs an SPA wavefunction through |G| electron pairs, with |G| denoting the number of edges in the graph. The total wavefunction can therefore be represented by at most e R e ≤ max (R e ) |G| real coefficients -for a minimal correlated orbital set that would for example be 2|G|. If second order orbital optimization is employed the overall procedure stays classically tractable, with the orbital optimization scaling at most quartic with system size. [45,46] The numerical results of Ref. [12] indicate that Graph-SPAs are a good approximation of the wavefunction if the chemical graph is dominant in its importance. The importance of chemical graphs usually cannot be quantified explicitly, but they can be ordered in their relative importance based on chemical intuition and arguments from Lewis theory (see Sec. 3 for explicit examples). It might be argued, that in fact the SPA wavefunction quantifies the graph (e.g. through its fidelity with respect to the exact ground state, or via the energy expectation value) and not the other way around.

Quantum Circuits from multiple Graphs
As SPAs are classically tractable, it is clear that they cannot represent all molecules with arbitrary precision. Obvious candidates are molecules with multiple chemical graphs necessary for a satisfactory description. In Sec. 3 explicit examples for such systems are given. One way to incorporate more than one chemical graph in the design of a quantum circuit is by rotating the orbitals that represent one graph into orbitals that represent another graph, apply electron correlators, and rotate back into the original frame. In a way this end for construct initial guess for optimal orbitals initialize Hamiltonian H while energy is changing do θ * ← arg min θ ⟨H⟩ U (θ) optimize orbitals update Hamiltonian end while return U , H, θ * can be seen as an autoencoder [47] with the orbital rotators as encoder and the corresponding adjoint as decoder. The main difference in spirit to [47] is that the goal here is not compression and that there is no "training" involved. The procedure is formalized in Heuristic 2 and its success depends on the right choice and initialization of orbital rotators, which will again be tackled through insight from the molecular graph. Although there is also a heuristic choice in electron correlators, in practice it often suffices to use electron pair-excitations (the simplest electronic correlators) with angles initialized to zero. In the next section we will walk through some explicit examples how Heuristic 2 can be applied in practice and how successful circuits can be transferred between different molecules with similar graph representation.

Heuristic 2 Multi-Graph Circuit
Require: : Chemical graphs G SelectG=(V, E) ∈ G Construct U SPA and H according to Heuristic 1.

Examples
In the following we will apply the concepts discussed in the previous section to explicit examples beyond (effective) two-electron systems. This will demonstrate explicitly how suitable orbital rotations and useful initial guesses can be constructed by simple intuition drawn from the chemical graphs of the molecules. We start with small strongly correlated multi-electron systems H 4 and H 6 that are challenging for classical methods and finally show how the circuits constructed for H 4 can be transferred to other effective four electron systems like BeH 2 and the π-system of C 4 H 6 .

Linear H 4
A simple example for a electronically non-trivial system is the linear H 4 where four hydrogen atoms are placed on a straight line with 1.5 Å inter-atomic distance .
The most intuitive chemical graph is to view the system as two hydrogen molecules .
We will now represent this chemical graph with a quantum circuit in a minimal orbital basis, i.e. H 4 /STO-3G (4,8), obtained by placing a single s-type orbital on each hydrogen atom, in the same way as in (1). A quantum circuit that represents this chemical graph on eight qubits is build from two separated circuits that represent the individual H 2 bonds where the individual U 4 2 circuit are identical to (3) that described an individual hydrogen molecule. As before, the atomic basis functions need to be rotated to optimal linear combinations for the wavefunction created by the U (4,8) SPA circuit. An intuitive initial guess is to rotate them to localized orbitals that represent individual hydrogen molecules as in (8). Starting from (orthonormalized) atomic orbitals, this rotation acts as With this initial guess for the orbitals, the wavefunction prepared by U (4,8) SPA with optimized angle has an error of 40 millihartree in energy with respect to exact diagonalization. If we further optimize the orbitals through standard procedures the error can be reduced to 16 millihartree (denoted as SPA in Tab. 1) with optimal orbitals that look like (20) (see the appendix for the explicit coefficients). As before we can generate a transformed Hamiltonian with the optimized orbitals instead of including the rotation in the circuit (see appendix B and Sec. C.a of Ref. [12] for details on the implementation), so in the following the orbitals in (20) will be our basis. Due to the equidistant placing of the hydrogen atoms, the chemical graph (17) where we could also connect the outer atoms to form a spatially separated singlet pair similar to the solution of the stretched H 2 discussed in the previous section. We will chose not to do so in order to form a shallower and more localized circuit. The strategy is now, to rotate the optimized orbitals in (20) to a suitable representation for the chemical graph in (21) with the corresponding unitary U R , then correlate the orbitals that describe the central bond (between orbitals 1 and 2) with a unitary U C = U C 2 1 , and finally rotate back to the original frame with U † R . The full circuit, simply abbreviated as SPA+ is then which is the same circuit as shown in Fig. 1. We use a pair-restricted double excitation that transfers spin-paired electrons between two spatial orbitals (e.g. from orbital 1 to orbital 2) denoted as where the circuit represents an optimized form of the Jordan-Wigner encoded unitary e −i θ 2 G with the generator G = i a † 2 a 3 a † 4 a 5 − h.c. . Note, that the correlator type is independent of the chemical graph. In this case the pair-restricted double excitations were chosen since they are the cheapest fermionic correlators that have been applied in several other works. [12,13,25,48] In order to get a good initial guess for the parameters of the rotation circuit U R we split its construction into two parts: First (approximately) rotating back to the local frame U † R 0 and second forming the desired orbitals with an initial guess formed from the atomic orbitals through the unitary U R 1 . This corresponds to the following pattern with right-hand-sides defined as in Eq. (7). With this strategy we observe rapid convergence and further reduction in the energy error to 8 millihartree displayed in Tab. 1 along with a small comparison to the prominent k-UpCCGSD [13] approach that uses the same building blocks. Here, we used fixed initial angles (set to zero) for the k-UpCCGSD optimization in order to have a comparable and reproducible setting. The SPA and SPA+ circuits outperform the k-UpCCGSD flavors in all metrics (energy, fidelity, number of parameters, number of cnots, circuit depth, runtimeindicated with iteration counts). Another interesting effect can be observed in the fidelities with the H 4 ground state which is better for UpCCD compared to UpCC(G)SD without being reflected in the corresponding energies. The reason for this being overlaps with high energy excited states in the UpCCD wavefunction that increase the expected value of the energy. The fidelity drop from UpCCD to UpCC(G)SD witnesses the increased complexity of the optimization landscape that make the H 4 system hard to converge without exploiting its chemical structure. Note that UpCCGSD and 2-UpCCGSD can achieve higher accuracies in energy and fidelity for better initial angles. Repeated optimization runs with randomized initial angles can increase the accuracy, UpCCGSD is however not able to come close to the classically tracable SPA while 2-UpCCGSD can achive better energies and fidelities as SPA+ with significantly increased computational cost and an unreliable procedure. ADAPT-VQE with an operator pool build from the UpCCGSD operations achieves the same accuracy as 2-UpCCGSD with a more efficient circuit at the price of more iterations in the parameter optimization and 12 operator screening rounds. Despite having access to the same set of operations as SPA and SPA+, ADAPT(UpCCGSD) does not end up with a comparable circuit with respect to all metrics. An intuitive explanation is, that the adaptive procedure cannot locally detect the U R U C U † R motif as, for example, U R alone will not improve the energy, let alone parts of U R . In the same manner, U C has not the same effect without the basis change before and after.

Linear H 6
A similar system to the linear H 4 is the linear H 6 with the same inter-atomic distances of 1.5Å. We will now construct circuits that prepare good approximations for the groundstate of H 6 /STO-3G(6,12) by applying the same strategy as for the H 4 /STO-3G (4,8). The two main chemical graphs are in this case (25) and the corresponding SPA+ circuit is constructed as U 12 6 = U 4 2 (θ a ) ⊗ U 4 2 (θ b ) ⊗ U 4 2 (θ c ) with orbitals that resembles the first graph in (25). As before further simplifications based on symmetry (θ a = θ c ) could be made but will not be included in the analysis of the associated computational cost. In order to represent the second graph in (25) we follow the same strategy as before by constructing a orbital rotation unitary U R in two parts: One rotating to a localized frame followed by a second one that forms the orbitals which represent the central H-H bond. The rotation is again followed by a correlating unitary U C given by a pair-excitation on the central H-H followed by a rotation to the original frame. The so constructed SPA+ circuit is shown in full in (44) in the appendix and the obtained results are given in Tab Other than the rather artificial H 4 toy model, the BeH 2 molecule and the π-system of C 4 H 6 come closer to real-life examples. In the following we will see how they can be treated analog to the H 4 system by using the same circuits and intuition. We start by looking at the chemical graphs that represent the BeH 2 system (26) where we have depicted the central beryllium atom with two vertices, representing the two valence electrons. The first graph in (26) represents the molecular system with two Be-H bonds and the second the isolated atomic systems; their relative importance will vary with the Be-H inter-atomic distances. Note that the graphs are identical to the H 4 graphs in (21) and (17). Using the frozen-core approximation and a minimally correlated orbital basis -meaning that we use the same number of spatial orbitals as we have active (i.e. non-frozen) electrons, we arrive at an eight qubit representation of BeH 2 /STO-3G (4,8). In this case, we can remove the p y and p x orbitals form the STO-3G set -assuming the Be-H bonds are aligned in z direction), as they will not contribute to the optimized orbitals due to symmetry reasons. As for H 4 the remaining orbitals are then optimized with respect to the U (4,8) SPA wavefunction. A suitable guess for the molecular and atomic case is , (27) where s and p z orbitals are depicted in equal weighted superpositions with the atomic case on the right and the molecular case on the left hand side of (27). The beryllium atom is the green circle and the shades of blue and red in the orbitals again denote positive and negative regions of the orbitals. The optimized coefficients can be found in the appendix.
In the same way, we can treat the π-system of C 4 H 6 . This denotes the four electrons in orbitals spanned by the 4 p z (assuming the molecule lies in the x, y plane) basis function of the   [12] while ADAPT uses standard compilation. The ADAPT(UpCCGSD) circuits in particular are therefore expected to have lower cnot counts and depths as similar optimization can be applied to parts of the circuit. Note that circuit-depths can be further reduced through more efficient compilation of the orbital rotations [49,50]   STO-3G set. The other 26 electrons doubly occupy their corresponding frozen orbitals determined through Hartree-Fock. We will denote this active space of the so called π-system as C 4 H 6 /STO-3G (4,8) where the optimal π-orbitals are, again, determined through the U (4,8) SPA wavefunction and will look similar to the H 4 orbitals (20) just with p z -orbitals instead of s-orbitals, so that the identidcal initial guesses can be used for the two systems. The coefficients of the optimized π-orbitals can be found in the appendix. Note that those π-orbitals differ from the Hartree-Fock description of the π-system. In Fig. 3 we show the performance of identical circuits on H 4 /STO-3G (4,8), C 4 H 6 /STO-3G (4,8) and BeH 2 /STO-3G(4,8) with equilibrium (1.5Å) and stretched (3.0Å) Be-H bond distances. As circuits we use the SPA circuit (that also determines the optimized orbitals) in (3), the SPA+ circuit displayed in Fig. 1 and two extensions where we replace the single central pair-correlator (23) by a more general block U 4 C containing four pair correlators with individual parameters which are then placed before and after the rotations. In a second step we also improve the freedom in the orbital rotations by replacing the sequences in Fig. 1 a more flexible block U RR , which we abbreviate in Fig. 3 with RR while using R for the original sequence in Fig. 1. The circuit templates for the modified correlator and rotation blocks look like This allows for a systematic improvement of the energy below the millihartree threshold for all four systems.

Limitations
As Heuristic 2 is relying on chemical insight in order to provide important chemical graphs, it is clear, that the approach will not be able to describe all molecules with sufficient accuracy. Besides from transition metal compounds -which are fairly unexplored in the context of variational quantum algorithms -a simple example is the nitrogen molecule N 2 where only one reasonable chemical graph can be constructed, but the corresponding Graph-SPA energy still differs significantly from the exact ground state energy even at the equilibrium bond distance. This holds true for a minimal basis (STO-3G) as well as for directly determined MRA-PNOs and can be observed already in the (6,12) active space that freezes out the lone-pairs and includes only the triple bond as it was shown for N 2 /MRA-PNO(6,12) before [12]. It can be speculated, that the energy can be brought down further with the right intuition about graphs and corresponding orbitals. Indeed a successful trial could be achieved reducing the energy error of N 2 /STO-3G(6,12) close to equilibrium bond distance from around 40 millihartree (single graph SPA) to 3 millihartree (using two additional graphs that mix the π and σ bonds equipped with additional correlators similar to Fig. 3). This is however far from a well defined procedure and at this point not suitable for automatization which is why N 2 /STO-3G(6,12) is considered a challenging system in the scope of this work. It is reasonable to assume that this will extend to The used energies are from quantum circuits that are identical for all three systems. The SPA circuit is given in Eq. (18) and SPA+ denotes the circuit in Fig. 1. The other two variants follow the same pattern as in Fig. 1 but use more general, locally restricted, correlator blocks abbreviated as C, and extended rotation layers abbreviated as RR, both shown in (28).
other systems with triple-bonds and more challenging cases like the hextuple bond of the chromium dimer [51]. Combined with circuits obtained via global-optimization of UCC type operators [52] -recently applied to N 2 -better heuristics might be discovered in the future.
Although well-suited for state-of-the-art research in variational quantum circuit design the techniques described in this work are currently restricted to neutral molecules with even numbers of electrons. Extensions that go beyond this restricted formulation can however be envisioned. One way forward could be through "restriced-open" SPAs, containing a single unpaired electron, in their formulation as starting states. In appendix G we give an explicit example.

Computational Details & Code Availability
The examples in the last section were computed using tequila [53,54] with optimized gradients and quantum chemistry framework described in Ref. [28]. We used pyscf [55,56] as backend for Gaussian integral evaluation and orbital optimization with more details on the latter described in Ref. [12]. Direct basis determination via MRA-PNOs was performed with madness [57] 3 as described in Ref. [39] with the diagonal approximation described in Ref. [12] and the MP2-PNO implementation of Ref. [40] without cusp regularization. The latter interfaces mean-field implementations described in [58,59]. Quantum circuits were processed and simulated automatically through tequila using qulacs [60] as backend and the Jordan-Wigner encoding implemented within openfermion [61]. Optimization of the circuit parameters was performed with the automatically differentiable framework [28] within tequila using jax [62] and scipy [63]. The tequila library is available on github and an explicit code example is provided in the appendix.

Conclusion & Outlook
In this work molecular circuit design heuristics based on chemical graphs were developed. The heuristics show a good balance between hardware-efficiency and physically motivated design principles by producing relatively shallow and low-parametrized circuits with good convergence properties. This provides a way forward in all mayor challenges of variational quantum circuit design: Initial state generation by reducing the first graph to a classically tractable separable pair approximation. This is furthermore used to get an improved set of orbitals (compared to standard Hartree-Fock). Operator ordering and parameter initialization is achieved through chemical insight guided by the structure of the chemical graphs. So far, the molecular circuits were constructed manually by selecting graphs, orbital rotators and the corresponding initial orbital guesses through chemical intuition. The construction therefore still contains a relatively high human factor and is potentially challenging without some background in electronic structure. The original work on SPA circuits [12] used directly determined pair-natural orbitals [40] that can be computed through a blackbox procedure without background knowledge. Equipped with further heuristics, the techniques from this work could lead to an automatized construction of molecular quantum circuits and orbitals in a similar manner. In combination with the basis-set-free framework of [39], existing automatization protocols [64], and explicitly correlated corrections [65] a path towards an automatized black-box approach with accurate numerical precision could be envisioned. In particular, future heuristics could benefit from modern correlation measures [66][67][68] in order to effectively place electronic correlators in the spirit of [69]. An interesting direction towards improved heuristics could also be towards machine learning approaches, as they are for example developed to generate molecular graphs with specific properties [70], in order to generate and rank important graphs for a given set of molecular coordinates as well as to generate suitable initial guesses for orbital rotations. The combination with recent developments of circuit re-compilation [71,72] offers an interesting path towards hardware-adapted circuits with improved properties.
in Fig. 1, and in 44) do not need to be applied as quantum gates necessarily. In Heuristic 1 this is already used by directly solving for the optimized orbital rotation angles θ through a second order expansion of U † HU where U contains all orbital rotations. We refer to Ref. [12] for more details on the implementation used in this work and to Refs. [45,46,78] (in particular the appendix of [46]) for other examples of orbital-optimized circuits using the same technique.
In Jordan-Wigner encoding the two orbitals are mapped to qubits 2i and 2j and there spin-down counterparts to qubits 2i + 1 and 2j + 1 respectively. An individual spin-orbital rotation between orbitals ϕ 0 ↑ (mapped to qubit 0) and ϕ 1↑ (mapped to qubit 2) in Jordan-Wigner representation takes the form of two three-qubit Pauli rotations which can be graphically represented by and where each multi-Pauli rotation can be decomposed like where the first X and Y rotations have static angles of π 2 and the trailing ones − π 2 . For more details on the decomposition of Pauli rotations into standard gates see the review [11] and Ref. [53] for the corresponding implementation within the tequila framework.
In the same manner, the spatial orbitals can be transformed by applying the spin-up and spin-down transformations in sequence i ↓ with their respective angles being identical.

C Details on the frozen-core errors of LiH
In Fig. 4 we show the frozen-core errors of LiH in the minimal basis STO-3G, a larger basis with specialized core-correlation functions cc-pCVDZ and a customized MRA-PNO basis that was constructed to include correlation in the core electrons (by optimizing one function for each PNO-MP2 pair [12,39,40]). The STO-3G basis is missing core polarization functions, so that there is no notable representation of that correlation in the basis, leading to a frozen-core error below the millihartree threshold for all bond distances. Here the frozen-core errors is identical to the SPA/STO-3G(2,12) error in Fig. 2 of the main text. The cc-pCVDZ basis admids a significant frozen-core error, due to the presence of specialized core-polarization functions in the basis. The cc-pCVDZ basis is inconvenient for the investigation of the SPA performance as it would require 40 qubits to represent. This is why we employed a system adapted MRA-PNO basis in Fig. 2 in order to demonstrate that the LiH molecule can be efficiently represented by separable pairs. As the frozen-core error for the MRA-PNOs is significant and comparable to cc-pCVDZ, the SPA wavefunctions  includes two electrons pairs. The opt-SPA results is again below the millihartree threshold for all bond distances. This shows, that, despite having significant core correlation in the system, the corresponding wavefunction is still separable into two electron pairs. In chemical terms this means, that the core electrons are not correlated with the valence ("bonding") electrons in a significant way.

D Optimized Orbitals
In the article we represent orbitals with small cartoons that depict them as linear combinations of their respective atomic basis orbitals. The optimized orbitals for H2/STO-3G(2,4) were for example shown as with the two atomic orbitals on the left, and the two optimized combination son the right. The corresponding matrix, holding the coefficients that define the optimized orbitals in a linear superposition of atomic basis orbitals is and U R represents the corresponding quantum circuit (7) that has the same effect (see Sec B) when acting on an Hamiltonian H formed from orthonormalized atomic basis orbitalsi.e. H ′ = U † R HU R is the same as a qubit Hamiltonian formed from orbitals obtained with the matrix above. In the following the optimized coefficients with respect to Graph-SPA circuits are displayed as linear combinations of their respective (non-normalized) atomic orbitals of the corresponding STO-3G basis set.

D.1 Optimized H 4 Orbitals
Coefficient matrix for the optimized orbitals of the linear H 4 system represented in (20) Rows represent the 5 optimized orbitals with the first row beeing the frozen core orbital. The columns denote the atomic orbitals from the STO-3G set in the order: Be-S, Be-S, Be-P z , right-H-S, left-H-S. The initial guess for the optimization where equal weighted superpositions depicted in (27) with no contribution of the first Be-S. The initial guess of the core orbital (not depicted in (27) was only the first Be-S orbital.

G Comment on systems with odd electron number
The main article focuses primarily on systems with even electron number, which are currently the dominating usecases in comparable literature. Extensions of the current framework to odd-numbered electrons is possible, but will require some work in the underlying details. In the following a sketch for a potential approach will be given, using the neutral H 3 /STO-3G(3,6) molecule -a 3-electron system in 6 qubits -as example.
This work used SPA circuits in order to represent the first molecular graph. SPAs were developed for even-numbered electrons as well, and although we could extend them to oddnumbered electrons by simply replacing one of the electron-pairs with a single electron, we will chose a simpler strategy where all graphs are represented with the same U R U C U R motif. As an initial state (which in the main article is included in the SPA construction), we then chose a basis state with the right amount of electrons. In the H 3 case this can for example be with two spin-paired electrons in the first, and a single electron in the third orbital. On this initial state we are then applying unitaries that correspond to specific molecular graphs (as in the main text) |Ψ(a, b, c⟩ = U (c)U (b)U (a) |Ψ 0 ⟩ with the graphs being (47) In Tab. 3 we show the performance of such circuits where we, as in the main article, started from atomic orbitals and either used a single orbital rotator per graph (corresponding to a localized orbital picture) or allowed more freedom by placing three rotators per graph (corresponding to full orbital optimization). Similar to the H 4 results in the main text we see that graphs a and b are dominant and lead to exact results when the circuits are equipped with enough freedom in the orbital rotations.