Shallow-circuit variational quantum eigensolver based on symmetry-inspired Hilbert space partitioning for quantum chemical calculations

Development of resource-friendly quantum algorithms remains highly desirable for noisy intermediate-scale quantum computing. Based on the variational quantum eigensolver (VQE) with unitary coupled cluster ansatz, we demonstrate that partitioning of the Hilbert space made possible by the point group symmetry of the molecular systems greatly reduces the number of variational operators by confining the variational search within a subspace. In addition, we found that instead of including all subterms for each excitation operator, a single-term representation suffices to reach required accuracy for various molecules tested, resulting in an additional shortening of the quantum circuit. With these strategies, VQE calculations on a noiseless quantum simulator achieve energies within a few meVs of those obtained with the full UCCSD ansatz for $\mathrm{H}_4$ square, $\mathrm{H}_4$ chain and $\mathrm{H}_6$ hexagon molecules; while the number of controlled-NOT (CNOT) gates, a measure of the quantum-circuit depth, is reduced by a factor of as large as 35. Furthermore, we introduced an efficient"score"parameter to rank the excitation operators, so that the operators causing larger energy reduction can be applied first. Using $\mathrm{H}_4$ square and $\mathrm{H}_4$ chain as examples, We demonstrated on noisy quantum simulators that the first few variational operators can bring the energy within the chemical accuracy, while additional operators do not improve the energy since the accumulative noise outweighs the gain from the expansion of the variational ansatz.


