# Digital quantum simulation of strong correlation effects with iterative
quantum phase estimation over the variational quantum eigensolver algorithm:
$\mathrm{H_{4}}$ on a circle as a case study

Dipanjali Halder1,∗, V. S. Prasannaa2,∗, Valay Agarawal1, Rahul Maitra1,† 1
Department of Chemistry,
Indian Institute of Technology Bombay, Powai, Mumbai 400076, India
2 Centre for Quantum Engineering, Research and Education, TCG CREST, Salt
Lake, Kolkata 700091, India
∗ These authors have contributed equally to this work
†<EMAIL_ADDRESS>
###### Abstract

The iterative quantum phase estimation algorithm, applied to calculating the
ground state energies of quantum chemical systems, is theoretically appealing
in its wide scope of being able to handle both weakly and strongly correlated
regimes. However, the goodness of the initial state that is sent as an input
to the algorithm could strongly decide the quality of the results obtained. In
this work, we generate the initial state by using the classical-quantum hybrid
variational quantum eigensolver algorithm with unitary coupled cluster ansätz.
We apply the procedure to obtain the ground state energies of the $H_{4}$
molecule on a circle, as the system exhibits an interplay of dynamic as well
as static correlation effects at different geometries. Furthermore, we argue
on the importance of static correlation in construction of the reference
determinant, and propose a minimally parametrized unitary coupled cluster
ansätz, which drastically reduces number of variational parameters while
incorporating the static correlation effects in the wavefunction. We
demonstrate that a carefully and appropriately prepared initial state can
greatly reduce the effects of noise due to sampling in the estimation of the
desired eigenphase.

##  I Introduction

Quantum computers offer great promise in efficiently handling certain problems
that are otherwise intractable on classical devices shor ; dj ; hhl . Over the
last two decades, the emerging field of quantum science and technologies has
witnessed remarkable progress with applications spread across physics,
chemistry, mathematics, biology, and finance math ; bio ; fin . In particular,
novel approaches to calculate ground state energies of quantum chemical
systems, and which promise an exponential speed-up over classical devices,
have made these class of problems a killer application of quantum computing
(for example, see Ref. ag ; ag1 ; ag2 ; Whitfield ).

In a seminal work, Aspuru Guzik et. al. carried out a digital quantum
simulation of two molecular systems, using the full configuration interaction
(FCI), and computed their ground state energies ag . The approach was based on
the famous quantum phase estimation (QPE) algorithm of Abrams and Lloyd Abrams
and Lloyd 1997 ; Abrams and Llyod 1999 . The QPE algorithm will lead to a
significant reduction in time complexity as compared to classical devices in
the future when many more qubits are available.

Based on a complementary paradigm, the variational quantum eigensolver (VQE)
algorithm is a hybrid quantum-classical algorithm, which aims to find the
lowest eigenstate of the interacting many-body Hamiltonian by varying the
energy expectation value function with respect to the parameters in the
wavefunction ansätz Peruzzo 2014 ; McClean 2016 . The method relies upon the
Rayleigh-Ritz variational principle Griffiths  and therefore provides an upper
bound to the exact eigenenergy of the many-body Hamiltonian. It is worth
adding at this point that the choice of a good ansätz can, in general, give
results with excellent precision. One such widely used variational form is the
chemically motivated unitary coupled cluster (UCC) ansätz. The coupled cluster
approach, which is widely regarded as the gold standard of electronic
structure calculations, is known to efficiently capture dynamical correlation
in a system. However, in the regions of molecular strong correlation, which is
attributed to the quasi-degeneracy of multiple constituent determinants, the
traditional VQE algorithm may not fare well, and as a result, suffer from
inconsistent description over the molecular potential energy surface (PES)
Sugi1 . On the other hand, the QPE algorithm does not possess the limitation
if appropriate care taken to prepare the initial state, as we shall discuss
next.

The QPE algorithm qpe_terhal  can be thought of as a projection of an input
state (called the reference state or the initial state interchangeably) onto
an eigenstate, whose eigenvalue we then compute with the algorithm. This
necessitates the initially prepared reference state having at least a non-
negligible overlap with the exact eigenstate of the Hamiltonian, and hence
preparation of a good initial state is of utmost importance. There have been a
number of works in literature that have addressed the issue of constructing an
initial state that has substantial overlap with the eigenstate of the system
Hamiltonian. In the molecular weak correlation regime, the simplest choice,
namely the single determinant Hartree-Fock (HF) function, often serves this
purpose quite well. However, this is often not the case while dealing with
strongly correlated systems, and in such scenarios, using a combination of
multiple Stater determinants as the initial reference function becomes very
important. This line of thought has led various authors to construct the
initial state by various means, for example, by employing the classically
computed Complete Active Space Self-Consistent Field (CASSCF) method Pittner ;
Hoffmann  or using Adaptive Sampling CI refhttps://arxiv.org/abs/1809.05523
techniques. Both of these approaches have shown substantial squared overlap of
the initial reference state with the eigenstate of the Hamiltonian, and thus a
high probability of success for the QPE algorithm. While the CASSCF method is
extremely powerful and useful for strongly correlated molecular systems, the
associated computational cost increases steeply against the size of the active
space. Thus, it is quite important to develop a wavefunction preparation
technique that should lead to handling the strong and weak correlation effects
in a balanced manner, and with relatively low computational scaling. Sugisaki
et. al. (doi/10.1021/acscentsci.8b00788)  have computed the reference state by
exploiting the diradical characters of the molecules from spin-projected
Unrestricted HF wavefunctions. This method entirely bypasses any post-HF
calculations for the reference state preparation, and drastically improves the
success probability of QPE, particularly for molecules with intermediate
diradical characters. A complementary approach for initial state preparation
is often followed, where the HF Hamiltonian is slowly and systematically
changed over to the FCI Hamiltonian, such that the system remains in its
instantaneous ground state. This is known as the Adiabatic State Preparation
(ASP) in literature ag ; Sugi2 ; Sugi3 .

In this work, we employ an iterative version of QPE algorithm (IQPE) iqpe ;
iqpe1 , and opt for an alternative way for preparing its input reference state
by using the VQE algorithm with the UCC ansätz. We refer to this approach
hereafter as IQPE over VQE. A schematic of the approach is given in Fig. 1. We
pick IQPE over traditional QPE for our simulations, in view of many ancilla
qubits that are required for the latter, which in turn makes computations very
expensive. To the best of our knowledge, there is one other work in
literature, where the authors independently suggest QPE over VQE, but they do
not study the effect of various VQE parameters in the initial state
preparation Romero . In this work, we systematically tune the number of VQE
iterations, and study its connection to the probability of successfully
landing on the right state when fed as an input to the IQPE algorithm. We will
demonstrate that a multi-determinantal reference state meticulously prepared
using the VQE algorithm is likely to have a superposition of all the dominant
determinants of the FCI wavefunction, and thus we anticipate that one can
systematically tune it to have considerable overlap with the exact eigenstate
of the Hamiltonian, even in strong correlation regimes. We then demonstrate
the strength and applicability of this method by carefully selecting a
physical system that displays varied electronic complexity across its PES,
that is, the system should be dominated by strong correlation effects in some
geometries and weak in others to serve as a good testbed for the algorithm.
After generating the PES and studying the precision of our final results with
the number of IQPE iterations, we show that our approach provides a good and
systematically improvable probability of success to land on the correct state.
We would also introduce a minimally parametrized CC ansätz, where the multi-
determinantal reference function is generated by incorporating the static
correlation effects via a unitary ansätz with internal excitations. This
approach promises to be a good compromise between the accuracy in terms of the
squared overlap and the number of variational parameters, while constructing
the reference state. We finally extend our analyses by studying the scenario
with Qiskit's QASM backend, where sampling comes into play.

