Geometric quantum adiabatic methods for quantum chemistry

Existing quantum algorithms for quantum chemistry work well near the equilibrium geometry of molecules, but the results can become unstable when the chemical bonds are broken at large atomic distances. For any adiabatic approach, this usually leads to serious problems, such as level crossing and/or energy gap closing along the adiabatic evolution path. In this work, we propose a quantum algorithm based on adiabatic evolution to obtain molecular eigenstates and eigenenergies in quantum chemistry, which exploits a smooth geometric deformation by changing bond lengths and bond angles. Even with a simple uniform stretching of chemical bonds, this algorithm performs more stably and achieves better accuracy than our previous adiabatic method [Phys. Rev. Research 3, 013104 (2021)]. It solves the problems related to energy gap closing and level crossing along the adiabatic evolution path at large atomic distances. We demonstrate its utility in several examples, including H${}_2$O, CH${}_2$, and a chemical reaction of H${}_2$+D${}_2\rightarrow$ 2HD. Furthermore, our fidelity analysis demonstrates that even with finite bond length changes, our algorithm still achieves high fidelity with the ground state.


I. INTRODUCTION
In recent years, quantum chemistry [1] has emerged as a promising domain science field where quantum computers can potentially lead to a breakthrough [2]. Even for ground-state total energy problems, many systems exhibit strong character of many-body correlation effects (e.g., bond-stretching and transition states in chemical reactions), which are beyond the scope of widely used mean-field methods. As wave-function-based correlated methods remain computationally intractable against system size on classical computers, various quantum algorithms, such as iterative quantum phase estimation (iQPE) [3], variational quantum eigensolver (VQE) [4,5], quantum adiabatic evolution (QAE) [6,7], and quantum annealing [8,9], have been proposed to tackle these quantum chemistry problems. Several quantum algorithms have been implemented and demonstrated for small molecules [2,[4][5][6][9][10][11][12]. Further improved methods such as the equation of motion [13] and the adaptive VQE [14] have also been designed. Some of us have previously proposed an adiabatic method that uses a particular interpolation of two Hamiltonians [15]: (1) the final one being the targeted full many-body molecular Hamiltonian and (2) the initial one being the maximally commuting (MC) Hamiltonian constructed from the final Hamiltonian. Such a quantum adiabatic evolution from the initial MC Hamiltonian (MC-QAE) approach works well for molecules near their equilibrium position, but the results become inaccurate at large atomic separations, due to energy crossing and the presence of * dlu@bnl.gov † tzu-chieh.wei@stonybrook.edu dense low lying levels. Previously we proposed a quantum Zeno approach to overcome this issue, which uses projection to instantaneous eigenstates of the discretized time-dependent Hamiltonian [15]. However, the spectral projection in the quantum Zeno approach at the present is less practical to implement than discretized Trotter evolution of the adiabatic Hamiltonian.
Here, we propose an adiabatic evolution method to compute ground state and low-lying excited state energies of the many-body Hamiltonian of molecular systems along an adiabatic geometric path, which we refer to as the geometric quantum adiabatic evolution (GeoQAE). The GeoQAE starts with an initial geometry where MC-QAE can solve the ground state and low-lying excited states accurately. Then the system evolves adiabatically following a smooth geometric deformation by, e.g., stretching or shrinking bonds and/or increasing or decreasing bond angles. We show that the GeoQAE approach does not suffer from the energy level acrossing issue in the MC-QAE approach at large atomic distances. We demonstrate its utility by computing the energies of the ground state and low-lying excited states of exemplary molecular systems, including H 2 O, CH 2 , and the potential energy surface of the chemical reaction of H 2 +D 2 → 2HD.
In this work, we found an important property of the MC Hamiltonian, proposed in our earlier work [15], that in the qubit basis the MC Hamiltonian is equivalent to keeping only diagonal terms of the full configuration interaction (CI) Hamiltonian. The MC Hamiltonian in fact is closely related to the Hartree-Fock approximation, because under the fermionic basis it consists of diagonal one-body terms and two-body Hartree and Fock terms. We show that the Hartree-Fock ground state is an eigenstate of the MC Hamiltonian and the Hartree-Fock ground-state energy is the corresponding eigenenergy.
The remaining of the paper is organized as follows. In Sec. II, we briefly review the setup in the molecular Hamiltonian and discuss our choice of the initial Hamiltonian in the adiabatic evolution. In Sec. III, we present our geometric adiabatic approach and describe the new protocol. There, we also present the analysis of the gap in the new adiabatic Hamiltonian to justify the approach. In Sec. IV, we show the results by the new GeoQAE approach applied to two molecules, H 2 O and CH 2 , yielding better results than the previous MC-QAE approach and the other quantum Zeno approach. In Sec. V, we present a study on a simplified chemical reaction of H 2 +D 2 → 2HD, via the ground-state energy versus the bond lengths. In Sec. VI, we discuss the errors and fidelity from the perspectives of the choice of the geometric path, the evolution time T and the discrete step number M in an evolution. We make some concluding remarks in Sec. VII.

II. MOLECULAR HAMILTONIAN AND THE ITS MAXIMALLY COMMUTING HAMILTONIAN
In this study, we focus on finding the ground-state energy and low lying excited states of a molecule at fixed atomic coordinates. For a given set of spin-orbitals, the many-body Hamiltonian can be written in the second quantized form, where i, j, k, l label the spin orbitals. The one-body hopping t ij terms and two-body interacting terms u ijkl are given by the following expressions, where Ψ i (x 1 ) are single particle spin wave functions of orbitals i. Note that [ij|kl] is expressed in the so-called chemists' notation and the same quantity is denoted as ik|jl in the physicists' notation [1]. These coefficients are calculated using the standard quantum chemistry package, PySCF [16], written in the Python language. The choice of the basis set needs to be made to compute Eqs. 2. Large basis sets can in principle lead to better the numerical convergence with respect to the basis set size, however they require more qubits to implement the quantum algorithm. In this work the Slater-type orbital (STO)-3G basis (see e.g. [1,17]) is used as a compromise between accuracy and computational cost, while our main conclusion does not rely on the choice of the basis set.
The goal is to find the eigenstates and eigenenergies of Hamiltonian H of molecules, i.e., H|ψ E = E|ψ E , in particular the ground state |ψ G and its energy E G , as well as low lying states. In order for the problem to be solved on quantum computers, we need to convert the above fermionic Hamiltonian to one composed of qubits. Existing methods include the Jordan-Wigner, parity, Bravyi-Kitaev, and superfast Bravyi-Kitaev transformations [18][19][20]. By using any of these methods, we can transform the fermion operators into Pauli operators, where P i 's are n-qubit Pauli operators and h i 's are the corresponding coefficients. For the results presented below, we use the parity and Jordan-Wigner transformation methods (see Appendix A).

A. Initial Hamiltonian and initial state
We will first use the Hatree-Fock procedure, briefly reviewed in Appendix B, which uses an initial set of atomic orbitals with the full Hamiltonian, such as shown in Eq. (1) and aims to obtain a new set of converged molecular orbitals, particularly indexed below by Greek letters α, β, γ, etc. We will emphasize the molecular basis by adding the superscript MO to the hopping and interaction coefficient in the resultant full Hamiltonian as t MO and u MO , respectively, The resultant Hartree-Fock state is then a single-Slater determinant in this molecular basis. For convenience, we denote the qubit version of H F full as H P full , where the subscript P means Pauli operators.
Our choice of the initial Hamiltonian H P I (in the MC-QAE approach) is the diagonal part of the qubit version of the full Hamiltonian H P full , which only consists of product of Z-terms and identities. For transformations mentioned above from the fermionic to the qubit representation, the corresponding eigenstates of H P I are in the computational basis, which corresponds to singledeterminant states in the fermionic picture. In the fermionic Hamiltonian, one Slater determinant is the simultaneous eigenvector of the corresponding number operators a † α a α for all α. If the transformation converts a † α a α into Pauli Z terms, and a † α a β non-Z terms for β = α, then the diagonal part of H P full is the summation of all possible combination of a † α a α , a † α a † β a β a α and a † α a † β a α a β . Therefore, the fermionic version of our initial Hamilto-nian is where the superscript F is added explicitly to emphasize the H I in terms of fermion operators. We remark that when written in the Hartree-Fock molecular basis, the Hartree-Fock ground-state wave function is an eigenstate of H F I with eigenenergy E HF . This can be easily verified by directly applying H F I to the Hartree-Fock ground state |ψ HF = |1..10..0 by filling up orbitals with lowest energies, and obtaining E HF to be the Hartree-Fock ground-state energy, αβγδ is the two-electron integral of spin-orbitals calculated in the Hartree-Fock molecular basis. More details are shown in the Appendix C.
We remark that the Hartree-Fock ground state is also an eigenstate of the so-called Fock operator (a.k.a. the Hartree-Fock Hamiltonian) F in Eq. (B7),but with an eigenenergyẼ different from the Hartree-Fock groundstate energy [1], due to double counting in the interaction energy (i.e. a factor of two difference in the second term), The Hamiltonian form (4), when converted to the qubit basis (denoted by H I ), is the one found from the previous work [15] constructed by finding the maximum commuting set of terms in the Hamiltonian. Prior results [15] show that for some molecules if one chooses a bond length near the equilibrium distance and uses the H P I as the initial Hamiltonian, it can drive the Hartree-Fock initial state to the final ground state of the full Hamiltonian via adiabatic evolution interpolating the two Hamiltonians, Below we shall see that the Hartree-Fock state does not always lead to the final ground state. However, we can choose other computational basis states or their superposition as the initial states. Note that the electron number is preserved during the evolution. So it is possible to choose an appropriate initial state to target the certain sectors with desired the electron number. In some particular cases where the system has additional symmetries and the direct evolution fails, we may use a symmetrized initial state instead to overcome the problem; see Sec. IV B below. For example, supposing the MC Hamiltonian has two degenerate ground states |ψ HF ≡ |ψ 1 = |01 , |ψ 2 = |10 and we have some prior knowledge that the final ground state is a singlet state, we can modify our initial state according to Due to the symmetry preserved by both initial and final Hamiltonians, choosing an initial state sharing the same symmetry as the expected finial ground state will restrict the evolution space within the subspace of desired symmetry. Thus, the evolution will likely end up with the final ground state. We note that the above state ψ 0 in the qubit basis can be initialized by a short-depth circuit.

III. GEOMETRIC ADIABATIC PATH AND THE NEW PROTOCOL
To overcome that problem that the MC-QAE approach fails at large bond lengths despite that it works well near the equilibrium distance, we propose a new idea, justify it and examine it in several examples. Suppose we have the qubit Hamiltonian H P (r 0 ) (and have obtained its ground state, e.g. using the previous MC-QAE method that works) for a molecule at r 0 near the equilibrium position. In order to obtain the ground state and energy of H P full (r) at large molecular distances, we construct a series of subsequent Hamiltonians H P full (r 1 ), H P full (r 2 ), ..., H P full (r N = r) so that the sequence r 0 , r 1 , . . . , r N represents a gradual change in the atomic positions. Then we let the system evolve piece-wise according to their interpolation as follows, For each step, the evolution from H P (r i ) to H P (r i+1 ) is likely to succeed due to the continuity with r of the Hamiltonian H P (r).

A. Gap analysis
To justify and illustrate that this alternative approach works, we compare the spectrum of two path dependent Hamiltonians leading to the same final Hamiltonian The second Hamtilonian begins with the full Hamiltonian at an earlier distance, H P (r = 1.758Å) and ends with the final Hamiltonian H P (r = 1.958Å), (12) Their spectra as a function of s are shown in Fig. 1. There are energy crossings in H 1 (s), as shown in Fig. 1a, and this is the reason that the direct MC-QAE fails. However, assuming one can arrive at the ground state of a Hamiltonian at, e.g., a shorter distance r = 1.758Å, one can use this Hamiltonian as the initial one to adiabatically evolve its ground state (as well as excited states) to that of H P full (r = 1.958Å) or one at a shorter distance. The resultant path-dependent Hamiltonian H 2 (s) has energy levels smoothly connected without any crossing, as shown in Fig. 1b. This is perhaps expected, as it corresponds to a small stretching of the chemical bonds. Thus, we see that the energy-crossing problem is solved by evolving from the nearby bond length r = 1.758Å instead of the maximum commuting Hamiltonian at r = 1.958Å.
The above is not an isolated example. In fact, for most cases, the energy gaps between levels during the evolution from H P full (r i ) to H P full (r i+1 ) remain finite and large, even in the region of large bond lengths where direct evolution (from the MC Hamiltonian) fails. This consequence relies on the continuity in the Hamiltonians when atoms are displaced gradually in changing their geometry. (See Appendix D for several methods we use to ensure continuity in the Hamiltonians.) For the cases where the Hartree-Fock method converges, the changing of oneand two-body coefficients is small when ∆r = r i+1 − r i is sufficiently small. This results in the separation of energy levels in the interpolated Hamiltonian. In our cal-culations, we found |∆r| ≈ 0.1r 0 is good enough for the evolving Hamiltonian to maintain gaps. However, when classical Hartree-Fock fails to converge, we may need to transform the Hamiltonian H P full (r i+1 ) back to the previous molecular basis {|ϕ µ ri }, in order to maintain numerical continuity of the coefficients. (As also seen below in Sec. VI, despite the small gap when ∆r is not small, the fidelity of the final ground state can still be large.) In our experience, by using the converged solution in a prior distance as the initial guess in the (classical) Hartree-Fock procedure usually results in the convergence, even if the random initial guess does not.

B. New protocol
In the following we assume that the Hamiltonians are all written in the converged Hartree-Fock molecular basis. With the assumption that the changing of the geometric configuration ∆r i = r i+1 − r i is much smaller than the atomic distance, we can ensure that the next Hamiltonian is not dramatically different from the previous one. We refer to this approach as the geometric quantum adiabatic evolution (GeoQAE) and the whole protocol can be summarized as follows.
1. Choose a proper bond length r 0 (near the equilibrium distance), construct a Hamiltonian of a molecule, and transform the fermion operators to Pauli operators; 2. Choose the all-Z terms to construct the initial Hamiltonian H P I (r 0 ) and the Hartree-Fock ground state or other appropriate computational states as the initial state ψ 0 ; 3. Discretize the time steps by t k = k M T and obtain a series of Hamiltonians H 0 (t k ) in Eq. 10; 4.
Perform the evolution operator e −iH0(t k )∆T on ψ 0 successively for k = 1, . . . , M , to aim at one eigenstate of the final Hamiltonian H P full (r 0 ); 5. Choose a suitable series of bond lengths r 1 , r 2 , ..., r N with r N equal to the desired bond length. Then repeat similar evolution procedure (3 and 4) with the initial Hamiltonian H P full (r i ) and final Hamiltonian H P full (r i+1 ) multiple times. We expect to obtain an approximated ground eigenstate of H P full (r N ). We can perform measurements to obtain its corresponding energy.
We remark that the above points 1-4 are exactly the MC-QAE approach, and the protocol for GeoQAE uses MC-QAE as its first step at an appropriate molecular distance r 0 (usually around the equilibrium position, e.g. such as determined the Hartree-Fock method) to arrive at the exact molecular eigenstates and subsequently undergoes the geometric change of the molecular configuration, as in point 5 above to reach the corresponding eigenstates.

IV. RESULTS ON MOLECULAR ENERGIES AND FIDELITIES
We apply our geometric adiabatic path approach to three different systems H 2 O, CH 2 , and a simple chemical reaction H 2 + D 2 −→ 2HD. We numerically calculate the following evolution, and evaluate the resultant energy, where ∆T = T /M and |ψ(0) is a suitably chosen initial state. The molecular Hamiltonian (i.e. its coefficients t's and u's) is calculated by the PySCF package [16] and transformed to Pauli operators by the Qiskit [21]. The transformation we employ in this work includes the Jordan-Wigner transformation and the parity mapping. We also freeze the 1s orbital of larger atoms, such as C and O. If not specified, we use the Hartree-Fock state as our initial state. The case of CH 2 requires the use of other initial states. We also compare our results by the GeoQAE method with those obtained with direct evolution (i.e. the MC-QAE) from H P I (r i ) to H P full (r i ). We now present several case studies.  In our GeoQAE calculation, we choose r 0 = 0.9584Å in the initial MC-QAE step with the Hartree-Fock state as the initial state to obtain the ground state of the molecular Hamiltonian at this position. Then for other positions, we follow the procedure outlined in Sec. III. Compared to the direct evolution from the MC Hamiltonian, the GeoQAE method gives very good results (e.g. with relative errors being 10 −5 or smaller) even in the large atomic distances where MC-QAE fails, as shown in   2. Nevertheless, due to the sequential evolution the error will accumulate from all steps of evolution. Therefore, to achieve the same fixed accuracy at larger atomic distances, we may need more discretization steps M and a larger evolution time T . But a recent study shows that the accumulation error in using Trotterization for the adiabatic evolution is not as severe as one would expect if the initial state is an eigenstate [22], as we have also observed in the bottom panel of Fig. 2.
Ground-state fidelity. We also compute the fidelity of the final evolved state |ψ(T ) with the exact molecular ground state |ψ g of H P , i.e. f (ψ g , ψ(T )) = | ψ g |ψ(T ) | [23]. In these simulations we use evolution time T = 40 and discretize the continuous time evolution to M = 20 time slices. We compare the resultant ground-state fidelity using two different evolutions, one with MC-QAE and the other with GeoQAE.
Even though the energy calculation by the MC-QAE (see Fig. 2) may only has at most 0.1% error in the range we consider, it is actually not targeting the ground state (for r > ∼ 1.5Å), but an excited state (or some combination of them) instead. As the result of the fidelity in Fig. 3 shows, the final state from MC-QAE at large r's has zero fidelity with the ground state. On the other hand the evolved final state via GeoQAE does have close to unity fidelity with the ground state of the final Hamiltonian. With the MC-QAE, for d > ∼ 1.76Å, the fidelity goes to nearly 0. This is because the direct adiabatic evolution from the MC Hamiltonian experiences level crossings and can only find excited states for large molecular distances d, which was noticed in our previous work [15] and illustrated in Fig. 1a. But with the geometric evolution, we can maintain the fidelity to nearly 1 even at large atomic distances. Alternatively, one may also explore other nearby atomic arrangement (to the target distance) as the starting point of the GeoQAE algorithm.

B. CH2
Energy of the lowest few levels and the choice of initial states. For the CH 2 molecule, we fix the H-C-H angle to be the equilibrium angle θ = 101.89 • and vary the C-H bond. The results of simulating our GeoQAE method are shown in Fig. 4. Here we freeze the 1s orbital of the C atom, similar to the H 2 O case, but use the Jordan-Wigner mapping to qubits. The final qubit number to represent the Hamiltonian is 12, as we do not have the qubit reduction as in the parity mapping.
We list the lowest 6 eigenstates of the MC Hamiltonian H P I (r 0 ) at r 0 = 1.1089Å,  In the above, the first 6 binary numbers in the kets represent the occupation of the molecular orbitals with spin up (usually labeled as α) in the order of increasing energy and the last 6 numbers represent the occupation of the same set of molecular orbitals but with spin down (usually labeled as β). The Hartree-Fock wave function turns out to be φ 6 , which is not the lowest eigenstate of the H P I . The wave functions φ 1 and φ 2 are triplets and degenerate ground states (w.r.t. H P I ). The wave functions φ 4 and φ 5 are also degenerate w.r.t. H P I but they differ by opposite spins in the last two occupied orbitals, i.e. ↓ 3 ↑ 4 vs. ↑ 3 ↓ 4 (where the subscripts indicate the molecular orbitals). Due to the fact that the MC and full Hamiltonians conserve the total spin angular momentum and its z component, we can consider combinations of two single-Slater determinants (having the same initial energy) to form singlet and/or triplet states of opposite spins as the initial states. For example, from φ 4 and φ 5 we can define a singlet ψ 3 ≡ (φ 4 − φ 5 )/ √ 2 and a triplet ψ 1,c ≡ (φ 4 + φ 5 )/ √ 2 wave functions for the initial states. For convenience, we also define ψ 1,a ≡ φ 1 , ψ 1,b ≡ φ 2 , ψ 2 ≡ φ 3 , and ψ 4 ≡ φ 6 .
The reason to define these four sets of wave functions ψ 1 's, ψ 2 , ψ 3 and ψ 4 is that their evolution under the combination of (1 − s)H P I + sH P full gives rise to the lowest four energy levels of the final Hamiltonian around the equilibrium position (with the ground-state degeneracy being 3). We have carried out this step of evolution (which is the MC-QAE procedure) numerically to verify that they indeed achieve the corresponding energies and levels with 99% fidelity or above. Subsequently, we proceed with the GeoQAE steps to obtain energies and wave functions at other molecular distances (both smaller and larger). Shown in Fig. 4, the three energy curves originating from ψ 1 's, ψ 2 and ψ 3 remain the lowest three energy curves, whereas the curve originating from ψ 4 evolves to a higher excited state at large molecular distances (as seen from its gap between the three lower curves). It is interesting that the lowest two energy curves cross, i.e. the three-fold degenerate ground states become the first excited states at larger molecular distances. Given that these different levels have different symmetries, the adiabatic evolution will not lead to any mixing.
We remark that the superposition states, such as ψ 1,c and ψ 3 , are similar to the Greenberger-Horne-Zeilinger states and can be easily created and initialized by a shortdepth circuit consisting of Hadamard and CNOT gates. In the above, φ i 's were calculated by direct diagonalizing the initial MC Hamiltonian, which is essentially classical (with Hamiltonian terms being products of Pauli Z and identity operators). Despite that finding the lowest configuration of a generic classical (e.g. Ising spin-glass) Hamiltonian can be NP, the time complexity of finding its low-lying states is still much smaller than those of the final quantum Hamiltonian for molecules. However, it is also not known whether the MC Hamiltonians from quantum chemistry problems are necessarily NP-hard.
Fidelity with lowest three levels and the other excited level. In addition to the energy curves, we also calculate the fidelity of the evolved states Ψ T with the corresponding eigenstates Ψ E (exactly solved numerically), f (ψ T , ψ E ) = | ψ E |ψ T |, and the results are shown in Fig. 5. Despite that the ground states (near equilibrium) are three-fold degenerate, their spin configurations are different and hence their states can be numerically separated and distinguished, the use of the fidelity expression f (ψ T , ψ E ) = | ψ E |ψ T | is appropriate. In case that the numerically solved degenerate states Ψ E,i can be arbitrary superposed, one can average over the degen- We can conclude from the fidelity results shown in Fig. 5 that the geometric evolution can find all the lowest three levels and the other excited level (the latter corresponding to the 4th level around the equilibrium) with high fidelity. In contrast, the direct evolution from Maximum commuting Hamiltonian (i.e. MC-QAE) encounters problem at larger distances (not shown explicitly) [15], similar to the case of H 2 O in the previous section.

V. CHEMICAL REACTION
The knowledge of the potential energy surface of a chemical reaction is critical to understand the reaction dynamics and kinetics. In particular, reaction energies and reaction barrier heights (i.e., the energy differences between reactants/products and the transition state) are key quantities that dictate the reaction energetics. It is important to develop quantum algorithms to compute the reaction potential energy surface.
In the previous section, we have seen that molecular energies can be accurately obtained by our GeoQAE, which extends the validity of the MC-QAE approach to a much greater extent. As a further application, we now apply this method to obtain the energy landscape of one chemical reaction. We use H 2 + D 2 −→ 2HD as an example. Such a system is also of interest because of experimental realization of cold controlled chemistry [24].
To reduce the complexity of the problem, we consider that the four atoms are initially set as the vertices of a rectangle with smaller horizontal edge; see Fig. 6. And the reaction is simplified to a process that the horizontal edge is getting longer, and it can be viewed effectively as the horizontal hydrogen molecules break up and form two new vertical molecules separately each with one deuterium. The highest energy of such process is when the four atoms become a square. Thus, in order to know the energy barrier of such reaction, it is crucial to know the accurate energy of the square configuration [25]. For those cases that are not square, we find that the MC-QAE method works fairly well. But when it comes close to a square configuration, the MC-QAE fails because of the energy crossing arising from the extra symmetry of a square. In this symmetric regime, the GeoQAE performs much better, as demonstrated in Figs. 7a and 7b.
In our calculations, we use the STO-3G basis for orbitals and the Jordan-Wigner mapping to convert the fermionic Hamiltonians to the corresponding qubit Hamiltonians. The total qubit number required is 8. For the off-diagonal points representing unequal distances in the two graphs, we use the MC-QAE for both cases. For the diagonal points, we use a two-step approach: first we let the system evolve to its nearest off-diagonal neighbor (via the MC-QAE), then do an adiabatic evolution to evolve to its ground state at the diagonal point via the GeoQAE, which is illustrated by the arrows in Fig. 8. For each run of evolution from one distance to another, we set T = 80 and M = 40.
The energy difference between the GeoQAE and the exact solution is shown Fig. 8 and the results are very accurate except in the region of larger r's, which can be improved by using larger M and T (see e.g. Sec VI). In Fig. 9, we explicitly compare the results of the MC-QAE and GeoQAE at r 1 = r 2 , i.e. the square configuration. It can be seen that the energies from the MC-QAE are visibly higher than those from the GeoQAE, with the latter having 10 −4 relative errors or smaller. This again shows the drastic improvement of the new approach by

VI. ERROR ANALYSIS AND COMPUTATIONAL COMPLEXITY
Apart from the possible noise and error from real quantum computers, which was discussed in our previous paper [15], the theoretical error of the geometric adiabatic path mainly depends on three parts: the choice of geometric path, the evolution time T and the discretizing step number M . We discuss these below.

A. Choices of geometric path
The geometric path is crucial for the whole procedure. While in principle when ∆r is small enough, the physical evolution is continuous, it is still important in practice to choose a proper ∆r that is neither too large that the adiabatic evolution fails nor too small that the computational time becomes too long and errors accumulate too much.
Here we explore the length of ∆r as an example. When ∆r gets larger, it requires fewer (accumulated) steps to evolve to our desired bond length. However, within each evolution from one distance to another, the minimum gap will get smaller. Thus, we need more evolution time for each step in order to achieve the same accuracy. As shown in Fig. 10, the energy gaps decrease exponentially with ∆r as expected. Despite this, the ground-state fidelity can still achieve a value of 0.93 with ∆r = 2.0Å with T = 40 and M = 20. This can be further improved to 0.99 larger parameter values T = 160 and M = 80 (four times more discretization steps). In addition, the geometric path itself also plays an important role here. We have only explored in this work changing few geometric parameter, i.e. one or two bond lengths. The vast freedom to choose the geometric path indicates potential advantage. However, there needs further investigation on how to design an efficient and practical geometric path for more sophisticated cases. One possible approach is to define a coordinate of collective variables and identify the difficult cases that need the geometric path method as critical points along the coordinate. For instance, the reaction coordinate is often used to study chemical reactions, and the transition state is the maximum (or saddle) point on a one-dimensional (or two-dimensional) potential energy surface. When we need the geometric path method to reach the transition state, a path along the reaction coordinate could be a good choice. But we leave these for future exploration. The two parameters T and M correspond to the the degree of adiabaticity and discretization approximating rate, respectively. If the adiabatic evolution does not fail (i.e. has a gap), we can theoretically achieve any accuracy by choosing sufficiently large T and M . However, such large choices are not necessarily practical due to accumulation of noisy gate errors. In our simulations, we usually choose T to be inverse proportional to the minimum gap of the evolution, and fixed the ratio of T /M to keep the dicreatization error small. (However, one may not necessarily have the information on the gaps, and may need to empirically test a few choices to see if the energy is lowered or converged as T increases.) We tested the influence of T and M on the fidelity with the ground state and the results are presented in the Fig. 11. In our study, we find that T has larger impact than M on the fidelity.

VII. CONCLUSION
We have substantially improved an adiabatic algorithm for molecular energies (MC-QAE) introduced in a previous work by utilizing a geometric adiabatic evolution path (GeoQAE). We showed by numerical simulations of the adiabatic evolution that this new approach gives significant improvement on both ground-state energies and the fidelity of wave functions (as well as lowest few excited states) for different molecules in the region where the MC-QAE fails due to the energy level crossing. We also demonstrated its potential application by simulating the potential surface for an exemplary chemical reaction involving two pairs of molecules. In our fidelity analysis, we find that our GeoQAE also gives high fidelity ground states using adiabatic evolution between two Hamiltonians with finite difference in bond lengths. The potential advantage of our method may rise from the vast freedom in choosing the adiabatic path.
The idea of choosing a suitable Hamiltonian that can be connected naturally and smoothly to the final Hamiltonian can be applied to other many-body problems. In our present focus on the molecular energy, the MC Hamiltonian associated with a particular final molecular Hamiltonian generalizes the Fock operator such that the Hartree-Fock ground state is an eigenstate with the eigenenergy being the Hartree-Fock ground-state energy.
Although on the present noisy quantum devices, our GeoQAE approach might not yet be implemented to yield accurate results due to the required large gate numbers for the evolution operators, hardwares are continuing to be improved and this may make our algorithm practical in the future. Moreover, the existence of gaps in the suitably chosen geometric path may allow alternative quantum approximate optimization algorithm (QAOA) like approaches or better variational ansätze to be developed. p i s. For example, the fermionic state |101110 f will be mapped to parity state |110100 p . Now consider a fermionic state with spin configuration |f α 1 ...f α n , f β n+1 ...f β 2n and its parity correspondence |p 1 ...p n , p n+1 ...p 2n . Note that the p n counts for the parity of all α-spin orbitals and p 2n counts for the parity of all orbitals. Therefore, if one further assume the conservation of the spin and particle number, the parity numbers p n , p 2n are constants. Thus, one can remove the two qubits and only consider the states in the subspace spanned by states with same electron number and spin configuration, which reduces the complexity of both quantum computations and classical simulations.
were introduced in the main text. The energy expression becomes In Hartree-Fock method, we assume the ground state wave function can be approximated by a single Slater determinant, and the Hamiltonian can be reduced to a one-body Hamiltonian, wheref is referred to as the Fock operator. In the Hartee-Fock calculation, we seek for an optimal molecular basis set (labeled by Greek letters in the main text) that can minimize the energy expectation value, starting from a chosen set of atomic orbitals (labelled by Roman letters), such as the STO-3G basis set. An iterative calculation can be performed to obtain the optimal molecular basis set. After we obtain one basis set for the Fock operator, we diagonalize the operator and rewrite the full Hamiltonian in the new basis set, then perform the diagonalization again until the result converges.
We use the quantum chemistry package PySCF [16] to perform this classical step of calculations.