I. INTRODUCTION
Quantum computers have been projected to be the ultimate solution to classically intractable problems owing to the exponential expansion of information that can be processed on quantum bits (qubits) compared with classical bits. However, to fully realize the advantage of quantum computing, quantum devices that integrate thousands or more qubits with sufficiently long coherent time have to be developed, which remains a significant challenge as of today. In the foreseeable future, one still has to work with so-called noisy intermediate-scale quantum devices (NISQ) [1], and practical quantum algorithms need to be made aware of this limitation. One group of such algorithms is the variational quantum eigensolver (VQE) [2], which uses a variational approach to optimize an objective function on a quantum/classical hybrid architecture. The preparation of parametrized quantum states and measurement of the expectation value of the objective function are performed on a quantum computer with relatively shallow circuits, while an optimization algorithm is implemented on a classical computer to find the optimal parameters.
Quantum chemistry has been one of the most active fields for quantum computing [3], realizing a proposal of solving quantum-chemical problems on a quantum architecture that Feynman made nearly 30 years ago [4]. * fzhang@ameslab.gov † ykent@iastate.edu There has been significant development of using VQE to solve quantum-chemical problems in both theory [5][6][7][8][9][10][11][12][13][14][15] and experiments on real quantum devices [2,7,10,[16][17][18][19]. The most commonly used ansatz is derived from the unitary coupled cluster (UCC) method [11,20,21], which is an extension of the well-known coupled cluster theory for describing the correlation effects in quantum systems [22]. In most applications, only single and double excitations are included, resulting in UCCSD. With this truncation, UCCSD in general cannot reach the true ground state energy. Another necessary step for implementing UCCSD is Trotterization [23], that is, the expansion of e A+B as e A/n e B/n n , where e A/n and e B/n can be efficiently implemented on quantum computers using available one-and two-qubit gates [3]. This expansion is exact only in the limit of n → ∞ when operators A and B do not commute. Conventionally, only a single Trotter step (n = 1) was used to represent UCCSD; more trotter steps barely improve the ansatz while significantly elongate the quantum circults [7,17]. A variation of UCCSD was also introduced recently, which, by successively adding operators to the ansatz one at a time in an adaptive process, can reach the accuracy of the true ground state with relatively shallow circults [14].
Despite the truncation and Trotterization, the implementation of UCCSD VQE on real devices has still been limited to small molecules, including H 2 , HHe + , LiH and BeH 2 [2,7,10,[16][17][18]. We note that VQE has also been applied to effective interacting few-site models that emerge from infinite lattice systems within Gutzwiller embedding theory [19,24,25]. In this paper, we use the intrinsic symmetry to further simplify UCCSD to meet NISQ requirements. The implementation of symmetry in constructing variational ansätze is not a new concept. In fact, the particle-number symmetry and Z 2 symmetry have been fully considered in selecting the excitation operators in UCCSD [11,14]. The point group symmetry of a molecule can also prohibit certain spin-orbital excitations [17,26]. More specifically, the point-group symmetry can further divide the Hilbert space preserving quantum numbers associated with the particle-number and Z 2 symmetries into several subspaces. Here, we will introduce an efficient graph clustering technique to identify these subspaces in the qubit representation. This general scheme allows us to systematically identify the most relevant excitation operators for the purpose of constructing a trial state with an energy close to the ground state energy. Based on the Hamiltonian matrix, this method is numerically cheap and does not require a sophisticated group theoretical analysis of the problem. It can be further combined with other strategies to reduce the gate complexity of the resulting circuit. Below, we will establish an importance ordering among the different excitation operators and combine all subterms associated with a particular operator. The resulting circuits are significantly shortened, allowing us to efficiently reach chemical accuracy for molecules H 4 and H 6 . The calculations were performed using the toolkit QISKiT developed by IBM [27], with both noise-free statevector and noisy QASM simulators.

II. FORMALISM AND RESULTS
The second quantization is applied to construct the Hamiltonian of the molecules. Atomic orbitals in the minimal basis (STO-3g) [28] are used in the calculations. Relevant spin-orbitals for constructing the basis of the Hilbert space are determined according to the Hartree-Fock (HF) calculations. In the second-quantized formulation, the electronic Hamiltonian is expressed as where h pq and h pqrs are one-electron and two-electron integrals, respectively, and σ and λ denote spins. h pq and h pqrs are calculated with the PySCF package [29]. The creation and annihilation operators in Eq. 1 are defined on 2N spin-orbitals (N is the total number of electrons). In order to solve the Hamiltonian on a qubit-based quantum computer, it is necessary to transform the Fock state |f 2N , f 2N −1 , · · · , f 1 , where f i is the occupation number of the i th spin-orbital (0 or 1), to a qubit state |q 2M , q 2N −1 , · · · , q 1 with M ≤ N . Accordingly, the ferminoic operators in Eq. 1 are transformed to qubit operators that can be realized in quantum circuits based on Pauli gates. We use the parity encoding method [30], in which the p th qubit stores the parity of the total occupation number of the first p spin-orbitals: (mod 2). If the spin-orbitals in the Fock state are arranged in such a way that the first N spin-orbitals describe spin-up states and the last N spin-orbitals describe spin-down states, then q N and q 2N are equal to the number of spin-up electrons (mod 2) and the number of electrons (mod 2), respectively. For non-relativistic molecules, these two numbers will be conserved. Consequently, the two qubits q N and q 2N will only be acted on by identity or Pauli Z operators, which can then be replaced by the corresponding eigenvalues, resulting in a Hamiltonian that only acts on 2M = 2N − 2 qubits [3]. In other words, with the parity encoding method, one can effectively save two qubits in quantum computing.
VQE employs the Rayleigh-Ritz variational principle where E 0 is the ground state energy and θ are variational parameters. A quantum circuit with moderate depth is used to apply a unitary operator U ( θ) on the initial state |ψ 0 , which is chosen to be the HF state in our calculations, creating a parametrized state |ψ( θ) : is measured on a quantum computer, while θ are varied to minimize the Rayleigh-Ritz quotient (left hand side of Eq. 2) on a classical computer. Unitary coupled cluster (UCC) is a chemistry-inspired ansatz that has been widely used in solving quantum-chemical problems [22]. In UCC, U ( θ) can be written as U (θ) = e θ(T −T † ) , where T can be any Hermitian excitation operator. However, only single and double excitations are usually selected, and UCC in this form is called UCCSD: where the subscripts αβ and ij denote occupied and virtual spin-orbitals, respectively. Using a single Trotter step, the UCCSD ansatz can be expressed as (4) The occupied and virtual spin-orbitals in Eq. 4 are selected in such a way that the net magnetization of the molecule is conserved. The total number of excitation operators in UCCSD for a non-magnetic system with N electrons can be calculated as: M = 2(N/2) 2 + N/2(N/2 − 1) + (N/2) 4 . M increases rapidly with N , making it challenging to implement the UCCSD ansatz for even moderate values of N on NISQ. For instance, M = 26 when N = 4, which already requires over a thousand CNOT gates to prepare the ansatz state ψ( θ).
An n-qubit system can span a Hilbert space with a dimension of 2 n . By construction, all the excitation operators included in UCCSD only act on a subspace preserving the total number of electrons and the net magnetization. This subspace separates into smaller subspaces if there are additional symmetry elements from the structure of the molecule. The ground-state of the Hamiltonian lies in one of these subspaces. Once this subspace is identified, one can confine the variational search within this subspace. That is, starting from an initial state |ψ 0 in this subspace, one only needs to apply the excitation operators that keep the ansatz states |ψ( θ) in the same subspace. In this way, the number of excitation operators in the variational ansatz can be greatly reduced.
Let us analyze the H 2 dimer as an example to illustrate this idea [7]. With two-qubit reduction applied in the parity encoding scheme, the fermionic system can be mapped onto a 2-qubit system, spanning a 4-dimensional Hilbert space. The four basis vectors in this qubit space |00 , |01 , |10 and |11 (parity encoding) correspond to the four Fock states |f 2↓ , f 1↓ , f 2↑ , f 1↑ = |0110 , |0101 , |1010 and |1001 , respectively. All the four Fock states conserve the total number of electrons (2) and the net spin (0). At a H-H distance of 0.725Å, the Hamiltonian in Eq. 1 can be represented by the following 4 × 4 matrix (in units of eV): By inspection, it is easily seen that the 4-dimensional Hilbert space can be separated into two subspaces S 1 and S 2 , spanned by {|00 , |11 } and {|01 , |10 }, respectively. The HF state is represented by the qubit state |01 (or the Fock state |0101 ). Starting from |01 , one needs to select excitation operators that flip both qubits so that the ansatz states remain in the same subspace. There are two single excitation operators and one double excitation operator in the full UCCSD for Here, X and Y are Pauli matrices, and the subscripts specify which qubit the Pauli matrix acts on. The two single hopping terms Y 1 and Y 2 transfer |01 onto |00 and |11 , respectively, both of which are out of the subspace where the initial state |01 is located. Therefore, the only relevant operator is the double excitation There is another simplification that can be made.
I is the identity matrix, and ǫ ijk is the parity of the permutation (ijk). Any state in the subspace spanned by |01 and |10 is an eigenstate of Z 2 Z 1 with an eigenvalue of -1. Thus, Z 2 Z 1 can be replaced with -1 when acting on this subspace: can be combined into one term iX 2 Y 1 , further reducing the circuit length by half.
We performed VQE calculations on the ground-state energy of H 2 molecule, using the QASM quantum simulator as implemented in the quantum computing toolkit QISKiT [27]. To simulate real NISQ devices, a noise model is implemented by including depolarizing gate errors for all qubits participating in the gate. To investigate the effect of circuit simplification, we implemented two different variational ansätze e iθ3/2(X2Y1−Y2X1) e iθ2Y2 e iθ1Y1 and e iθX2Y1 , corresponding to the full UCCSD and its simplified form, respectively. Since UCCSD is exact for H 2 , both ansätze are expected to give the exact diagonalization (ED) results without noise. The fermionic Hamiltonian in Eq. 1 is mapped onto a sum of tensor products of Pauli matrices (Pauli strings). The expectation value of each Pauli string was measured separately by averaging over 1024 shots. The error associated with the imperfect averaging is also included in the QASM simulator. The variational parameters are updated classically using the simultaneous perturbation stochastic approximation (SPSA) algorithm with 200 maximal iterations. The QASM simulations were performed 10 times independently to obtain an estimation of the error bar. In Fig. 1 (a), we plot the dissociation curve for H 2 , calculated by QASM simulations and ED. In Fig. 1 (b), we show the error for the two different ansätze. The simplified ansatz containing a single term iX 2 Y 1 gives smaller error at all bond lengths. While noticeable fluctuations exist in different trials as can be seen from the error-bar size, the average values for the single-term ansatz give acceptable accuracy, with an averaged error of 11.7 meV over all all measured bond lengths. For the full UCCSD ansatz, the error bars are larger, and the averaged error is increased to 27.1 meV, clearly demonstrating that the shortened quantum circuit results in a significant noise reduction.
The first excited state of H 2 lies in the subspace spanned by |00 and |11 . Thus, one can also obtain the energy of the first excited state by implementing VQE in this subspace [12,13]. We performed QASM simulations to verify this. Here, only a single term iX 2 Y 1 is included in the variational ansatz, but we choose |00 as the initial state. In Fig. 2 (a), we plot the first excitation state energy as a function of r, obtained from QASM simulations and ED; while in (b), we show the simulation error as a function of r. Again, accurate results can be obtained, with the averaged error being 6.6 meV.
It is a non-trivial task to identify the Hilbert space separation, especially when the ferminoic Hamiltonian is transformed to a qubit representation. We use the following algorithm based on graph clustering. A graph is created by denoting each basis vector of the Hilbert space as a node and connecting any two nodes i and j with an edge if the qubit Hamiltonian matrix element H ij is larger than a cutoff value: |H ij | > ǫ c . Here, ǫ c is set to 10 −6 eV. The problem to be solved is to separate the total space of connected nodes into isolated clusters so that any two nodes within the same cluster can be linked with a continuous path (not necessarily directed connected with an edge), while such a path does not exist for any two nodes belonging to different clusters. We provide the algorithm for solving this clustering prob- lem in pseudocodes in the Appendix. Each cluster will then represent a separate subspace. To ensure the VQE calculation is confined in one of these subspaces, we use the following procedure to select the fermionic excitation operators. Assume the m qubit basis vectors of the subspace correspond to m Fock states: |ψ 0 , |ψ 1 , · · · , and |ψ m−1 . |ψ 0 is pre-selected as the initial state for VQE calculations. Every other state ψ i (0 < i < m) is related to |ψ 0 by an excitation operator T i : |ψ i = T i |ψ 0 . The operators will be selected from all U i = T i − T † i (so that e θU is unitary), provided they satisfy the following two conditions: • U i contains only single or double excitation.
• U i does not transfer any of the basis vectors |ψ 0 , . . . , |ψ m−1 out of the subspace.
Each fermionic operator U i can be mapped onto a sum of tensor products of Pauli matrices [11]. While we have shown that all the terms in the sum can be exactly grouped into a single term in the example of H 2 , it is not achievable mathematically in general cases. On the other hand, a full treatment of all the subterms requires a significantly elongated quantum circuit. For instance, a unitarized double excitation operator transforms into 8 subterms in qubit representation [11], which requires a quantum circuit 8 times as long as that for a single term. On a noisy device or simulator where the noise grows cu- mulatively with the circuit length, it is likely that the extra noise associated with the elongated circuit outweighs the extra accuracy it gains from the extended variational degrees of freedom. Thus, it is helpful to compare the performance of the variational ansatz containing all the subterms for each excitation operator with the one that only contains a single term.
We performed calculations on H 4 molecule in the square configuration for a wide range of nearestneighboring H-H distance r. The Hamiltonian was constructed based on spin-orbital basis obtained from restricted HF calculations. To ensure that the HF calculation converges to states that are smooth with continuously varying H-H distances, the converged one-particle density matrix at one r was used as the starting point for the next r. The 8 spin-orbitals for H 4 were mapped onto a 6-qubit system with parity encoding and subsequent twoqubit reduction. The relevant Hilbert space for H 4 with zero net magnetization has a dimension of 4 2 · 4 2 = 36. Using the graph clustering algorithm, this space can be further separated into 4 subspaces, with dimensions of 8, 8, 10 and 10, respectively. The same partitioning of the Hilbert space can also be obtained for the Hamiltonian constructed based on symmetrized superposition of natural atomic orbitals (i.e., without the self-consistent HF calculations), showing the partitioning comes from the intrinsic point-group symmetry of the molecule, which is D 4h in this case. Since our method depends solely on the Hamiltonian matrix, a detailed group theoretical study is not necessary to identify the partitioning of the Hilbert space. Fig. 3 shows the VQE calculations in which each excitation operator was represented by the first term based on the lexicographical order. For example, among the 8 Pauli strings resulting from the double excitation operator a † 8 a † 4 a 1 a 5 : i , only the first term iY 6 X 5 X 4 X 3 X 2 X 1 was included in the variational ansatz. VQE calculations were performed in all four subspaces separately, using the noiseless statevector simulator implemented in QISKiT [27]. The statevector simulator uses matrices, rather than Pauli gates, to represent the qubit operators. Thus, it does not involve any noises associated with gate infidelity or imperfect averaging. In each subspace, a reasonable choice of the initial state is the basis state with the lowest diagonal element of the Hamiltonian matrix. An optimized energy close to the ground-state energy calculated by ED can be obtained in the subspace containing |001011 for all H-H distances, clear demonstrating that the ansatz with single-term representation for excitation operators is sufficient to give satisfactory accuracy. We also verified that the final answer is not sensitive to which term was selected. In the following, only single-term operators will be considered unless otherwise noted.
A closer inspection of the ED solution can reveal that there is an energy level crossing between two states with different total angular momentum S = 1 and S = 0. In Fig. 4 (a), we plot the energy of the two states as a function of the H-H distance. The energy difference ∆E = E(S = 0) − E(S = 1) is shown in Fig. 4 (b). At bond lengths r < r c = 0.8Å, the S = 1 state has the lowest energy. The level crossing occurs at r c = 0.8Å, where the S = 0 state becomes more stable. In our VQE calculations, only E is minimized without conserving the S 2 quantum numbers. Therefore, due to near degeneracies of S = 0 and S = 1 at low energies, the energy optimized VQE wavefunction is spin contaminated. Since both these two eigenstates are located in the same subspace in which the VQE calculations were performed, and the third eigenstate in this subspace with S = 2 is distantly separated from the first two states in the vicinity of the equilibrium bond length, the optimized VQE state (blue triangles) is essentially a superposition of the S = 0 state and the S = 1 state in this range. In Fig. 4 (c), we show the overlap between the VQE state and the ED ground state defined as | Ψ VQE |Ψ ED | 2 . Initially, the overlap drops as the energy difference between the S = 0 and S = 1 states decreases. The smallest overlap occurs at r = 0.8Å, which reflects the level crossing when the S = 0 and S = 1 states become degenerate. Then, the overlap increases as the S = 0 and S = 1 states are separated again, and reaches a maximum at r = 1.4Å, which also coincides with the largest energy difference between the S = 0 and S = 1 states. The separation between the two eigenstates shrinks when r further increases. For large r approaching the atomic limit, the energy of the S = 2 state decreases and eventually becomes degenerate with the first two eigenstates. Consequently, the VQE state becomes a superposition of the three states when approaching this limit: the overlap of the VQE optimized state and the S = 0 eigenstate decreases to as low as 0.34 at r = 3.0Å. Since our variational ansatz is designed to minimize the energy, it is acceptable that the optimized VQE state is spin contaminated. However, if one wants to preserve the spin quantum numbers, further constraints can be applied in the selection of excitation operators (T ) to enforce T, S 2 = 0 [31].
It is also worth noting that the VQE leading to the lowest-energy solution was performed in a subspace orthogonal to the restricted HF state |001001 , using |001011 as the initial state. In fact, the qubit state |001011 can be transformed back to the Fock state |00110101 , in which the spin-up electrons (the right half) and spin-down electrons (the left half) do not occupy the same spatial orbitals; while in the restricted HF state |00110011 (|001001 in qubit representation), each molecular orbital is doubly occupied by a spin-up and a spin-down electrons. Without partitioning the Hilbert space, the natural choice is to run VQE with the full set of UCCSD operators using the restricted HF state as the initial state. We showed the results of such calculations on the statevector simulator in Fig. 4 (a) as the black circles, where one can see that VQE failed to reach a satisfying estimation of the ground-state energy. This shows that the restricted HF state is not always a good choice as the initial state in VQE calculations, and the partitioning of Hilbert space as performed in the current work can help identify a suitable alternative choice.
The order in which the operators are included in the ansatz also matters, since it is preferable to add operators that create relatively high energy reductions first when the noise is still under control. A related idea of building an effective variational ansatz was illustrated in the ADAPT-VQE algorithm [14]. Here, we design the following process based on the second-order perturbation theory to efficiently rank the operators. For each U i connecting |ψ 0 to |ψ i , we denote the Hamiltonian matrix elements H 00 = ǫ 0 , H ii = ǫ i and H 0i = H i0 = ǫ 0i . Then, we define a "score" s i = min(|ǫ 0i |, ǫ 2 0i /|ǫ 0 − ǫ i |). The operators will be ranked in the descending order of s i . In Fig. 5, we show the VQE results by adding excitation operators to the variational ansatz, one at a time according to a lexicographical order as well as decreasing order of s i . The H-H distance was kept at 1.2Å, which is the equilibrium distance as seen in Fig. 3. |001011 was set as the initial state. It can be clearly seen that by adding the operators with large s i first, one can achieve a faster drop of the variational energy at the beginning. Since the noise-free statevector simulator was used, the two different sequences eventually lead to the same energy when all operators are included. The final energy is 0.014 eV above the ED result, which is well within the chemical accuracy of 1 kcal per mole, or 0.043 eV per molecule. Also shown in Fig. 5 are the VQE results by including all subterms for each excitation operator, following the decreasing order of s i . When only one operator was in- cluded, the results with a single term or all subterms are exactly the same. On the other hand, when additional operators were added, the VQE energy with all subterms included is slightly lower than that with only a single term. This is expected since more variational degrees of freedom are allowed with the additional Pauli terms. However, the improvement is limited. With all 6 variational operators added, the energy difference between all terms and a single term is only 4.7 meV. Calculations on noisy QASM simulators are also shown in green triangles in Fig. 5. The error bar was determined based on 10 independent calculations. When 4 or less number of excitation operators were included in the variational ansatz, results consistent with statevector simulators can be achieved on the noisy simulator. However, when the number of operators exceeds 4, the error bar significantly increases, and no further improvement of the energy can be obtained.
Two other geometric configurations were studied: an equidistant linear H 4 chain and hexagonal H 6 . For the H 4 chain, partitioning of the Hilbert space results in two relevant subspaces with the dimension of 16 and 20, respectively. The ground state can be approached by performing VQE in the 20-dimensional subspace using the HF state as the starting point. A total number of 14 excitation operators were selected and ranked according to the score parameter introduced in the above. As shown in Fig. 6 (a), on a statevector simulator, 10 excitation operators added in the decreasing order of the score can reach the ED ground-state energy within 1 meV. The remaining 4 excitation operators essentially have no effect on the optimized VQE energy. Therefore, only the first 10 operators were considered in QASM simulations. Again, 10 independent runs were performed for statistical analysis. Similar to the case of H 4 -square, only the first 7 excitation operators resulted in reduction of the VQE energy, while further addition of excitation operators did not help because the accumulative noise became too big. Nevertheless, QASM simulations can still reach the lowest energy only 10 meV above the ED result, which is well within the chemical accuracy.
The hexagonal H 6 can be mapped onto a 10-qubit system with parity encoding and subsequent application of 2-qubit reduction enabled by the Z 2 symmetry. The relevant subspace in the qubit system has a dimension of 6 3 · 6 3 = 400 for 6 total electrons and zero net magnetization. This subspace which can be further partitioned into 4 smaller subspaces with dimensions of 96, 96, 104 and 104, respectively. Similar to the H 4 -chain case, we identified that the subspace containing the restricted HF state is where the ground state is located. 41 operators up to double excitation can be selected in this subspace. In Fig. 6 (b), we plot the change of the VQE energy with successive addition of these excitation operators following the decreasing order of the score. One can see that in this case, the score parameter does not fully describe the "importance" of the operator, since operators 13-16 generate larger energy reductions than the operators immediately preceding them. This is not surprising because underlying perturbation theory can fail when the amplitudes of excitation operators become large. Nevertheless, our scheme still identifies most operators that cause significant energy drop with little extra computational cost. In fact, nearly half the the operators that essentially have no effects on energy reduction were found and put to the end of the list. Alternatively, the ADAPT-VQE algorithm [14] uses iteratively evaluated gradients of the energy with respect to the operator amplitudes to rank the operators. Since a separate optimization is required in each iteration, this treatment, while being more accurate, is considerably more expensive computationally.
The CNOT gate is an essential component in gatebased quantum computers; and the number of CNOT gates is a good indicator of the depth of the quantum circuit. In Fig. 7, we plot the number of CNOT gates in VQE with three different variational ansätze: single-term representation in subspace, all-subterm representation in subspace, and the full USSCD. For the hexagonal H 6 , the full UCCSD circuit contains 9,600 CNOT gates, which is out of the working range of the current NISQ. On the other hand, implementing the ansatz with the single-term representation for all operators in the relevant subspace only requires 260 CNOT gates, a reduction by a factor of 35. Since no quantum speedup can be gained on a classical simulator, it is still time-consuming to simulate H 6 on the noisy QASM simulator, even with the signifi- cant simplification of the circuit: it takes a few hours to prepare the parametrized variational state and measure its expectation value with 1024 shots. For this reason, we chose not to include the QASM calculations on H 6 without hurting the main conclusions.

III. CONCLUSION
We introduce a graph clustering algorithm to partition the Hilbert space into subspaces that is made possible by the intrinsic point group symmetry of the molecular systems. This step significantly reduces the number of variational operators since VQE ansätze can be confined to act within a particular subspace. Besides, it helps to obtain excitation energies, as shown for the case of H 2 ), or identify the correct initial state, as demonstrated for H 4 . Each excitation operator in UCCSD can be transformed into multiple Pauli terms, requiring a lengthy circuit to represent. We demonstrate with various examples that a single-term representation of excitation operators can reach required accuracy, while dramatically shortening the quantum circuit. VQE calculations on noiseless statevector quantum simulators achieve energies within a few meVs of those obtained with the full USSCD ansatz for H 4 square, H 4 chain and H 6 hexagon molecules. A "score" parameter was introduced at little extra computational cost, which allows us to rank the excitation operators so that the operators causing larger energy reduction can be applied first. Using H 4 -square and H 4chain as examples, we demonstrate on noisy quantum simulators that only the first few variational operators identified with this strategy are effective in reducing the energy to within the chemical accuracy.