The paper is organised as follows: we discuss in Section II the VQE, QPE, and
the IQPE algorithms, followed by the hybrid IQPE over VQE approach. This is
followed by the details of our computations and the rationale behind our
choice of system. Finally, we present our obtained results and the analyses
that follow from our data in Section III, where we systematically expound
details such as the PES, effect of IQPE iterations in precision of our
results, the effect of VQE iterations on state preparation for IQPE, he
minimally parametrized ansätz, and finally analysis of our results by taking
into account sampling. We conclude in Section IV.

##  II Theory

###  II.1 The Variational Quantum Eigensolver algorithm

As mentioned earlier, the classical-quantum hybrid VQE algorithm relies upon
the variational principle , and hence the ground state energy, $E(\theta)$,
given by

$\displaystyle E(\theta)$ $\displaystyle=$
$\displaystyle\frac{\langle\Psi(\theta)|H|\Psi(\theta)\rangle}{\langle\Psi(\theta)|\Psi(\theta)\rangle}$
(1) $\displaystyle=$
$\displaystyle\langle\Psi_{HF}|U(\theta)^{\dagger}HU(\theta)|\Psi_{HF}\rangle,$
(2)

is, in principle, guaranteed to be an upper bound to the true ground state
energy. Here, the state, $|\Psi(\theta)\rangle$, is expressed in terms of a
unitary operator, $U(\theta)$, which is parametrized by a set of parameters
concisely denoted as $\theta$, and $U(\theta)$ acts upon a reference state,
$|\Psi_{HF}\rangle$, which is the HF state in this work. We choose $U(\theta)$
to be the UCC ansätz with double excitations (UCCD) unless mentioned
otherwise, and it is given by

$\displaystyle\ket{\Psi(\theta)}$ $\displaystyle=$ $\displaystyle
e^{{{T_{2}}({\theta})}-{T_{2}}^{\dagger}({\theta})}\ket{\psi_{HF}};$ (3)
$\displaystyle{T_{2}}({\theta})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum\limits_{i,j;a,b}\theta_{ij}^{ab}a_{b}^{\dagger}a_{a}^{\dagger}a_{j}a_{i}.$
(4)

Here, the subscripts $i,j$ refer to the occupied spin-orbitals in the HF
reference determinant, while $a,b$ are the unoccupied/virtual spin-orbitals.
The t-amplitudes, $\theta_{ij}^{ab}$s, play the role of the VQE parameters.
The choice to not consider single excitations has to do with the system that
we work with, where we expect the double excitations to be dominant, and we
shall explain and demonstrate this point in Section III. One may also choose
only a selection of double excitations involving the chemically active
orbitals. This important point will be elaborated towards the end of the
manuscript, in the context of the minimally parametrized ansätz. Note that in
Eq. (2), the many-body molecular electronic Hamiltonian, $H$, can also be
expressed in second quantized form as

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum\limits_{p,q}h_{pq}a_{p}^{\dagger}a_{q}+\frac{1}{2}\sum\limits_{p,q,r,s}h_{pqrs}a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}$
(5) $\displaystyle=$ $\displaystyle\sum_{\alpha}h_{\alpha}P_{\alpha},$ (6)

where $h_{pq}$ and $h_{pqrs}$ refer to one- and two- electron integrals,
respectively. The integrals are supplied as inputs to the VQE algorithm, using
a program run on a classical computer. In Eq. (6), the Hamiltonian is
compactly described using a single index, $\alpha$, and the $P_{\alpha}$s no
longer refer to the second quantized creation and annihilation operators, but
their mapped version into their qubit operator form, that is, a string of
tensor product of Pauli operators. This can be done by, for example, the
Jordan-Wigner transformation map . $T_{2}(\theta)$ can also be appropriately
mapped into its qubit operator form, and this allows $|\Psi(\theta)\rangle$ to
be recast in its circuit form Whitfield , with appropriately chosen initial
guess parameters, $\theta_{init}$. The expectation value is calculated by
either using matrix operations (statevector (SV) backend) or the circuit is
evaluated by measuring $H$ on $|\Psi(\theta_{init})\rangle$ (QASM backend) and
the corresponding energy extracted. This sub-part constitutes the quantum
module of the classical-quantum VQE approach. The energy is then minimized
with respect to the parameters by using a suitable optimization algorithm on a
classical computer, over a series of iterations, where in each VQE iteration,
the parameter $\theta$ incrementally changes, and is fed back into the circuit
to obtain an updated value of energy for that iteration, and the process is
repeated until a minimum is found.

The excitation operators that occur in the exponent of Eq. (3), when acting on
the reference HF determinant, capture the physical effects arising from the
resulting determinants that differ from their HF counterpart by double
excitations. With a sufficient number of VQE iterations for energy
minimization, the HF reference function thus evolves to a correlated ground
state wavefunction, which is expected to have a substantial overlap with the
exact ground eigenstate. While this expectation is indeed the case for
molecules in weakly correlated regimes, the same is not strictly true for
strong correlation, which we shall explain later with a concrete example,
$H_{4}$ on a circle. This leads VQE to provide a poor estimation of the energy
in the presence of molecular strong correlation, as we will demonstrate later.
However, such a correlated wavefunction resulting as a byproduct of the
iterative optimization of the energy may be used for QPE with a high
probability of its success, even in in the regions of molecular strong
correlation.

In the following subsection, we will briefly discuss the general principles of
the QPE algorithm starting from an arbitrary reference function, and highlight
the importance of reference state preparation in this regard. We will then
demonstrate how the evolved state resulting from VQE energy optimization can
suitably be fed as the reference state into the QPE algorithm.

###  II.2 Quantum phase estimation algorithm

The algorithm aims at estimating the energy eigenvalues of a physical
Hamiltonian. The objective is to estimate the phase, $\phi$, of an eigenvector
of a unitary operator, ${U}$, such that

$\displaystyle{U\ket{\Psi}=e^{2\pi{i}\phi}\ket{\Psi}}.$ (7)

The algorithm uses two registers. One of them contains a set of $k$ ancilla
qubits, prepared initially in a state, $\otimes_{i=1}^{k}\ket{0}$, while the
other one contains a set of qubits required to define the initial state of the
system. This arrangement enables phase kickback via a series of controlled
unitary gates. The inverse quantum Fourier transform is then employed to
change bases. This step is followed by extracting the phase, $\phi$. The
precision of the algorithm depends on $k$. However, in practice, the number of
ancilla qubits which can be implemented is very limited, in view of the
current Noise Intermediate Scale Quantum (NISQ) era that we live in, thus
restricting the applicability of the promising approach to small systems
and/or limiting the precision. An iterative version of QPE, namely the IQPE
algorithm, uses a single ancilla qubit to solve the eigenvalue problem.

###  II.3 The iterative quantum phase estimation algorithm

The IQPE algorithm carries out the conventional QPE with only a single ancilla
qubit. The precision of the obtained results no longer depends on the number
of ancilla qubits, but on the number of iterations. IQPE too requires two
registers. However, the first register contains only a single ancilla qubit,
initially in the state $\ket{0}$, and the second register contains the qubits
necessary for defining the initial state of the system, $\ket{\Psi_{init}}$.

The phase, $\phi$, which is to be estimated, can be written down as a series
of bits where each bit, $\phi_{i}\epsilon\\{0,1\\}$.

$\displaystyle\phi=\frac{\phi_{1}}{2}+\frac{\phi_{2}}{4}+\ldots+\frac{\phi_{m}}{2^{m}}=0.\phi_{1}\phi_{2}\ldots\phi_{m}.$
(8)

The algorithm starts by estimating the least significant bit, $\phi_{m}$, for
which a Hadamard gate is applied to the ancilla qubit, $\ket{q_{0}}$, so that
it is transformed to $\ket{+}$. This is followed by a controlled-$U^{2^{m-1}}$
operation between the ancilla qubit register ($q_{0}$) and the state register
($q_{1}$), such that the combined state of the ancilla and the state register
becomes

$\displaystyle\ket{q_{0}}\ket{\Psi_{init}}\longrightarrow\ket{0}\ket{\Psi_{init}}+e^{i2\pi
0.\phi_{m}}\ket{1}\ket{\Psi_{init}}.$ (9)

Hence, the phase that is kicked back into the ancilla is $0.\phi_{m}$. As
$\phi_{m}$ is either 0 or 1 , a measurement on ancilla will lead to $\ket{+}$
if $\phi_{m}=0$, and $\ket{-}$ if $\phi_{m}=1$. This would produce outcomes
$\ket{0}$ and $\ket{1}$, respectively, when measured in the x-basis. After the
measurement of the least significant bit, $\phi_{m-1}$ is measured, using the
value of $\phi_{m}$. For estimation of $\phi_{m-1}$, $2^{m-2}$ controlled-U
operations are applied between $q_{0}$ and $q_{1}$ such that

$\displaystyle\ket{q_{0}}\ket{\Psi_{init}}\longrightarrow\ket{0}\ket{\Psi_{init}}+e^{i2\pi
0.\phi_{m-1}}e^{i2\pi\frac{\phi_{m}}{4}}\ket{1}\ket{\Psi_{init}}.$

Thus, the phase that is kicked back is $0.\phi_{m-1}+\frac{\phi_{m}}{4}$. To
extract the bit $\phi_{m-1}$, a phase correction of $(-2\pi\phi_{m}/{4})$ is
applied with the help of rotation gates, which transform the ancilla into
$\ket{0}+e^{i2\pi 0.\phi_{m-1}}\ket{1}$. The ancilla is then measured in the
x-basis to extract $\phi_{m-1}$. Similarly, the more significant phase bits
are measured using the knowledge of the previously measured lesser significant
phase bits, and the process is iterated backwards till all the phase bits are
recovered.

The above algorithm works even if the initial state $\ket{\Psi_{init}}$ on the
$q_{1}$ register is not an eigenfunction of the unitary operator, U but a
linear superposition of various eigenfunctions Cruz ,

$\displaystyle\ket{\Psi_{init}}=\sum\limits_{n}c_{n}\ket{u_{n}}.$ (11)

In such a situation, the algorithm will lead us to one of the eigenfunctions,
with a probability $|{c_{n}}|^{2}$ and provide an estimate of the local phase
of the eigenfunction $\ket{u_{n}}.$ In our present work, we would explore two
different unitary ansätze to construct multi-determinantal reference state via
VQE, and we will show that a careful construction of the reference state is
essential for the success of the IQPE algorithm.

An important point to note is that the initial state does not influence the
precision of the phase, rather it only affects the probability with which the
phase of a particular eigenstate is measured. The precision of the algorithm
itself is solely determined by the number of iterations in IQPE.

##  III Results and Discussions

|
---|---
(a) | (b)
Figure 1: An overview of the algorithm employed in this work. The choice of
system is $H_{4}$ on a circle is shown in subfigure (a). Our scheme to obtain
the ground state energy involves carrying out a variational quantum
eigensolver (VQE) computation, followed by iterative quantum phase estimation
(IQPE), as shown in subfigure (b). In the former, the unitary coupled cluster
variational form with double excitations, denoted by $U_{UCCD}(\theta)$, is
used. $\theta$ refers to the variational parameters, with the initial guess
parameters given by $\theta_{init}$. $U_{UCCD}(\theta)$ acts on the HF initial
state. $h_{\alpha}$ refers to the one- and two- electron integrals. The
parameters are updated in each evaluation, and the multi-determinantal state
thus prepared after optimization is fed into the IQPE algorithm. The figure
shows a schematic of the algorithm for the $n^{th}$ iteration.

###  III.1 Methodology

We obtained the one- and two- electron integrals from PySCF pyscf , and
carried out the VQE as well as the IQPE computations on Qiskit 0.26 qiskit .
We chose 1.738 Å for the radius of the $H_{4}$ circle as shown in Fig. one(a),
the same as that from Ref. h4 . We used the GAMESS gamess  package
(specifically, the determinant-based ALDET approach aldet1 ; aldet2 ; aldet3 )
for the relevant FCI coefficients to obtain success probabilities.

For our computations, we chose the contracted STO-3G basis sto  that contains
eight spinorbitals in total, COBYLA optimizer cobyla , Jordan-Wigner mapping
scheme map , and the UCCD variational form ucc . It is worth adding at this
point that UCCD with an untruncated double excitation space for $H_{4}$ on a
circle in the STO-3G basis, leads to an 18-parameter VQE computation. We
considered two options for $\theta_{init}$: zero initial guess, and the
parameters built from many-body perturbation theory taken to second order in
energy (MP2). We also worked with two choices for the backend simulator: the
statevector and the QASM backend. We set a Trotter step of one for all our
calculations. Also, we have employed the direct spinorbital to qubit mapping:
that is each spinorbital is represented as one qubit, independent of the
occupancy of the spinorbitals.

###  III.2 Application to $H_{4}$ on a circle

The case of four $H$ atoms on a circle is a standard and challenging test bed
to verify the efficacy of newly developed electronic structure algorithms. The
system is characterized by the angle, $\beta$, as shown in Fig. 1(a). For
values of $\beta$ away from 90°, for example, 80°, the system behaves as two
non-interacting $H_{2}$ molecules, and is hence a strongly single reference
case, with the correlated ground state wavefunction predominantly
characterized by the HF state. This is clearly seen in the following FCI
ground state generated via a determinant-based CI approach:

$\displaystyle\ket{\Psi_{FCI}^{80^{\degree}}}$ $\displaystyle=$ $\displaystyle
0.642\ket{00110011}+0.350\ket{11001100}$ (12) $\displaystyle-$ $\displaystyle
0.317\ket{01010101}+0.251\ket{01101001}$ $\displaystyle+$ $\displaystyle
0.251\ket{10010110}-0.241\ket{10101010}$ $\displaystyle+$ $\displaystyle
0.223\ket{01011010}+0.223\ket{10100101}$ $\displaystyle-$ $\displaystyle
0.200\ket{10011001}-0.194\ket{01100110}.$

In the above equation (and the subsequent Eqs. (13), (LABEL:vqe_80_500) and
(LABEL:vqe_89p8_3000)), only those determinants are included which have
absolute coefficients greater than 0.05. Also, throughout this manuscript,
each determinant has been written in its occupancy number representation, and
the indices are arranged in the block spin format, where all the $\alpha$
spins are first counted, followed by the $\beta$ spins. From Eq. (12), we
immediately see that the HF configuration dominates, followed by two doubly
excited determinants.

Beyond 85°, the system starts showing the signature of molecular strong
correlation, with one or more determinants becoming quasi-degenerate to the HF
determinant. Very close to 90°, the HF determinant and a doubly excited
determinant become equally dominant, as shown in the FCI wavefunction below:

$\displaystyle\ket{\Psi_{FCI}^{89.8^{\degree}}}$ $\displaystyle=$
$\displaystyle 0.519\ket{00110011}-0.509\ket{01010101}$ (13) $\displaystyle+$
$\displaystyle 0.345\ket{11001100}-0.340\ket{10101010}$ $\displaystyle+$
$\displaystyle 0.280\ket{01101001}+0.280\ket{10010110}$ $\displaystyle+$
$\displaystyle 0.143\ket{10100101}+0.143\ket{01011010}$ $\displaystyle-$
$\displaystyle 0.137\ket{11000011}-0.137\ket{00111100}.$

From the above equation, we see that at $89.8\degree$, the HF state and the
doubly excited determinants are equally dominant in the description of the
correlated exact ground state with $C_{HF}^{89.8\degree,GS}=0.519$ and
$C_{01010101}^{89.8\degree,GS}=0.509$, and hence one needs to treat all the
relevant dominant determinants in the same footing in order to have a balanced
description of dynamical and static correlation effects. Furthermore, there
are some excited states (denoted as $ES_{n}$) with same symmetry as the ground
state (denoted as GS), where the coefficient of the HF determinant is
sufficiently large in the corresponding FCI wavefunction, along with other
doubly excited determinants. For example, for the lowest excited state
$ES_{1}$, the coefficient of the HF determinant is -0.297 in the FCI
wavefunction, whereas there exist two other excited states with the same
symmetry (denoted as $ES_{2}$ and $ES_{3}$) where the HF determinant
contributions are 0.543 and 0.431, respectively. Thus, starting from the HF
state, the IQPE algorithm is probabilistically more likely to land up in the
excited roots, as we will show in the subsequent section. This is quite unlike
the scenario when $\beta=80\degree$. In our approach to follow, we would
employ iterative optimization via the VQE-UCCD ansätz to construct the
reference determinant, and we would demonstrate that the overlap of the
resultant wavefunction with the FCI wavefunction described in Eq. (12) and Eq.
(13) can systematically be improved, which leads to increase the probability
of success for the subsequent IQPE. Note that the ground state of $H_{4}$
system is of $A_{g}$ symmetry. It is clear that there is no single excitation
operator belonging to the $A_{g}$ symmetry, and therefore, the inclusion of
single excitation operators has no effect on the quality of the UCC wave
function.

The multi-determinantal state thus prepared using the VQE algorithm is
subsequently supplied as an input to the IQPE algorithm for phase estimation
as shown in Fig. 1(b). One may note that with a larger number of VQE
iterations, the initial reference state for IQPE gets a more substantial
overlap with the exact target function. Thus, with more VQE iterations, the
success probability of the IQPE algorithm to extract the correct phase is
expected to increase. For a fixed number of IQPE iterations, the precision of
an extracted phase remains unchanged.

###  III.3 The potential energy surface and the effect of the number of IQPE
iterations

As a first step, we determine the PES of $H_{4}$ on a circle. In particular,
we are interested in investigating the performance of our approach in both
weak and strong correlation regimes. For this purpose, we plot, as shown in
Fig. 2, the ground state energies against the angle, $\beta$, which was
defined in Fig. 1(a). We vary $\beta$ from 80$\degree$ to 100$\degree$. Note
that the plot is expected to be symmetric on both sides of $\beta=90\degree$;
however, due to inaccessibility of the electronic configurations at
$\beta=90\degree$, we stop at a value that is reasonably close to it h4 . The
plot also shows the FCI energies in the chosen range of angles, as it acts as
the exact result obtainable within a single particle basis, and hence a
benchmark for comparison. We will elucidate on two aspects: precision (with
respect to FCI) and non-parallelity error (NPE). It is evident from the
results shown in the Fig. 2 that VQE with UCCD ansätz performs very well and
agrees to $\sim$ 0.1 milliHartree ($mE_{h}$) with respect to FCI, between
$\beta=$ 80$\degree$ and $\beta=$ 85$\degree$, but starts deviating beyond
that point. This behaviour shows that in geometries where strong correlation
effects begin to become important, VQE no longer can make predictions of
ground state energies with the kind of precision that it demonstrates for
weakly correlated domains. In line with this observation, NPE provides a
proper quantifier to assess the efficacy of a theory to treat the entire PES
with same level of precision. NPE, across a PES, is defined as the difference
between the maximum and minimum deviation of the results from a theory with
respect to FCI Sherrill1 . That is, it quantifies the spread of the
(non-)parallelity error across the PES. We find from our data that the NPE is
more than 1 $mE_{h}$ for VQE when UCCD ansätz is employed, which is not within
acceptable bounds, given that the chemical accuracy itself is $\sim$ 1
$mE_{h}$. To check if IQPE over VQE circumvents the large NPE obtained from
VQE, we employed the IQPE algorithm with 16 iterations, with the initial
reference state prepared using VQE in the UCCD ansätz with 1000 iterations,
and using MP2 initial guess parameters. Our calculation show that the NPE was
only about 0.2 $mE_{h}$ for IQPE when started with the reference function
generated by VQE with UCCD ansätz. This is not surprising as a 16-bit IQPE is
supposed to provide energy, which is precise up to $\sim 10^{-4}E_{h}$, and
thus NPE is, in principle, controlled only by the number of iterations in
IQPE. However, with VQE-UCCD, the quality of the initial guess wave function
greatly affects the probability to obtain the correct ground state in the
subsequent IQPE stage. We will elucidate this point in the next subsection.

We should also note that IQPE over a HF reference determinant is in most
occasions unable to extract the correct phase that corresponds to the ground
state eigenfunction for all angles beyond 80$\degree$ (vide infra). Since the
HF determinant has significant contribution to the excited state wavefunctions
as compared to the ground state wavefunction, as described in the previous
section and also shown in Fig. 4, starting from the HF reference determinant
in most occasions leads to incorrectly landing on an excited state
eigenfunction. It is also worth adding that such a quantitatively accurate PES
obtained from IQPE over VQE with high probability of success can later be
extended to evaluating properties such as vibrational frequencies,
dissociation energies, etc. This attests the strength of IQPE over VQE,
provided the initial state from VQE is chosen well.

Figure 2: Potential Energy Surface of $H_{4}$ on a circle, versus the angle
$\beta$. The energies are given in units of Hartree ($E_{h}$). The blue curve
gives the FCI values over the chosen range of angles, and the converged VQE
results are provided in red. The results with 16 IQPE iterations over a state
prepared by 1000 VQE iterations (via COBYLA optimizer and MP2 initial guess)
are given in green color. The figure clearly shows that the non-parallelity
error is very small in the IQPE over VQE method, whereas the VQE approach
shows a marked deviation in the strong correlation regimes.  Figure 3:
Difference between FCI and the IQPE energies (in $E_{h}$) versus the number of
IQPE iterations, for three geometries: (a) $\beta=80\degree$, (b)
$\beta=85\degree$, and (c) $\beta=89.8\degree$. The initial reference function
was optimised via VQE algorithm with 1000 COBYLA iterations, starting from MP2
guess amplitudes.
---
(a)
(b)
(c)
Figure 4: Plots illustrating our results for three chosen angles, $\beta=$
80°(subfigure (a)), $\beta=$ 85°(subfigure (b)), and $\beta=$ 89.8°(subfigure
(c)), reflecting the increasing importance of strong correlation effects. Each
of the sub-figures contain the overlaps of both the ground and excited states
with the FCI state, defined by $\langle\Psi\arrowvert\Psi_{FCI}\rangle$,
plotted against the number of VQE iterations for a given angle. The blue and
black curves denote the ground (GS) and an excited state (ES1) overlaps
respectively, with zero initial guess for the VQE parameters. The cyan and
yellow curves present the same quantities, but for an MP2 initial guess for
the VQE parameters. Additionally, we also provide the overlaps for two other
excited states (ES2 and ES3), starting with MP2 guess, and they are shown as
red and green curves, respectively.

To understand the number of IQPE iterations at which the calculated energy
converges for $H_{4}$ on a circle, we have plotted in Fig. 3 the change in the
difference between the FCI energy and observed energy with number of IQPE
iterations, for $\beta=80\degree,\ 85\degree,\ \mathrm{and}\ 89.8\degree$. We
have used MP2 initial guess for VQE parameters, and set maximum number of VQE
iterations to 1000 for the state preparation. It is worth adding at this point
that we have not gone beyond 16 IQPE iterations in view of the increasingly
steep computational requirements. The figure shows that one needs to go at
least upto 13 IQPE iterations to reach sub-$mE_{h}$ precision, independent of
whether we are in weak or strong correlation regime. This is not surprising as
a 13-bit IQPE calculation leads to precision of $\pm(2\pi\times
2^{-13})E_{h}\sim 10^{-4}E_{h}$; it is simply the probability of correctly
extracting the phase to a pre-defined precision, while noting that the value
on which we land can greatly be enhanced in a controlled manner via tuning
various VQE parameters. The importance of an appropriate state preparation is
far more obvious for strong correlation regimes than single reference regimes,
which we would demonstrate later.

###  III.4 The effect of the number of VQE iterations: state preparation for
IQPE

The success of the IQPE algorithm depends on the closeness of the initial
reference state to the exact eigenstate of the many-body Hamiltonian. Although
the final precision of the estimated energy by the IQPE algorithm is
independent of the initial reference state, an appropriately prepared
reference state largely enhances the probability of extracting the correct
phase corresponding to the target state. In our approach, starting from the HF
determinant, we have constructed a multi-determinantal reference state by
including dynamical correlation through iterative optimization via the VQE
algorithm. Understandably, a sufficiently large number of VQE iterations would
ensure that the resulting multi-determinantal state has substantial overlap
with the exact ground eigenstate. Furthermore, it is observed that if the
initial state is poorly prepared (e.g. a HF state, or a multi-determinantal
state generated with a very low number of VQE iterations with zero or random
initial guess amplitudes), there is a non-negligible probability that IQPE
extracts a phase which corresponds to an excited state. We find that this is
particularly the situation for molecules in strongly correlated regime since
the UCCD ansätz often fails to account for the short range static correlation.
To this end, we have presented the squared overlap of the prepared state, with
the FCI ground (denoted as GS) as well as with lowest excited (denoted by ESn)
states, as functions of the number of VQE iterations for three different
angles arranged in increasing order of strong correlation, as shown in Fig. 4.
In the weak correlation regime, for which we have chosen 80$\degree$ as the
representative case, the HF determinant has a squared overlap with the FCI
ground state of nearly
$|\langle\Psi_{HF}|\Psi_{FCI}^{80\degree,GS}\rangle|^{2}=0.41$, while its
squared overlap with the lowest excited state,
$|\langle\Psi_{HF}|\Psi_{FCI}^{80\degree,ES_{1}}\rangle|^{2}=0.004$. We note
that we have introduced an additional index on the superscript wherever
necessary, to distinguish between GS and ES FCI states. With only a few VQE
iterations, the overlap between the reference state and the target FCI ground
state can significantly be improved: In fact, with only about 500 iterations
and a single VQE calculation starting from MP2 amplitudes, the correlated
state takes the form

$\displaystyle\ket{\Psi_{VQE(500)}^{80^{\degree}}}$ $\displaystyle=$
$\displaystyle 0.639\ket{00110011}+0.356\ket{11001100}$ $\displaystyle-$
$\displaystyle 0.322\ket{01010101}+0.211\ket{01101001}$ $\displaystyle+$
$\displaystyle 0.211\ket{10010110}-0.244\ket{10101010}$ $\displaystyle+$
$\displaystyle 0.255\ket{01011010}+0.254\ket{10100101}$ $\displaystyle-$
$\displaystyle 0.188\ket{10011001}-0.184\ket{01100110}$ $\displaystyle-$
$\displaystyle 0.073\ket{00111100}-0.074\ket{11000011},$

having a squared overlap of nearly one with the FCI ground state (shown in Eq.
(12)), that is,
$|\langle\Psi_{VQE(500)}^{80\degree}|\Psi_{FCI}^{80\degree,GS}\rangle|^{2}\rightarrow
1.0$, while
$|\langle\Psi_{VQE(500)}^{80\degree}|\psi_{FCI}^{80\degree,ES_{n}}\rangle|^{2}\rightarrow
0$. This ensures guaranteed estimation of the correct phase corresponding to
the ground state. On the other hand, in the strongly correlated regime
(89.8$\degree$), the HF determinant has significant squared overlap with both
the ground and a couple of excited states of the same symmetry, and hence the
phase estimation with HF state as the initial reference often incorrectly
leads to one of the excited states.

This issue can greatly be alleviated by optimising the reference determinant
via the inclusion of correlation effects through VQE; the overlap between the
VQE state and the FCI ground state steadily increases with the number of VQE
iterations and it saturates around 1000 VQE iterations when we begin with
either the MP2 or the zero initial guess amplitudes. Further optimization does
not lead to any improvement of the VQE energy and the prepared state. As an
example,

$\displaystyle\ket{\Psi_{VQE(3000)}^{89.8^{\degree}}}$ $\displaystyle=$
$\displaystyle 0.565\ket{00110011}-0.407\ket{01010101}$ $\displaystyle+$
$\displaystyle 0.378\ket{11001100}-0.270\ket{10101010}$ $\displaystyle+$
$\displaystyle 0.262\ket{01101001}+0.262\ket{10010110}$ $\displaystyle+$
$\displaystyle 0.193\ket{10100101}+0.193\ket{01011010}$ $\displaystyle-$
$\displaystyle 0.071\ket{11000011}-0.071\ket{00111100}$ $\displaystyle-$
$\displaystyle 0.124\ket{10011001}-0.120\ket{01100110}.$

This resulting multi-determinantal state, generated with 3000 iterations by a
single VQE calculation starting from the MP2 amplitudes, has about 95% overlap
with the FCI ground state:
$|\langle\Psi_{VQE(3000)}^{89.8\degree}|\Psi_{FCI}^{89.8\degree,GS}\rangle|^{2}=0.9$.
On the other hand, its overlap with the lowest excited state reduces
significantly to
$|\langle\Psi_{VQE(3000)}^{89.8\degree}|\psi_{FCI}^{89.8\degree,ES_{1}}\rangle|^{2}=0.057$.
Starting from zero guess amplitudes, we find that
$|\langle\Psi_{VQE(3000)}^{89.8\degree}|\Psi_{FCI}^{89.8\degree,GS}\rangle|^{2}=0.88$,
and
$|\langle\Psi_{VQE(3000)}^{89.8\degree}|\psi_{FCI}^{89.8\degree,ES_{1}}\rangle|^{2}=0.096$.
One may also note that the overlap of the VQE generated state with $ES_{2}$
and $ES_{3}$ sharply decays to nearly zero with very few VQE iterations, as
shown in Fig. 4(c). We argue that the construction of the multi-determinantal
reference state via VQE iterations leads to significantly enhanced overlap
with the target ground state, ensuring a quicker and secured evolution via
IQPE.

Figure 5: Figure illustrating the dependence of the quality of our results on
the VQE iterations (a measure of the goodness of the initial state for IQPE
procedure) as well as on the number of shots for 89.8°. The energies are
appropriately color coded, as shown on the right panel. The exact FCI energy
is also marked appropriately in the panel. Each circle denotes one repetition
of the experiment for a given number of shots and a given initial state
preparation. The number of same colored circles within a block, therefore,
reflects the frequency of occurrence of a given value of energy. Note the
appearance of two different colours (in subfigure (b) with 25, 50 and 100
shots, corresponding to two different closely spaced energy values obtained by
IQPE when the initial reference state is constructed via full double
excitation inclusive VQE. For the extreme case of 10000 shots, IQPE predicts a
single value for the energy, irrespective of the truncation scheme in UCCD
ansätz for the VQE state preparation.

We briefly comment on the error due to Trotter step, and found it to be well
within 0.1 $mE_{h}$. We also checked the error from the choice of mapping
(Jordan-Wigner, Parity, and Bravyi-Kitaev schemes), and found that too to be
well within 0.1 $mE_{h}$. This is not too surprising, as the transformed qubit
Hamiltonians from each of these mappings are expected to be isospectral to the
electronic Hamiltonian iso .

###  III.5 Effect of IQPE state preparation on noise due to sampling

The results discussed until now employ the SV backend. We now study the
influence of the number of individual bit sampling (also referred to as shots)
for IQPE in our results, by using IBM's QASM backend. It is well known that in
IQPE, the probability of correctly extracting any individual bit is, in
practice, always less than one due to inexact phase expansion. One way to
alleviate this is to repeat the sampling of each bit multiple times, and
obtain the correct phase via the classical error correction technique. If the
probability of extracting a given bit is greater than 0.5, repeated sampling
would reduce the probability of the measurement error exponentially Chernoff .
In the following paragraphs, we investigate the effect of sampling noise in
IQPE, with the state preparation using HF state and the VQE algorithm. The VQE
stage uses the UCCD ansätz, which is set via two different truncation schemes.
In the first scheme, we use an untruncated double excitation space, whereas
for the second, we truncate the double excitation space, and term the latter
as minimally parametrized UCCD ansätz. For all the computations discussed
below, we choose the representative case of 89.8°with 14-bit IQPE based on our
observations from Fig. 3, and vary the number of VQE iterations and the number
of shots together. We use MP2 initial guess for the VQE calculations. Each
such computation is repeated 40 times to build reasonable statistics.

####  III.5.1 State preparation via the HF state

While starting with the HF state as the initial state for IQPE (the
traditionally followed IQPE approach) as shown in Fig. 5(a), the IQPE
algorithm almost always estimates the wrong eigenphase. Starting from the HF
determinant and with as low as 25 shots, the measured eigenphases show severe
spread ($\Delta$, which we define as the difference between the highest and
the lowest observed values) of 1.60 $E_{h}$ when the experiment is repeated 40
times. As we increase the number of shots, the spread of the measured
eigenvalues reduces ($\Delta=$0.985 $E_{h}$ and 0.788 $E_{h}$ with 50 and 100
QASM shots, respectively), and eventually at a sufficiently large number of
shots (10000 shots), leads to a single measured value of the eigenstate, which
corresponds to one of the excited roots. This result is consistent with our
findings from calculations using the SV backend.

####  III.5.2 State preparation via the UCCD ansätz with untruncated double
excitation space

In this section, we briefly discuss our findings when IQPE is repeated with
the initial reference state generated via VQE-UCCD with complete double
excitation manifold. The reference state preparation thus involves a
18-parameter VQE computation. The results shown in Fig. 5(b) reflect that an
appropriately prepared multi-determinantal state with complete inclusion of
two-body excitation operators in the UCCD ansätz can greatly reduce the
requirement of repeating the experiment with large number of QASM shots (in
the IQPE stage of the algorithm). In contrast to the trend observed when the
initial state for IQPE is the HF determinant, the current choice of initial
state displays a much lesser spread in energies. Even with a low number of
QASM shots, but starting with a reasonable state (prepared with 200 VQE
iterations), we observe a spread of only a few $mE_{h}$ ($\Delta=$0.00076
$E_{h}$) around the exact ground state energy (Fig. 5(b)). This is already
comparable to the desired tolerance limit, given the 14-bit IQPE achieves the
target precision of $\pm(2^{-14}\times 2\pi)E_{h}\sim 10^{-4}E_{h}$. With
further increase in the number of shots for IQPE, the estimated eigenphase
shows high bias towards the modal value, which corresponds to the correct
ground state energy that one observes using the SV backend. Increase in the
number of VQE iterations shows little effect on the spread of the observed
energy values, although it enhances the bias towards the modal value,
particularly when the number of QASM samplings per individual bit is
increased. We note that one may obtain the correct eigenphase with much lower
number of VQE iterations and the number of individual bit samplings in IQPE
(QASM shots) when the VQE calculation is performed starting from the MP2 guess
amplitudes, than starting from zero. Nonetheless, even starting from zero
guess amplitudes, upon 40 repetitions of the experiment, one obtains a
sub-$mE_{h}$ spread with very low numbers of QASM shots and VQE iterations.
Thus, we numerically demonstrate that an appropriately prepared state leads to
significant increase in the probability of correctly obtaining the desired
phase to a given precision.

####  III.5.3 The reference state preparation with minimally parametrized UCCD
ansätz

While the UCCD ansätz employed above for the reference state preparation is
quite promising, it may often be plagued by incomplete inclusion of static
correlation in the regions of molecular strong correlation. Furthermore, UCC
in the entire double excitation space involve variational optimization of a
large number of parameters, which results in higher cost. Since the static
correlation often plays a dominant role over dynamical correlation for
generating the wave function having larger overlap with the FCI wavefunction,
one may selectively include only those cluster operators in the unitary ansätz
which are labelled by the chemically active orbitals. This reduces the number
of variational parameters drastically while includes the static correlation
effects to a high extent. We thus propose:

$U_{min}=\exp{(T_{act}(\theta)-T_{act}^{\dagger}(\theta))}$ (16)

where $T_{act}$ is the subset of two-body excitation operators, which are
labelled only by the active orbitals. Below, we shall demonstrate that the
minimally parametrized ansätz proposed above is quite accurate for the purpose
of reference state preparation with efficient inclusion of the static
correlation effect for future study of QPE-based approaches.

Now, we systematically study the effect of noise due to sampling when the
initial reference state is prepared via VQE but within a truncated excitation
space (minimally parametrized ansätz). For this purpose, we again chose to
work in the strong correlation regime ($\beta=$89.8°). We note from Eq. (13)
that the exact ground state wavefunction has large contribution from two
dominant determinants: the HF determinant, and the doubly excited determinant,
where two electrons from the highest occupied molecular orbital (one from the
alpha spinorbital and one from beta spinorbital) are excited to the lowest
unoccupied molecular spinorbitals. Truncation of the excitation space thus
leads to the reduction of the number of VQE parameters to only one. In such a
case, the UCCD ansätz reduces to
$U_{min}=\exp(T_{15}^{26}-(T_{15}^{26})^{\dagger})$, with
$T_{15}^{26}=\theta_{15}^{26}\\{a_{2}^{\dagger}a_{6}^{\dagger}a_{5}a_{1}\\}$,
where $a_{2}^{\dagger}$ $(a_{1})$ and $a_{6}^{\dagger}$ $(a_{5})$ are the
creation (annihilation) operators with up and down spinorbitals in our block
spin arrangement. Note that the iterative optimization of the single VQE
parameter generates the multi-determinantal state which is a linear
superposition of $\ket{00110011}$ and $\ket{01010101}$. The obtained results
from Fig. 5(c) show that IQPE phase estimation starting from the multi-
determinantal state prepared within the truncated active excitation space
shows similar statistics upon repeated experiments as we previously obtained
when started from a multi-determinantal reference state generated via the
complete space VQE optimization.

In Fig. 5(c) bottom panel, we show that with the converged VQE parameters
(where convergence was reached at 24 iterations; note that this is in contrast
to full excitation space case, where convergence is reached at well over 3000
VQE iterations) and as few as 25 QASM shots for IQPE bit measurement, the
final energy shows significantly large spread ($\Delta=0.98E_{h}$) compared to
when the initial reference state was prepared with untruncated doubles
excitation UCCD. However, the modal value (with frequency 25 times out of 40
repetitions) is same as one observes with SV or using a full excitation space
UCCD. As expected, with more number of shots, the spread in the observed
energy decreases, with a sharper gradient than one observes with untruncated
UCCD states. Moreover, the modal value shows stronger bias towards the exact
one observed with SV backend (with frequency 35, 39 and 40 times out of 40
repetitions with 50, 100 and 10000 QASM shots respectively). Furthermore, as
shown in Fig. 5(c) top panel, a similar statistics is observed when the
optimization of the reference state over the single parameter is restricted to
only 10 iterations. Thus, due to the computational gain and given the
constraint in resources in the current era, we recommend for a state
preparation with truncated excitation space UCCD ansätz over a complete space
UCCD for a subsequent evolution by IQPE.

##  IV Conclusion

In this work, we presented a hybrid IQPE over VQE algorithm for the digital
quantum simulation of molecular strong correlation effects in $H_{4}$ on a
circle. While the VQE algorithm using the UCCD ansätz does not perform
satisfactorily around geometries where molecular strong correlation becomes
important, as exemplified by the non-parallelity of the PES as compared to its
FCI counterpart, the hybrid IQPE over VQE algorithm reproduces molecular PES
with very good precision regardless of the electronic complexity. In
particular, we find a whole order of magnitude improvement at strongly
correlated regimes.

With a reasonably high number of initial VQE iterations, one may prepare a
multi-determinantal reference state that has substantial overlap with the
exact ground eigenstate of the many-body Hamiltonian. This results in a much
faster and safer evolution via the subsequent IQPE algorithm. We find that
starting from a well prepared multi-determinantal state, the IQPE algorithm
takes more or less the same number of iterations for an acceptable precision
in energy, across different geometries. While the final precision of the
estimated ground state energy is determined by the number of IQPE iterations,
the success probability of landing in the correct eigenstate is controlled by
the number of VQE iterations. We also checked that one may further tune the
initial guess parameters at the VQE step to further increase the success
probability.

We also studied the effect of sampling noise by using the QASM backend for a
geometry that displays strong correlations in $H_{4}$ on a circle. We reported
results on three distinct cases for the choice of initial state for IQPE,
namely HF state, full excitation space UCC, and minimally parametrized UCC.
The HF result lands on an excited root, and was consistent with the SV
results. Although the trends with increasing number of shots are dissimilar
for the cases that use full excitation space UCC and minimally parametrized
UCC, they produce the same results when VQE is converged and IQPE is run with
a large number of shots. Our findings point that the bias of the modal value
can systematically be tuned with goodness of the prepared state, where
relevant excitations can be chosen with physically/chemically motivated
arguments, and the estimated modal value can show strong bias towards the
exact ground state energy upon repetition of the experiment multiple times.
Our detailed pilot study with the hybrid IQPE over VQE algorithm can serve as
a stepping stone to future explorations of strong correlation effects using
quantum computers.

## Acknowledgements

We are very thankful to Prof. K. Sugisaki for stimulating discussions on the
conceptual underpinnings of the IQPE algorithm. We thank Prof. B. P. Das and
Prof. D. Mukherjee for useful discussions on quantum computing in the initial
stages of this work. Most of the calculations were performed on the
computational facilities at IIT Bombay. The rest of the calculations were
carried out on National Supercomputing Mission's (NSM) computing resource,
`PARAM Siddhi-AI', at C-DAC Pune, which is implemented by C-DAC and supported
by the Ministry of Electronics and Information Technology (MeitY) and
Department of Science and Technology (DST), Government of India. RM thanks
Science and Engineering Research Board DST, Government of India and IRCC, IIT
Bombay for funding.

## References

  * (1) P. W. Shor, Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer, SIAM J. Comp. 26, 1484 (1997).
  * (2) D. Deutsch and R. Jozsa, Rapid Solution of Problems by Quantum Computation, Proc. Royal Soc. London Series A 439, 553 (1992).
  * (3) A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum Algorithm for Linear Systems of Equations, Phys. Rev. Lett. 15, 150502 (2009).
  * (4) S. Y. Yan, Quantum Computational Number Theory, Springer International Publishing Switzerland (2015).
  * (5) P. S. Emani et. al., Quantum Computing at the Frontiers of Biological Sciences, Nat. Methods 18, 701 (2021).
  * (6) R. Orus, S. Mugel, and E. Lizaso, Quantum computing for Finance: Overview and Prospects, Rev. Phys. 4, 100028 (2019).
  * (7) A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Simulated Quantum Computation of Molecular Energies, Science 309, 1704 (2005).
  * (8) J. Du, N. Xu, X. Peng, P. Wang, S. Wu, and D. Lu, NMR Implementation of a Molecular Hydrogen Quantum Simulation with Adiabatic State Preparation, Phys. Rev. Lett. 104, 030502 (2010).
  * (9) L. Veis and J. Pittner, Adiabatic State Preparation Study of Methylene, J. Chem. Phys. 140, 214111 (2014).
  * (10) Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, A. Aspuru-Guzik, Quantum Chemistry in the Age of Quantum Computing, Chem. Rev. 119, 10856 (2019).
  * (11) S. McArdle, S. Endo, A. Aspuru-Guzik, S. Benjamin, X. Yuan, Quantum Computational Chemistry, Rev. Mod. Phys. 92, 15003 (2020).
  * (12) J. D. Whitfield, J. Baltimore, and A. Aspuru-Guzik, Simulation of Electronic Structure Hamiltonians using Quantum Computers , Mol. Phys. 109, 735 (2011).
  * (13) D. S. Abrams and S. Lloyd , Simulation of Many-Body Fermi Systems on a Universal Quantum Computer , Phys. Rev. Lett. 79, 2586 (1997).
  * (14) D. S. Abrams and S. Lloyd , Quantum Algorithm Providing Exponential Speed Increase for Finding Eigenvalues and Eigenvectors , Phys. Rev. Lett. 83, 5162 (1999).
  * (15) A. Peruzzo, J. McClean, P. Shadbolt, M- H. Yung, X- Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien , A Variational Eigenvalue Solver on a Photonic Quantum Processor, Nat. Comm. 5, 4213 (2014).
  * (16) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The Theory of Variational Hybrid Quantum-Classical Algorithms, New J. Phys. 18, 023023 (2016).
  * (17) D. J. Griffiths, Introduction to Quantum Mechanics, 2nd Ed., Chapter 7 (2005).
  * (18) K. Sugisaki, T. Kato, Y. Minato, K. Okuwaki, and Y. Mochizuki, Towards Accurate Description of Chemical Reaction Energetics by Using Variational Quantum Eigensolver: A Case Study of the C2v Quasi-Reaction Pathway of Beryllium Insertion to H2 Molecule, chemRxiv, https://doi.org/10.33774/chemrxiv-2021-w7n78 (2021).
  * (19) T. E. O'Brien, B. Tarasinski, B. M. Terhal, Quantum Phase Estimation of Multiple Eigenvalues for Small-Scale (noisy) Experiments New J. Phys. 21, 023022 (2019)
  * (20) L. Veis and J. Pittner, Quantum Computing Applied to Calculations of Molecular Energies: $\mathrm{CH_{2}}$ Benchmark, J. Chem. Phys. 133, 194106 (2010).
  * (21) H. Wang, S. Kais, A. Aspuru-Guzik and M. R. Hoffmann, Quantum Algorithm for Obtaining the Energy Spectrum of Molecular Systems, Phys. Chem. Chem. Phys. 10, 5388 (2008).
  * (22) N. M. Tubman, C. Mejuto-Zaera, J. M. Epstein, D. Hait, D. S. Levine, W. Huggins, Z. Jiang, J. R. McClean, R. Babbush, M. Head-Gordon, and K. B. Whaley, Postponing the Orthogonality Catastrophe: Efficient State Preparation for Electronic Structure Simulations on Quantum Devices , preprint at https://arxiv.org/abs/1809.05523(2018).
  * (23) K. Sugisaki, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi, and T. Takui, Quantum Chemistry on Quantum Computers: A Method for Preparation of Multiconfigurational Wave Functions on Quantum Computers without Performing Post-Hartree–Fock Calculations, ACS Cent. Sci. 5, 167 (2019).
  * (24) A. Yu. Kitaev, Quantum Measurements and the Abelian Stabilizer Problem, Elec. Coll. Comput. Complex. TR96-003 (1996).
  * (25) M. Dobsicek, G. Johansson, V. Shumeiko, and G. Wendin, Arbitrary Accuracy Iterative Quantum Phase Estimation Algorithm using a Single Ancillary Qubit: A two-qubit benchmark, Phys. Rev. A 76, 030306(R) (2007).
  * (26) J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Strategies for Quantum Computing Molecular Energies using the Unitary Coupled Cluster Ansatz, Quantum Sci. Technol. 4, 014008 (2019).
  * (27) P. M. Q. Cruz, G. Catarina, and R. Gautier, Ronan and J. Fernández-Rossier, Optimizing Quantum Phase Estimation for the Simulation of Hamiltonian Eigenstates, Quantum Sci. Technol. 5, 044005 (2020).
  * (28) Q. Sun et. al., PySCF: The Python‐Based Simulations of Chemistry Framework, Wiley Interdisciplinary Reviews: Computational Molecular Science, 8, e1340 (2017).
  * (29) H. Abraham et. al., Qiskit: An Open-source Framework for Quantum Computing, 10.5281/zenodo.2562110 (2019).
  * (30) J. A. Gomez, T. M. Henderson, and G. E. Scuseria, Recoupling the Singlet- and Triplet-Pairing Channels in Single-Reference Coupled Cluster Theory, J. Chem. Phys 145, 134103 (2016).
  * (31) G. M. J. Barca et. al., Recent developments in the General Atomic and Molecular Electronic Structure System, J. Chem. Phys. 152, 154102 (2020).
  * (32) J.Ivanic and K.Ruedenberg, Identification of Deadwood in Configuration Spaces through General Direct Configuration Interaction, Theoret.Chem.Acc. 106, 339(2001).
  * (33) J.Ivanic, Direct Configuration Interaction and Multi-Configurational Self-Consistent-Field Method for Multiple Active Spaces with Variable Occupancies, Part I.Method, J. Chem. Phys. 119, 9364 (2003).
  * (34) J.Ivanic, Direct Configuration Interaction and Multi-Configurational Self-Consistent-Field Method for Multiple Active Spaces with Variable Occupancies, Part II. Applications, J. Chem. Phys. 119, 9377 (2003).
  * (35) W. J. Hehre, R. F. Stewart, and J. A. Pople, Self‐Consistent Molecular‐Orbital Methods. I. Use of Gaussian Expansions of Slater‐Type Atomic Orbitals, J. Chem. Phys. 51, 2657 (1969).
  * (36) M. Powell, Direct Search Algorithms for Optimization Calculations, Acta Numerica 7,287 (1998).
  * (37) J. T. Seeley, M. J. Richard, and P. J. Love, The Bravyi-Kitaev Transformation for Quantum Computation of Electronic Structure, J. Chem. Phys. 137, 224109 (2012).
  * (38) W. Kutzelnigg, Quantum chemistry in Fock Space. I. The Universal Wave and Energy Operators, J. Chem. Phys. 77, 3081 (1982).
  * (39) A. Dutta and C. D. Sherrill, Full Configuration Interaction Potential Energy Curves for Breaking Bonds to Hydrogen: An Assessment of Single-Reference Correlation Methods, J. Chem. Phys. 118, 1610 (2003).
  * (40) T-C. Yen, R. A. Lang, and A. F. Izmaylov, Exact and Approximate Symmetry Projectors for the Electronic Structure Problem on a Quantum Computer, J. Chem. Phys. 151, 164111 (2019).
  * (41) H. Chernoff, A Career in Statistics, in Past, Present, and Future of Statistics, Ed. by X. Lin, C. Genest, D. L. Banks, G. Molenberghs, D. W. Scott, J-L. Wang, Chapman and Hall/CRC (2014).