# Enhancing initial state overlap through orbital optimization for faster
molecular electronic ground-state energy estimation

Pauline J. Ollitrault  0000-0003-1351-7546<EMAIL_ADDRESS>QC
Ware Corporation, Palo Alto, California 94301, USA  Cristian L. Cortes
0000-0002-1163-2981 QC Ware Corporation, Palo Alto, California 94301, USA
Jérôme F. Gonthier  0000-0002-2933-4085 QC Ware Corporation, Palo Alto,
California 94301, USA  Robert M. Parrish  0000-0002-2406-4741 QC Ware
Corporation, Palo Alto, California 94301, USA  Dario Rocca
0000-0003-2122-6933 QC Ware Corporation, Palo Alto, California 94301, USA
Gian-Luca Anselmetti  0000-0002-8073-3567 Quantum Lab, Boehringer Ingelheim,
55218 Ingelheim am Rhein, Germany  Matthias Degroote  0000-0002-8850-7708
<EMAIL_ADDRESS>Quantum Lab, Boehringer Ingelheim,
55218 Ingelheim am Rhein, Germany  Nikolaj Moll  0000-0001-5645-4667 Quantum
Lab, Boehringer Ingelheim, 55218 Ingelheim am Rhein, Germany  Raffaele
Santagati  0000-0001-9645-0580 Quantum Lab, Boehringer Ingelheim, 55218
Ingelheim am Rhein, Germany  Michael Streif  0000-0002-7509-4748 Quantum Lab,
Boehringer Ingelheim, 55218 Ingelheim am Rhein, Germany

###### Abstract

The phase estimation algorithm is crucial for computing the ground state
energy of a molecular electronic Hamiltonian on a quantum computer. Its
efficiency depends on the overlap between the Hamiltonian’s ground state and
an initial state, which tends to decay exponentially with system size. We
showcase a practical orbital optimization scheme to alleviate this issue.
Applying our method to four iron-sulfur molecules, we achieve a notable
enhancement, up to two orders of magnitude, compared to localized orbitals.
Furthermore, our approach yields improved overlaps in cytochrome P450 enzyme
models.

Determining the electronic ground state is at the heart of making
computational predictions about molecular properties. While traditional
classical algorithms can accurately solve this problem for many simple
molecules, challenges persist for more complex systems [1]. Consequently,
significant research efforts are devoted to developing alternative algorithms
to elucidate the ground-state properties of increasingly complex molecules.
Quantum computing emerges as a particularly promising avenue in this pursuit.
Quantum phase estimation (QPE) [2] stands out as the prototypical algorithm
within this context, offering a pathway to efficiently determine the molecular
electronic ground state.

The total cost of QPE depends on the inverse of the weight (overlap squared),
$p_{0}=|\braket{\Psi_{\text{initial}}}{\Psi_{0}}|^{2}$ of the initial state
$\ket{\Psi_{\text{initial}}}$ with the targeted ground state $\ket{\Psi_{0}}$.
This $p_{0}^{-1}$ scaling reflects the number of times the circuit must be
repeated. Moreover, it was demonstrated that, because of statistical errors in
a single circuit run, the evolution time in each repetition also holds a
$\mathcal{O}(p_{0}^{-1})$ dependence in the worst case, raising the total
algorithm cost to $\mathcal{O}(p_{0}^{-2})$. [3, 4] While recent approaches,
going beyond QPE, achieve better scalings in $p_{0}$ [5], concern remains in
view of the orthogonality catastrophe, i.e., the overlap of a classically
tractable approximate wavefunction with the targeted exact ground state can
decrease exponentially as the system size grows. This phenomenon was adeptly
demonstrated by Lee et al. [6] in their investigation of iron-sulfur cluster
molecules spanning from 2 to 8 metal centers. In the quantum computing
literature, the conventional choice for the initial state is a single Slater
determinant (SD), e.g., the Hartree-Fock (HF) state. Lee et al. [6] showed
that the overlap of a single SD with a density matrix renormalization group
(DMRG) approximation of the ground state of the renowned FeMoco is on the
order of $10^{-7}$, i.e., $10^{7}$ QPE repetitions. State-of-the-art
implementations require $\sim 10^{10}$ Toffoli gates per QPE circuit for this
system, which for optimistic error rate considerations, would take more than a
day to complete [7]. This starkly illustrates the pressing need for enhanced
initial state preparation methods.

Several endeavors have been dedicated to this task, notably introducing
quantum circuits for preparing a configuration state function (CSF), i.e., a
linear combination of SDs defined as an eigenfunction of both the Hamiltonian
and the spin operators, [8, 9, 10] or other sums of SDs [11, 12, 13, 14]
instead of a single one. A similar idea involves the preparation of matrix
product states (MPSs). [15, 16] Alternatively, adiabatic [17] or low-depth
state preparation methods, such as the variational quantum eigensolver [18] or
imaginary time evolution [19, 20], are considered as possible candidates for
preparing a reasonable initial state. They can also be augmented with boosters
to effectively amplify their overlap with the ground state. [21] Moreover, a
compelling proposal for evaluating the quality of initial states via the
computation of an energy distribution has recently emerged. [22]

A key factor in preparing a good initial state, which has attracted less
attention so far, is carefully selecting the molecular orbital (MO) basis.
Previous studies have shown some enhancement in the overlap by working with
the natural orbitals (NOs). [11, 6, 23] Interestingly, it was also
demonstrated that $p_{0}$ converges faster than the energy, suggesting the
potential for achieving accurate overlap estimates classically without
requiring exhaustive knowledge of the ground state. [11] This observation
prompted us to explore the feasibility of leveraging approximate ground states
to refine the orbital basis, potentially surpassing the overlap achieved with
NOs. We are particularly interested in studying this effect in iron-sulfur
clusters as they are found at the active sites of numerous proteins, where
they play key roles in vital biological processes such as photosynthesis or
nitrogen fixation. These molecules exhibit strong correlations, posing
challenges for theoretical comprehension using conventional quantum chemistry
approaches. [24] Hence, they serve as ideal candidates for exploration via
quantum computing methodologies. Notably, they were one of the foci of the
aforementioned work of Lee et al. [6] on initial state preparation. To advance
our understanding in this area, we aim to closely examine their findings and
leverage them as foundational pillars for our research endeavors.

Even within the fault-tolerant quantum computing regime, and particularly in
its early stages when qubit availability remains constrained, the mapping of
fermionic systems to qubits necessitates the careful selection of a specific
set of MOs to encapsulate electronic correlation, known as the active space.
[25, 26, 27] In this study, we avoid delving into the intricate process of
active space selection. This topic is as important as it is complex, with
ample dedicated research within the classical quantum chemistry literature.
[28, 29, 30] The selection of the active space can impact both the value of
$p_{0}$ and the resulting ground state energy, underscoring the need for
meticulous consideration to uphold the chemistry of the target molecule. Here
we leverage the extensive research history surrounding iron-sulfur clusters to
select the active space. [31, 32] Our focus remains on performing orbital
rotations solely within this designated space, maintaining a fixed core energy
throughout our computations.

_Orbital optimization._ Let us consider an approximation,
$\ket{\tilde{\Psi}_{0}}=\sum_{I}^{N_{\text{SD}}}c_{I}\ket{I}$, of the targeted
ground state, $\ket{\Psi_{0}}$ with real $c_{I}$ amplitudes.
$\ket{\tilde{\Psi}_{0}}$ is given to us by a classically tractable procedure
that we will discuss shortly. Our goal today is to produce a single SD,
$\ket{J^{\prime}}$, in a rotated orbital basis that can be used as an
efficient initial state in QPE. Numerically, our task is to maximize the
overlap between $\ket{J^{\prime}}$ and $\ket{\tilde{\Psi}_{0}}$ by optimizing
the MOs, $\\{\ket{\psi^{\prime}_{j}}:0\leq j\leq 2L\\}$, in which
$\ket{J^{\prime}}$ is expressed. Here and throughout we use a spin-restricted
formalism where the spatial parts of the $\alpha$ and $\beta$ orbitals are
identical. $2L$ is the total number of orbitals in the active space since we
include $L$ $\alpha$-spin MOs and $L$ $\beta$-spin MOs. The SDs $\ket{I}$ in
the expansion of $\ket{\tilde{\Psi}_{0}}$ are expressed in the original
computational MO basis $\\{\ket{\psi_{i}}:0\leq i\leq 2L\\}$. As detailed in
appendix A, our task can be realized by optimizing the entries of a $L\times
L$ anti-Hermitian matrix, $\hat{\kappa}$, in order to maximize

$p_{0}\equiv|\braket{J^{\prime}}{\tilde{\Psi}_{0}}|^{2}=\Big{|}\sum_{I}^{N_{\text{SD}}}c_{I}\det\big{(}\hat{\mathrm{M}}^{\alpha}_{J^{\prime}I}\big{)}\det\big{(}\hat{\mathrm{M}}^{\beta}_{J^{\prime}I}\big{)}\Big{|}^{2}$
(1)

where $\hat{\mathrm{M}}=e^{-\hat{\kappa}}$ and
$\hat{\mathrm{M}}^{\sigma}_{J^{\prime}I}$ is a sub-matrix of
$\hat{\mathrm{M}}$ obtained by taking its rows and columns of indices
$J^{\prime}_{\sigma}$ and $I_{\sigma}$, respectively. $J^{\prime}_{\sigma}$
and $I_{\sigma}$ are sets of indices corresponding to occupied orbitals of
spin $\sigma$ in SD $\ket{J^{\prime}}$ and $\ket{I}$, respectively.

_Overlap accuracy._ The $p_{0}$ values we report in this work contain two
sources of error compared to the exact weight that rules the cost of QPE. The
first error comes from the fact that we only have access to an approximate
ground state $\ket{\tilde{\Psi}_{0}}$ and not the exact one, $\ket{\Psi_{0}}$.
Indeed, we recall that the ultimate goal is to run QPE to find the energy of
$\ket{\Psi_{0}}$ and that our only previous knowledge is a classically
tractable approximation of it. The second source of error comes from
truncating the sum in Eq. 1 to make the computation tractable.

Accurate classical quantum chemistry methods rely on a compressed
representation of the wavefunction to avoid storing the full amplitude
statevector. One such example is the use of MPSs as in, e.g., DMRG. In this
case, the wavefunction reads

$\ket{\Psi_{\text{MPS}}}=\sum_{n_{1}...n_{2L}}\sum_{a_{1}...a_{2L-1}}^{D}M^{n_{1}}_{1a_{1}}M^{n_{2}}_{a_{1}a_{2}}...M^{n_{2L}}_{a_{2L-1}1}\ket{n_{1}...n_{2L}}\,.$
(2)

Here, the statevector is expressed as a product of $2L$ rank 3 tensors with
index $n_{i}$ labeling the possible occupations of orbital $i$ and the other
indices running up to $D$, the so-called bond dimension. This bond dimension
determines the accuracy of the representation and the cost of the algorithm.
Another way of compressing the information is to work with a basis of CSFs
rather than SDs. A CSF is a linear combination of SDs that allows to produce
simultaneous eigenfunctions of both the Hamiltonian and the spin operator
$\hat{\text{S}}^{2}$ (whereas the SDs are not necessarily eigenfunctions of
$\hat{\text{S}}^{2}$). On the other hand, Eq. 1 implies a decompression of the
information since $\det(\hat{\mathrm{M}}^{\sigma}_{J^{\prime}I})$ has to be
calculated for all terms in the sum. Therefore, to make the computation
efficient, it is desirable to truncate the sum in Eq. 1 at the cost of
introducing an error in the overlap value.

To formally understand the significance of the above errors, let us define two
general wavefunctions $\ket{\Psi}$ and $\ket{\tilde{\Psi}}$, the residual
$\ket{r}=\ket{\Psi}-\ket{\tilde{\Psi}}$ and its norm
$\epsilon_{r}=||\ket{r}||_{2}=\sqrt{\braket{r}{r}}$. Given a single SD,
$\ket{J^{\prime}}$, with overlaps $\eta=\braket{J^{\prime}}{\Psi}$ and
$\tilde{\eta}=\braket{J^{\prime}}{\tilde{\Psi}}$ we want to find the relation
between $\eta$, $\tilde{\eta}$ and $\epsilon_{r}$. Note that we have,
$\eta-\tilde{\eta}=\braket{J^{\prime}}{r}$. Since
$|\braket{J^{\prime}}{r}|\leq||\ket{J^{\prime}}||_{2}\,\,||\ket{r}||_{2}$, we
find

$|\eta-\tilde{\eta}|\leq\epsilon_{r}\,.$ (3)

The bound on the first source of error, coming from the fact that we only know
an approximation of the true ground state, becomes

$\epsilon_{r,1}=\sqrt{2(1-\braket{\tilde{\Psi}_{0}}{\Psi_{0}})}\,.$ (4)

For the second type of error we note that truncating Eq. 1 is equivalent to
replacing the normalized statevector $\ket{\tilde{\Psi}_{0}}$ by a truncated
version of it $\ket{\tilde{\Psi}_{0,\text{tr}}}$, i.e., keeping only
amplitudes above a certain threshold. This means that
$\braket{\tilde{\Psi}_{0,\text{tr}}}{\tilde{\Psi}_{0}}=\braket{\tilde{\Psi}_{0,\text{tr}}}{\tilde{\Psi}_{0,\text{tr}}}<1$.
Therefore

$\epsilon_{r,2}=\sqrt{1-\braket{\tilde{\Psi}_{0,\text{tr}}}{\tilde{\Psi}_{0,\text{tr}}}}.$
(5)

Clearly, for the type of systems we are targeting with QPE, these bounds are
not tight enough since $\braket{\tilde{\Psi}_{0}}{\Psi_{0}}$ and
$\braket{\tilde{\Psi}_{0,\text{tr}}}{\tilde{\Psi}_{0,\text{tr}}}$ would tend
to zero. However, if the orbital bases of $\ket{J^{\prime}}$ and the reference
state are similar, then $e^{-\hat{\kappa}}$ is close to the identity matrix.
In other words, $\ket{J^{\prime}}$ overlaps with a limited part of the
reference state, i.e. only a few terms in Eq. 1 survive. This means that we
can rewrite $\eta=\braket{J^{\prime}}{\Psi_{\text{tr}}}$ and
$\tilde{\eta}=\braket{J^{\prime}}{\tilde{\Psi}_{\text{tr}}}$ where
$\ket{\Psi_{\text{tr}}}$ ($\ket{\tilde{\Psi}_{\text{tr}}}$) is a truncated
$\ket{\Psi}$ ($\ket{\tilde{\Psi}}$) with norm much below 1. This directly
reduces the bound to
$\epsilon_{r}=\sqrt{\braket{\Psi_{\text{tr}}}{\Psi_{\text{tr}}}+\braket{\tilde{\Psi}_{\text{tr}}}{\tilde{\Psi}_{\text{tr}}}-2\braket{\Psi_{\text{tr}}}{\tilde{\Psi}_{\text{tr}}}}$.
In the next paragraphs we demonstrate this bound in several numerical examples
before showing the performance of orbital optimization for iron-sulfur
clusters.

_Numerical results._ To calculate the overlaps of interest and optimize the
orbital basis, we need to find a reasonable approximation of the ground state.
In this regard, we employ either spin-adapted DMRG from Block 2 [33, 34] or
semistochastic heat bath configuration interaction (SHCI) from Dice [35, 36].
While the wavefunction resulting from SHCI can be directly used to calculate
the overlap according to Eq. 1, the MPS ensued by DMRG, expressed in CSF
basis, needs to be post-processed. In this case, we employ a sampling routine
to first save all CSFs with absolute amplitude above $t_{\text{CSF}}$. Then,
we transform them to a superposition of SDs following the algorithm of Ref.
[37], as detailed in Appendix B, with a threshold $t_{\text{SD}}$ under which
the SDs are discarded. We denote the resulting statevector as
$\ket{\tilde{\Psi}_{0}}$. To obtain a new orbital basis efficiently, we first
sort $\ket{\tilde{\Psi}_{0}}$ by order of absolute amplitude, from the largest
to the smallest. We choose $\ket{J^{\prime}}$ as the SD with the largest
amplitude in $\ket{\tilde{\Psi}_{0}}$. We optimize the orbital basis of
$\ket{J^{\prime}}$ using the GradientDescent optimizer of JAXopt [38] to
minimize $f(\hat{\kappa})=1-|\braket{J^{\prime}}{\tilde{\Psi}_{0}}|^{2}$
according to Eq. 1, with a given value $N_{\text{SD}}^{\text{opt}}$ to
truncate the sum. We always start the optimization by setting all
$\hat{\kappa}$ entries to zero, i.e., we start by expressing
$\ket{J^{\prime}}$ in the same basis as $\ket{\tilde{\Psi}_{0}}$. All orbital
rotations are strictly applied inside the active space. The core energy always
remains unchanged.

_Acenes._ As an initial exploration, we study the family of acenes, which are
polycyclic aromatic hydrocarbons with a linear arrangement of fused benzene
rings (see Figure 1). Acenes have unique electronic properties, including high
electron mobility and excellent charge transport characteristics. As they
exhibit semiconducting behavior, they are valuable in developing organic
semiconductors used in various electronic devices, such as organic field-
effect transistors, organic light-emitting diodes, and organic photovoltaics.
They are strongly-correlated systems and, therefore, serve as interesting
models for our study.

Following the work of Sharma et al. [39] (using their geometries, see Appendix
C) we obtain the MOs with B3LYP/6-31G(d,p). The active space corresponds to
the $\pi$ electrons distributed in all valence $\pi$ orbitals leading to
$2(2n+1)$ electrons in $2(2n+1)$ orbitals, where $n$ is the number of rings.

A preliminary study on naphtalene ($n=2$), detailed in Appendix C, illustrates
our previous discussion on the first type of error, i.e., it indicates that
the true $p_{0}$ error is substantially lower than $\epsilon_{r,1}$. In
particular, we obtain two approximations of the singlet ground state,
$\ket{\tilde{\Psi}_{0}^{\text{MO}}}$ and $\ket{\tilde{\Psi}_{0}^{\text{NO}}}$,
in the MO and NO basis, respectively. We show that, given the proximity
between the two orbital bases,
$|\braket{J^{\prime}}{\tilde{\Psi}_{0}^{\text{MO}}}-\braket{J^{\prime}}{\tilde{\Psi}_{0}^{\text{NO}}}|\ll||\ket{\tilde{\Psi}_{0}^{\text{MO}}}-\ket{\tilde{\Psi}_{0}^{\text{NO}}}||_{2}$.
This is repeated for ground states approximated with both spin-adapted DMRG
and SHCI, and for several single SD $\ket{J^{\prime}}$, leading, in all cases,
to the same observation.

Figure 1: Convergence of $p_{0}$ with the number of SDs included in the
reference statevector, $\ket{\tilde{\Psi}_{0}}$. The full lines correspond to
the error between $p_{0}$ calculated using $N_{\text{SD}}$ with respect to the
value of the overlap obtained from all sampled SDs. The error bound,
$\epsilon_{r}$, is also shown in dashed lines. The results are given for
acenes of different sizes $n$ and for singlet states.

We also approximate the ground states of larger acenes with $3\leq n\leq 6$,
in both singlet and triplet states, with spin-adapted DMRG. We obtain, in each
case, an optimized orbital basis and show that we can always increase $p_{0}$.
In all these cases the resulting enhancement is, however, very modest.
Nonetheless, this allows us to make several interesting observations. First,
we obtain the same basis by running the optimization using 10% or 50% of the
total number of SDs in the reference statevector as a value for
$N^{\text{opt}}_{\text{SD}}$ to truncate the sum in Eq. 1. We also show that
for all systems studied here, once the orbital basis is optimized and fixed,
the resulting value of $p_{0}$ converges quickly with $N_{\text{SD}}$. This is
highlighted in Figure 1 for the singlet states. We see that, for each system
size, the error made in the calculated $p_{0}$ is substantially lower than the
residual norm (corresponding to Eq. 5). We also show that the optimized
orbital basis differs from the NOs in all cases. The optimized basis is closer
to the starting MO basis compared to the NOs. Further detail is provided in
Appendix C.

_Iron-sulfur clusters._ As mentioned earlier, one of our focuses is on
determining whether optimizing the orbital basis can significantly enhance the
low overlaps observed in Ref. [6]. Here, we focus on the clusters with the
number of iron atoms $N_{\text{Fe}}=2$ and $N_{\text{Fe}}=4$. The active
spaces are then constructed as in Ref. [6], from the full valence space of the
Fe $3d$, S $3p$ and bonding ligand orbitals around each Fe atom. The resulting
active space sizes are shown in Figure 2. Further details, such as the
oxidation and spin states, are given in Appendix Table 4. In the rest of this
paragraph, we will denote these initial orbitals as MOs. Although these
orbitals deviate from the canonical HF MOs, we adopt this terminology to
maintain simplicity in reading and refrain from introducing additional
acronyms.

We approximate the ground states of the four clusters with spin-adapted DMRG
at bond dimension $D=8000$. We employ a singlet embedding [34] scheme for
[Fe2S2]-3. The resulting energies, given in Table 5, are in excellent
agreement with Ref. [6]. We use the converged state to optimize the orbitals
which we refer to as the OPT basis. In all cases this allows us to increase
the overlap by more than one order of magnitude, as shown in Figure 2.
Moreover, for all these clusters, optimizing the orbitals allows us to reach a
higher $p_{0}$ than starting with a CSF in the original MO basis. Further
detail is given in Appendix D.

Figure 2: (Top) Weight of a single SD or CSF expressed either in the original
MO basis or in the optimized (OPT) one for the four iron-sulfur clusters.
(Bottom) Improvement in $p_{0}$ obtained by starting either from a CSF in the
original MO basis or from a single SD in the optimized (OPT) basis compared to
starting with a single SD in the original MO basis. The $x$-axis labels give
the cluster name and its active space size.

To show that the $p_{0}$ values reported in Figure 2 are accurate, let us
focus on [Fe4S4]-2. Its structure is shown in Figure 3a. First, we compare the
weights obtained from different approximations of the ground state. Reference
states, $\ket{\tilde{\Psi}_{0}(D)}$ are found with spin-adapted DMRG at bond
dimensions $D\in\\{200,400,600,800,1000\\}$. We find an optimized orbital
basis, OPT$D$, using each of the $\ket{\tilde{\Psi}_{0}(D)}$. In Figure 3b
(top) we show $\Delta E=|E(8000)-E(D)|$, the difference between the energies
of $\ket{\tilde{\Psi}_{0}(8000)}$ and $\ket{\tilde{\Psi}_{0}(D)}$ as well as
the most important amplitude squared in $\ket{\tilde{\Psi}_{0}(D)}$. After the
optimization, we calculate the final value of
$|\braket{J^{\prime}}{\tilde{\Psi}_{0}(8000)}|^{2}$ with $\ket{J^{\prime}}$ in
OPT$D$, see Figure 3b (middle). Interestingly, these results clearly show that
statevectors obtained at a smaller $D$, while inaccurate in energy, already
allow to find an orbital basis which substantially improves $p_{0}$. In
addition, we study the convergence of
$|\braket{J^{\prime}}{\tilde{\Psi}_{0}(8000)}|^{2}$, where $\ket{J^{\prime}}$
is in the OPT(8000) basis, with the truncation $N_{\text{SD}}$ of the sum in
Eq. 1. As can be seen in Figure 3b (bottom), $p_{0}$ converges very quickly
with $N_{\text{SD}}$. In Appendix D, we show the same convergence behavior for
all four clusters. These observations indicate that having a qualitatively
correct description of the dominant part of the ground state wavefunction is
enough to find an orbital basis in which $p_{0}$ is improved.

Figure 3: a. Structure of [Fe4S4]-2. All the results in this Figure are
obtained for this cluster. b. (Top) Convergence of the energy and weight with
the DMRG bond dimension, $D$. $\Delta E=|E(8000)-E(D)|$ is the difference
between the reference energy and the one obtained at reduced $D$. Circles
represent the most important absolute amplitude squared in
$\ket{\tilde{\Psi}_{0}(D)}$. (Middle) Circles correspond to
$|\braket{J^{\prime}}{\tilde{\Psi}_{0}(8000)}|^{2}$ where $\ket{J^{\prime}}$
is in the orbital basis optimized using $\ket{\tilde{\Psi}_{0}(D)}$. The
dashed line is $p_{0}$(OPT), i.e.,
$|\braket{J^{\prime}}{\tilde{\Psi}_{0}(8000)}|^{2}$ where $\ket{J^{\prime}}$
is in the orbital basis optimized using $\ket{\tilde{\Psi}_{0}(8000)}$.
(Bottom) Convergence of $p_{0}$(OPT) with the number of SDs included in the
reference statevector, $\ket{\tilde{\Psi}_{0}(8000)}$. c. Overlap matrix
between the starting MOs and either the NOs or the optimized orbitals (OPT).

Finally, we compare the NO and optimized basis. The NOs are obtained by
diagonalizing the 1-particle density matrix of $\ket{\tilde{\Psi}_{0}(8000)}$.
In Figure 3c, we show the overlap matrix between the starting MOs and the NOs
as well as between the starting MOs and the optimized orbitals. As we observed
in the acene case, the NOs and optimized orbitals are different, and the
latter are much closer to the initial MOs. Here, the direct calculation of
$p_{0}$(NO) between the Aufbau filled SD in the NO basis, and the MO reference
state leads to values much lower than $p_{0}$(MO). However, the density of
non-diagonal elements in the MO-NO overlap matrix is large. This translates to
having more non-vanishing terms in the sum of Eq. 1. We therefore expect the
errors defined in Eqs. 4 and 5 to be large, hence resulting in unreliable
overlap estimates. To allow for a more fair comparison between the OPT and NO
bases, we run DMRG again in the NO basis and get $p_{0}$(NO) as the amplitude
squared of the most important determinant in the resulting state. This is only
done for the two $N_{\text{Fe}}=2$ clusters, as we can afford a large enough
bond dimension to compensate for the loss of locality in the orbitals. As the
energy and weight converge at $D=500$, we use this bond dimension to obtain
our states in the NO basis. The resulting energies are $-5092.8710$ Ha and
$-5092.7187$ Ha in agreement with their counterparts in the MO basis. We
obtain $p_{0}$(NO) values of $3.35\cdot 10^{-2}$ and $2.46\cdot 10^{-2}$ for
[Fe2S2]-2 and [Fe2S2]-3, respectively. In comparison the $p_{0}$(OPT) values
are $4.53\cdot 10^{-2}$ and $1.80\cdot 10^{-1}$ for [Fe2S2]-2 and [Fe2S2]-3,
respectively. In both cases we obtain better weights in the optimized basis.
_Cytochrome P450._ As a final example, we apply our method to a cytochrome
P450 enzyme, a system of notable relevance to the pharmaceutical industry
[40]. It has been recently demonstrated that such systems exhibit a complex
electronic structure, thereby positioning them as promising candidates for
exploration with future quantum computers [41]. Using our method, we
successfully enhanced the initial overlap of the $S=5/2$ system described in
an 43-orbital active space from $p_{0}=0.53$ to $p_{0}=0.58$. Similar results
are obtained for different active space sizes and spin states. We report
further details in Appendix E.
_Discussion and Conclusion._ In this work, we used orbital rotations to
optimize the overlap between a single SD and the ground state of molecular
electronic Hamiltonians. Specifically, we showed that for four iron-sulfur
clusters, optimizing the orbital basis can lead up to $\sim$2 orders of
magnitude improvement in comparison to previously predicted $p_{0}$ values
[6]. Interestingly, expressing a single SD in an optimized basis also improves
upon the previously reported $p_{0}$ values between a CSF and the targeted
state. Furthermore, the application of our method to the cytochrome P450
enzyme which already displays a substantially larger overlap than the iron-
sulfur clusters resulted in a further 10% increase in overlap.

Our optimization scheme leads to orbitals that resemble the starting ones,
i.e., the ones in which the targeted state is approximated. This translates to
a good convergence of the weight, i.e., truncating the sum in Eq. 1 to
relatively low $N_{\text{SD}}$ compared to the total number of SDs in the
approximated state already leads to accurate $p_{0}$. Moreover, even states
with largely unconverged energies yield improved orbital bases. This indicates
that a qualitatively correct description of the dominant part of the
wavefunction is sufficient to substantially enhance $p_{0}$.

The optimized orbitals also differ from the NOs. This is interesting as the
NOs are known to minimize the orbital correlation and the entropy of the
amplitudes of the statevector. [42] We do not expect the optimized orbital
basis to lead compact states, but this is irrelevant from the quantum
computing perspective, where only the overlap between the initial and the
target state matters.

While the shown method delays the orthogonality catastrophe, it does not fully
solve it. The weights calculated for the iron-sulfur clusters still
demonstrate rapidly decaying behavior. Therefore, the methods discussed in
this work must be combined with additional techniques to reach $p_{0}$ values
closer to unity. One possibility would be to optimize the orbitals, such as to
maximize the overlap of a CSF with the targeted state. However, classically
this becomes progressively more expensive as the number of SDs within a CSF
grows exponentially with the number of its unpaired electrons. While rotating
a general quantum state into a different orbital basis is a difficult task for
a classical computer, it can be done efficiently on a quantum computer by
applying a tessellation of Givens rotation [43]. Such basis rotations were
leveraged to reduce the runtime of QPE by encoding factorized Hamiltonians
[44, 45, 46, 7, 41, 47, 48]. More recently, work from Marti-Dafcik et al. [49]
introduces a compact representation of the wavefunction where different CSFs
are expressed in distinct orbital bases. As highlighted in their manuscript,
this ansatz is particularly well suited for preparing high overlap initial
states for QPE. Similarly, the orbital optimization presented here could be
efficiently performed with a quantum processor.

## References

  * Santagati  _et al._ [2024] R. Santagati, A. Aspuru-Guzik, R. Babbush, M. Degroote, L. González, E. Kyoseva, N. Moll, M. Oppel, R. M. Parrish, N. C. Rubin, M. Streif, C. S. Tautermann, H. Weiss, N. Wiebe, and C. Utschig-Utschig, Drug design on quantum computers, Nature Physics 10.1038/s41567-024-02411-5 (2024).
  * Kitaev  _et al._ [2002] A. Y. Kitaev, A. Shen, and M. N. Vyalyi,  _Classical and quantum computation_, 47 (American Mathematical Society, 2002).
  * Ge  _et al._ [2019] Y. Ge, J. Tura, and J. I. Cirac, Faster ground state preparation and high-precision ground energy estimation with fewer qubits, Journal of Mathematical Physics 60, https://doi.org/10.1063/1.5027484 (2019).
  * Lin and Tong [2022] L. Lin and Y. Tong, Heisenberg-limited ground-state energy estimation for early fault-tolerant quantum computers, PRX Quantum 3, 010318 (2022).
  * Dong  _et al._ [2022] Y. Dong, L. Lin, and Y. Tong, Ground-state preparation and energy estimation on early fault-tolerant quantum computers via quantum eigenvalue transformation of unitary matrices, PRX Quantum 3, 040305 (2022).
  * Lee  _et al._ [2023] S. Lee, J. Lee, H. Zhai, Y. Tong, A. M. Dalzell, A. Kumar, P. Helms, J. Gray, Z.-H. Cui, W. Liu, _et al._ , Evaluating the evidence for exponential quantum advantage in ground-state quantum chemistry, Nature Communications 14, 1952 (2023).
  * Lee  _et al._ [2021] J. Lee, D. W. Berry, C. Gidney, W. J. Huggins, J. R. McClean, N. Wiebe, and R. Babbush, Even more efficient quantum computations of chemistry through tensor hypercontraction, PRX Quantum 2, 030305 (2021).
  * Carbone  _et al._ [2022] A. Carbone, D. E. Galli, M. Motta, and B. Jones, Quantum circuits for the preparation of spin eigenfunctions on quantum computers, Symmetry 14, 624 (2022).
  * Sugisaki  _et al._ [2016] K. Sugisaki, S. Yamamoto, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi, and T. Takui, Quantum chemistry on quantum computers: a polynomial-time quantum algorithm for constructing the wave functions of open-shell molecules, The Journal of Physical Chemistry A 120, 6459 (2016).
  * Sugisaki  _et al._ [2019] K. Sugisaki, S. Yamamoto, S. Nakazawa, K. Toyota, K. Sato, D. Shiomi, and T. Takui, Open shell electronic state calculations on quantum computers: A quantum circuit for the preparation of configuration state functions based on serber construction, Chemical Physics Letters 737, 100002 (2019).
  * Tubman  _et al._ [2018] 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, arXiv preprint arXiv:1809.05523 https://doi.org/10.48550/arXiv.1809.05523 (2018).
  * Ortiz  _et al._ [2001] G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Quantum algorithms for fermionic simulations, Physical Review A 64, 022319 (2001).
  * Parrish  _et al._ [2019] R. M. Parrish, E. G. Hohenstein, P. L. McMahon, and T. J. Martínez, Quantum computation of electronic transitions using a variational quantum eigensolver, Physical Review Letters 122, 230401 (2019).
  * Babbush  _et al._ [2015] R. Babbush, J. McClean, D. Wecker, A. Aspuru-Guzik, and N. Wiebe, Chemical basis of trotter-suzuki errors in quantum chemistry simulation, Physical Review A 91, 022311 (2015).
  * Schön  _et al._ [2005] C. Schön, E. Solano, F. Verstraete, J. I. Cirac, and M. M. Wolf, Sequential generation of entangled multiqubit states, Physical Review Letters 95, 110503 (2005).
  * Malz  _et al._ [2024] D. Malz, G. Styliaris, Z.-Y. Wei, and J. I. Cirac, Preparation of matrix product states with log-depth quantum circuits, Physical Review Letters 132, 040404 (2024).
  * Aspuru-Guzik  _et al._ [2005] A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon, Simulated quantum computation of molecular energies, Science 309, 1704 (2005).
  * Peruzzo  _et al._ [2014] 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, Nature Communications 5, 4213 (2014).
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nature Physics 16, 205 (2020).
  * McArdle  _et al._ [2019] S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, Variational ansatz-based quantum simulation of imaginary time evolution, npj Quantum Information 5, 75 (2019).
  * Wang  _et al._ [2022] G. Wang, S. Sim, and P. D. Johnson, State preparation boosters for early fault-tolerant quantum computation, Quantum 6, 829 (2022).
  * Fomichev  _et al._ [2023] S. Fomichev, K. Hejazi, M. S. Zini, M. Kiser, J. F. Morales, P. A. M. Casares, A. Delgado, J. Huh, A.-C. Voigt, J. E. Mueller, _et al._ , Initial state preparation for quantum chemistry on quantum computers, arXiv preprint arXiv:2310.18410 https://doi.org/10.48550/arXiv.2310.18410 (2023).
  * Ratini  _et al._ [2023] L. Ratini, C. Capecci, and L. Guidoni, Natural orbitals and sparsity of quantum mutual information, arXiv preprint arXiv:2308.08056 https://doi.org/10.48550/arXiv.2308.08056 (2023).
  * Reiher  _et al._ [2017] M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, and M. Troyer, Elucidating reaction mechanisms on quantum computers, Proceedings of the National Academy of Sciences 114, 7555 (2017).
  * Takeshita  _et al._ [2020] T. Takeshita, N. C. Rubin, Z. Jiang, E. Lee, R. Babbush, and J. R. McClean, Increasing the representation accuracy of quantum simulations of chemistry without extra quantum resources, Physical Review X 10, 011004 (2020).
  * Vorwerk  _et al._ [2022] C. Vorwerk, N. Sheng, M. Govoni, B. Huang, and G. Galli, Quantum embedding theories to simulate condensed systems on quantum computers, Nature Computational Science 2, 424 (2022).
  * Rossmannek  _et al._ [2023] M. Rossmannek, F. Pavosevic, A. Rubio, and I. Tavernelli, Quantum embedding method for the simulation of strongly correlated systems on quantum computers, The Journal of Physical Chemistry Letters 14, 3491 (2023).
  * Roos  _et al._ [1980] B. O. Roos, P. R. Taylor, and P. E. Sigbahn, A complete active space scf method (casscf) using a density matrix formulated super-ci approach, Chemical Physics 48, 157 (1980).
  * Stein and Reiher [2019] C. J. Stein and M. Reiher, autocas: A program for fully automated multiconfigurational calculations, Journal of Computational Chemistry 40, 2216 (2019).
  * Veryazov  _et al._ [2011] V. Veryazov, P. Å. Malmqvist, and B. O. Roos, How to select active space for multiconfigurational quantum chemistry?, International Journal of Quantum Chemistry 111, 3329 (2011).
  * Li  _et al._ [2019] Z. Li, S. Guo, Q. Sun, and G. K.-L. Chan, Electronic landscape of the p-cluster of nitrogenase as revealed through many-electron quantum wavefunction simulations, Nature Chemistry 11, 1026 (2019).
  * Li and Chan [2017] Z. Li and G. K.-L. Chan, Spin-projected matrix product states: Versatile tool for strongly correlated systems, Journal of Chemical Theory and Computation 13, 2681 (2017).
  * Zhai  _et al._ [2023] H. Zhai, H. R. Larsson, S. Lee, Z.-H. Cui, T. Zhu, C. Sun, L. Peng, R. Peng, K. Liao, J. Tölle, J. Yang, S. Li, and G. K.-L. Chan, Block2: A comprehensive open source framework to develop and apply state-of-the-art DMRG algorithms in electronic structure and beyond, The Journal of Chemical Physics 159, 234801 (2023).
  * Sharma and Chan [2012] S. Sharma and G. K. Chan, Spin-adapted density matrix renormalization group algorithms for quantum chemistry, The Journal of Chemical Physics 136, https://doi.org/10.1063/1.3695642 (2012).
  * Holmes  _et al._ [2016] A. A. Holmes, N. M. Tubman, and C. J. Umrigar, Heat-bath configuration interaction: An efficient selected configuration interaction algorithm inspired by heat-bath sampling, Journal of Chemical Theory and Computation 12, 3674 (2016).
  * Sharma  _et al._ [2017] S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi, and C. J. Umrigar, Semistochastic heat-bath configuration interaction method: Selected configuration interaction with semistochastic perturbation theory, Journal of Chemical Theory and Computation 13, 1595 (2017).
  * Fales and Martínez [2020] B. S. Fales and T. J. Martínez, Fast transformations between configuration state function and slater determinant bases for direct configuration interaction, The Journal of Chemical Physics 152, https://doi.org/10.1063/5.0005155 (2020).
  * Blondel  _et al._ [2021] M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares-López, F. Pedregosa, and J.-P. Vert, Efficient and modular implicit differentiation, arXiv preprint arXiv:2105.15183 https://doi.org/10.48550/arXiv.2105.15183 (2021).
  * Sharma  _et al._ [2019] P. Sharma, V. Bernales, S. Knecht, D. G. Truhlar, and L. Gagliardi, Density matrix renormalization group pair-density functional theory (dmrg-pdft): singlet–triplet gaps in polyacenes and polyacetylenes, Chemical Science 10, 1716 (2019).
  * Nelson [2018] D. R. Nelson, Cytochrome p450 diversity in the tree of life, Biochimica et Biophysica Acta (BBA)-Proteins and Proteomics 1866, 141 (2018).
  * Goings  _et al._ [2022a] J. J. Goings, A. White, J. Lee, C. S. Tautermann, M. Degroote, C. Gidney, T. Shiozaki, R. Babbush, and N. C. Rubin, Reliably assessing the electronic structure of cytochrome p450 on today’s classical computers and tomorrow’s quantum computers, Proceedings of the National Academy of Sciences 119, e2203533119 (2022a).
  * Aliverti-Piuri  _et al._ [2024] D. Aliverti-Piuri, K. Chatterjee, L. Ding, K. Liao, J. Liebert, and C. Schilling, What can quantum information theory offer to quantum chemistry?, arXiv preprint arXiv:2403.08045 https://doi.org/10.48550/arXiv.2403.08045 (2024).
  * Anselmetti  _et al._ [2021] G.-L. R. Anselmetti, D. Wierichs, C. Gogolin, and R. M. Parrish, Local, expressive, quantum-number-preserving vqe ansätze for fermionic systems, New Journal of Physics 23, 113010 (2021).
  * Berry  _et al._ [2019] D. W. Berry, C. Gidney, M. Motta, J. R. McClean, and R. Babbush, Qubitization of arbitrary basis quantum chemistry leveraging sparsity and low rank factorization, Quantum 3, 208 (2019).
  * von Burg  _et al._ [2021] V. von Burg, G. H. Low, T. Häner, D. S. Steiger, M. Reiher, M. Roetteler, and M. Troyer, Quantum computing enhanced computational catalysis, Physical Review Research 3, 033055 (2021).
  * Motta  _et al._ [2021] M. Motta, E. Ye, J. R. McClean, Z. Li, A. J. Minnich, R. Babbush, and G. K.-L. Chan, Low rank representations for quantum simulation of electronic structure, npj Quantum Information 7, 83 (2021).
  * Oumarou  _et al._ [2022] O. Oumarou, M. Scheurer, R. M. Parrish, E. G. Hohenstein, and C. Gogolin, Accelerating quantum computations of chemistry through regularized compressed double factorization, arXiv preprint arXiv:2212.07957 https://doi.org/10.48550/arXiv.2212.07957 (2022).
  * Rocca  _et al._ [2024] D. Rocca, C. L. Cortes, J. Gonthier, P. J. Ollitrault, R. M. Parrish, G.-L. Anselmetti, M. Degroote, N. Moll, R. Santagati, and M. Streif, Reducing the runtime of fault-tolerant quantum simulations in chemistry through symmetry-compressed double factorization, arXiv preprint arXiv:2403.03502 https://doi.org/10.48550/arXiv.2403.03502 (2024).
  * Marti-Dafcik  _et al._ [2024] D. Marti-Dafcik, N. Lee, H. G. Burton, and D. P. Tew, Spin-coupled molecular orbitals: chemical intuition meets quantum chemistry, arXiv preprint arXiv:2402.08858 https://doi.org/10.48550/arXiv.2402.08858 (2024).
  * Ollitrault  _et al._ [2024] P. J. Ollitrault, C. L. Cortes, J. F. Gonthier, R. M. Parrish, D. Rocca, G.-L. Anselmetti, M. Degroote, N. Moll, R. Santagati, and M. Strei, Enhancing initial state overlap through orbital optimization for faster molecular electronic ground-state energy estimation (2024).
  * Sun  _et al._ [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, _et al._ , Pyscf: the python-based simulations of chemistry framework, Wiley Interdisciplinary Reviews: Computational Molecular Science 8, e1340 (2018).
  * Sharma  _et al._ [2014] S. Sharma, K. Sivalingam, F. Neese, and G. K.-L. Chan, Low-energy spectrum of iron–sulfur clusters directly from many-particle quantum mechanics, Nature Chemistry 6, 927 (2014).
  * Goings  _et al._ [2022b] J. J. Goings, A. White, C. S. Tautermann, M. Degroote, C. Gidney, T. Shiozaki, R. Babbush, and N. C. Rubin, Data for ”Reliably assessing the electronic structure of cytochrome P450 on today’s classical computers and tomorrow’s quantum computers” (2022b).

##  Appendix A Orbital optimization

An orbital basis is a collection of single-particle spin MOs defined by the
product of a spatial and a spin part

$\ket{\psi_{i_{\sigma}}}=\ket{\bar{\psi}_{i_{\sigma}}}\ket{\sigma}\,.$ (A.1)

The spatial part is expressed, in turn, as a linear combination of basis
functions, or atomic orbitals (AOs),

$\ket{\bar{\psi}_{i_{\sigma}}}=\sum_{p}^{N_{\text{AO}}}\mathcal{C}_{pi_{\sigma}}\ket{\chi_{p}}$
(A.2)

while the spin part is either $\ket{\alpha}$ or $\ket{\beta}$. There are then
$L_{\alpha}+L_{\beta}=2L$ MOs in the basis. We can write the overlap between
two MOs as

$\braket{\psi^{\prime}_{j_{\tau}}}{\psi_{i_{\sigma}}}=\delta_{\tau\sigma}\braket{\bar{\psi}^{\prime}_{j_{\tau}}}{\bar{\psi}_{i_{\sigma}}}=\delta_{\tau\sigma}\sum_{pq}^{N_{\text{AO}}}\mathcal{C}^{\prime}_{pj_{\tau}}\mathcal{C}_{qi_{\sigma}}S_{pq}$
(A.3)

where $S_{pq}=\braket{\chi_{p}}{\chi_{q}}$. The MO basis,
$\\{\ket{\psi^{\prime}_{i}}\\}$ can be expressed as a rotation of the
$\\{\ket{\psi_{i}}\\}$ basis. In other words, there exists an anti-Hermitian
matrix $\hat{\kappa}$ such that

$\mathcal{C}^{\prime}_{pj_{\tau}}=\sum_{k_{\tau}}^{L_{\tau}}[e^{-\hat{\kappa}}]_{j_{\tau}k_{\tau}}\mathcal{C}_{pk_{\tau}}\,.$
(A.4)

Hence, the overlap between two MOs becomes

$\braket{\psi^{\prime}_{j_{\tau}}}{\psi_{i_{\sigma}}}=\delta_{\tau\sigma}\sum_{pq}^{N_{\text{AO}}}\sum^{L_{\tau}}_{k_{\tau}}[e^{-\hat{\kappa}}]_{j_{\tau}k_{\tau}}\mathcal{C}_{pk_{\tau}}S_{pq}\mathcal{C}_{qi_{\sigma}}=\delta_{\tau\sigma}\sum^{L_{\tau}}_{k_{\tau}}[e^{-\hat{\kappa}}]_{j_{\tau}k_{\tau}}\mathcal{S}_{k_{\tau}i_{\sigma}}.$
(A.5)

$\hat{\mathcal{S}}$ is the matrix containing the MOs $\\{\ket{\psi_{i}}\\}$
mutual overlaps. Since the MOs are orthogonal, $\hat{\mathcal{S}}$ is simply
the identity matrix and hence, Eq. A.5 simplifies to

$\braket{\psi^{\prime}_{j_{\tau}}}{\psi_{i_{\sigma}}}=\left\\{\begin{array}[]{ll}[e^{-\hat{\kappa}}]_{j_{\sigma}i_{\sigma}}&\mbox{if
}\tau=\sigma\\\ 0&\mbox{otherwise.}\end{array}\right.$ (A.6)

A SD is an antisymmetrized product of these orbitals. For a
$N=N_{\alpha}+N_{\beta}$ particle system, the $N$ MOs included in the SD are
given by two sets of indices $I_{\alpha}=\\{i_{\alpha}^{\nu}:1\leq\nu\leq
N_{\alpha}\\}$ and $I_{\beta}=\\{i_{\beta}^{\mu}:1\leq\mu\leq N_{\beta}\\}$.
The corresponding SD can then be defined as

$\ket{I}\equiv\mathcal{A}\big{[}\prod_{\nu}\ket{\psi_{i_{\alpha}^{\nu}}}\prod_{\mu}\ket{\psi_{i_{\beta}^{\mu}}}\big{]}$
(A.7)

where $\mathcal{A}$ is the antisymmetrizer. A statevector is then defined as a
linear combination of SDs

$\ket{\Psi}\equiv\sum_{I}^{N_{\text{SD}}}c_{I}\ket{I}$ (A.8)

with real $c_{I}$ coefficients. The overlap between $\ket{\Psi}$ and a SD
defined in another orbital basis,

$\ket{J^{\prime}}\equiv\mathcal{A}\big{[}\prod_{\nu}\ket{\psi^{\prime}_{j_{\alpha}^{\nu}}}\prod_{\mu}\ket{\psi^{\prime}_{j_{\beta}^{\mu}}}\big{]},$
(A.9)

becomes

$\braket{J^{\prime}}{\Psi}\equiv\sum_{I}^{N_{\text{SD}}}c_{I}\braket{J^{\prime}}{I}$
(A.10)

The overlap between two SDs is given by the determinant of the matrix
containing all mutual orbital overlaps and thus,

$\begin{split}\braket{J^{\prime}}{I}&=\begin{vmatrix}\braket{\psi^{\prime}_{j_{\alpha}^{1}}}{\psi_{i_{\alpha}^{1}}}&\cdots&\braket{\psi^{\prime}_{j_{\alpha}^{1}}}{\psi_{i_{\alpha}^{N_{\alpha}}}}\\\
\vdots&\ddots&\vdots\\\
\braket{\psi^{\prime}_{j_{\alpha}^{N_{\alpha}}}}{\psi_{i_{\alpha}^{1}}}&\cdots&\braket{\psi^{\prime}_{j_{\alpha}^{N_{\alpha}}}}{\psi_{i_{\alpha}^{N_{\alpha}}}}\end{vmatrix}\times\begin{vmatrix}\braket{\psi^{\prime}_{j_{\beta}^{1}}}{\psi_{i_{\beta}^{1}}}&\cdots&\braket{\psi^{\prime}_{j_{\beta}^{1}}}{\psi_{i_{\beta}^{N_{\beta}}}}\\\
\vdots&\ddots&\vdots\\\
\braket{\psi^{\prime}_{j_{\beta}^{N_{\beta}}}}{\psi_{i_{\beta}^{1}}}&\cdots&\braket{\psi^{\prime}_{j_{\beta}^{N_{\beta}}}}{\psi_{i_{\beta}^{N_{\beta}}}}\end{vmatrix}\\\
&=\det\Big{(}\sum_{\nu^{\prime},\nu=1}^{N_{\alpha}}\ket{\nu^{\prime}}\bra{\nu}\braket{\psi_{j_{\alpha}^{\nu^{\prime}}}^{\prime}}{\psi_{i_{\alpha}^{\nu}}}\Big{)}\det\Big{(}\sum_{\mu^{\prime},\mu=1}^{N_{\beta}}\ket{\mu^{\prime}}\bra{\mu}\braket{\psi_{j_{\beta}^{\mu^{\prime}}}^{\prime}}{\psi_{i_{\beta}^{\mu}}}\Big{)}\\\
&=\det\Big{(}\sum_{\nu^{\prime},\nu=1}^{N_{\alpha}}\ket{\nu^{\prime}}\bra{\nu}[e^{-\hat{\kappa}}]_{j_{\alpha}^{\nu^{\prime}}i_{\alpha}^{\nu}}\Big{)}\det\Big{(}\sum_{\mu^{\prime},\mu=1}^{N_{\beta}}\ket{\mu^{\prime}}\bra{\mu}[e^{-\hat{\kappa}}]_{j_{\beta}^{\mu^{\prime}}i_{\beta}^{\mu}})\Big{)}\end{split}$
(A.11)

where we have made use of the fact that the overlap between MOs of different
spins is zero. In order to simplify the notation, we note
$\hat{\mathrm{M}}=e^{-\hat{\kappa}}$ the full $L\times L$ rotation matrix of
Eq. A.11. The size of $\hat{\mathrm{M}}$ is $L\times L$ because we employ a
spin-restricted formalism, i.e., the same matrix elements of
$\hat{\mathrm{M}}$ are used for both the $\alpha$ and $\beta$ terms. We also
introduce $\hat{\mathrm{M}}^{\sigma}_{J^{\prime}I}$ the sub-matrix of
$\hat{\mathrm{M}}$ obtained by taking its rows of indices
$J^{\prime}_{\sigma}$ and columns of indices $I_{\sigma}$. Eq. A.11 can be re-
written in matrix form as

$\braket{J^{\prime}}{I}=\det\big{(}\hat{\mathrm{M}}^{\alpha}_{J^{\prime}I}\big{)}\det\big{(}\hat{\mathrm{M}}^{\beta}_{J^{\prime}I}\big{)}$
(A.12)

and thus

$\braket{J^{\prime}}{\Psi}=\sum_{I}^{N_{\text{SD}}}c_{I}\det\big{(}\hat{\mathrm{M}}^{\alpha}_{J^{\prime}I}\big{)}\det\big{(}\hat{\mathrm{M}}^{\beta}_{J^{\prime}I}\big{)}.$
(A.13)

Finally, our goal is to optimize the entries of the $\hat{\kappa}$ to minimize

$f(\hat{\kappa})=1-|\braket{J^{\prime}}{\Psi}|^{2}.$ (A.14)

##  Appendix B Spin adaptation

Taking advantage of spin symmetry is important for both accuracy and
efficiency when the targeted molecules contain transition metals. Hence, we
employ a spin-adapted version of the DMRG algorithm. [34]. After sampling, the
resulting wavefunction is given as a superposition of CSFs. At this point we
have access to the overlap of a single CSF with the ground state which is
simply the maximum amplitude in the wavefunction. The sampled CSFs are then
transformed to SDs following the algorithm of Ref. [37]. In the following we
give further insight into our implementation of this algorithm. Let us start
by saying a few words about the indexing of the SD and CSF in the
statevectors. We call $N_{e}$ and $N_{o}$ the total number of electrons and
orbitals in the active space, respectively. First, the statevector is divided
into seniority blocks. The seniority number, $z$, gives the number of unpaired
electrons in the determinants of the corresponding block. Given a total and
projected spin, $S$ and $M$, respectively, within each seniority block there
are

$N_{\text{SD}}^{z}={z\choose{z/2+M}}$ (A.15)

possible SDs and

$N_{\text{CSF}}^{z}=\frac{2S+1}{z/2+S+1}{z\choose{z/2-S}}$ (A.16)

possible CSFs. We call $k_{z}$ the index running along the $N_{\text{SD}}^{z}$
possible SDs or, equivalently, along the $N_{\text{CSF}}^{z}$ possible CSFs.
Given the $k_{z}$ configuration, the $z$ singly occupied orbitals can be
distributed in $N_{o}\choose z$ ways among the $N_{o}$ orbitals of the active
space. The index along this axis is called $s_{z}$. Finally, given $k_{z}$ and
$s_{z}$, the paired electrons are distributed in
${N_{o}-z}\choose{\frac{N_{e}-z}{2}}$ ways among the remaining orbitals. The
index along this axis is $d_{z}$. In summary, given a seniority block $z$, a
SD or a CSF is fully determined by its indices $(k_{z},s_{z},d_{z})$.
Our algorithm starts by assigning these indices to each CSF in the DMRG
output. For each of the existing $z,s_{z},d_{z}$, all $N_{\text{SD}}^{z}$
determinants are created and their amplitude is calculated as

$A_{i}=\sum_{j=0}^{N^{z}_{\text{CSF}}}\tilde{A}_{j}D_{ij}$ (A.17)

where $i$ and $j$ replace the index $k_{z}$ used above to differentiate
between SDs and CSFs, respectively. $\tilde{A}_{j}$ is the amplitude of CSF
$j$ and $D_{ij}$ is given by the Clebsch-Gordan coefficients (see for instance
Ref. [37] for the explicit definition of $D_{ij}$). Determinants with
$|A_{i}|>t_{\text{SD}}$, where $t_{\text{SD}}$ is a threshold of choice, are
saved in the final statevector.

One subtlety arises when considering non-singlet systems. Since the spin-
adapted DMRG algorithm is better suited to singlet systems, a singlet
embedding scheme is used to treat higher spin states [34]. In this case, a
number (which depends on the molecule’s spin) of auxiliary singly occupied
orbitals need to be added to each determinant in the DMRG output such that the
overall spin becomes a singlet. The CSF to SD transformation is then performed
as explained above and the auxiliary orbitals are finally traced away in the
output statevector.

##  Appendix C Acenes

In this section we study the reliability of $p_{0}$ estimates calculated from
approximate wavefunctions. In this pursuit, we focus on acene molecules of
different sizes. The geometries of the different acenes are taken from the
work of Sharma et al. [39] and can be found in [50]. We obtain the MOs with
B3LYP/6-31G(d,p). The active space corresponds to the $\pi$ electrons
distributed in all valence $\pi$ orbitals leading to $2(2n+1)$ electrons in
$2(2n+1)$ orbitals, where $n$ is the number of rings.

We first focus on naphtalene ($n=2$). We find an approximation of its singlet
ground state with spin-adapted DMRG at bond dimension $D=500$. We make sure
that after sampling the CSFs and transformation to the SD basis the
wavefunction is normalized up to an error of $5\cdot 10^{-7}$. We diagonalize
the resulting 1-particle density matrix to obtain the NOs. The overlap matrix
between the NOs and MOs is shown on Fig. 4 (right). We can see that they are
almost identical with off-diagonal elements (2,6), (3,7), and their transpose,
slightly deviating from 0. We run the spin-adapted DMRG algorithm again with
the same bond dimension in the NOs basis. In Table 1 we show the resulting
energies, number of SDs and the weight $p_{0}$ of the determinant of highest
amplitude, i.e., its amplitude squared. Despite the high similarity between
the two orbital sets, $p_{0}$ differs in the two wavefunctions by $\sim
1.6\cdot 10^{-3}$, indicating that a change in the orbital basis directly
reflects in the value of $p_{0}$. The overlap squared between the wavefunction
converged in the MOs and NOs basis respectively is
$|\braket{\tilde{\Psi}_{0}^{\text{MO}}}{\tilde{\Psi}_{0}^{\text{NO}}}|^{2}=0.997521$
and indicates that DMRG converges to slightly different states depending on
the orbital basis. Let us write
$\ket{\tilde{\Psi}_{0}^{\text{MO}}}=\sum_{i}c_{i}^{\text{MO}}\ket{i}$ and
$\ket{\tilde{\Psi}_{0}^{\text{NO}}}=\sum_{i}c_{i}^{\text{NO}}\ket{i^{\prime}}$
using SDs $\ket{i}$ and $\ket{i^{\prime}}$. As per Eqs. 3 and 4, and
considering that the wavefunctions are normalized, we have

$|c_{i}^{\text{NO}}-\braket{i^{\prime}}{\tilde{\Psi}_{0}^{\text{MO}}}|\leq\epsilon_{r}$
(A.18)

and

$|c_{i}^{\text{MO}}-\braket{i}{\tilde{\Psi}_{0}^{\text{NO}}}|\leq\epsilon_{r}$
(A.19)

where

$\epsilon_{r}=\sqrt{2(1-\braket{\tilde{\Psi}_{0}^{\text{MO}}}{\tilde{\Psi}_{0}^{\text{NO}}})}.$
(A.20)

In Figure 4 (left), we show the errors,
$|c_{i}^{\text{NO}}-\braket{i^{\prime}}{\tilde{\Psi}_{0}^{\text{MO}}}|$ and
$|c_{i}^{\text{MO}}-\braket{i}{\tilde{\Psi}_{0}^{\text{NO}}}|$, for the first
10 determinants in $\ket{\tilde{\Psi}_{0}^{\text{NO}}}$ and
$\ket{\tilde{\Psi}_{0}^{\text{NO}}}$, respectively, together with
$\epsilon_{r}$. As anticipated in the main text, because the two orbital bases
are close, this bound is very loose. Actually, the error is much closer to
$\epsilon_{r}^{2}$ here.

Orbital basis | Method | Energy [Ha] | $N_{\text{SD}}$ | $p_{0}$
---|---|---|---|---
MOs | DMRG | -383.38958 | 1595 | 0.828714
NOs | DMRG | -383.38958 | 1578 | 0.830282
MOs | SHCI | -383.38957 | 1798 | 0.828957
NOs | SHCI | -383.38957 | 1802 | 0.830229
Table 1: Energy, number of SD and amplitude squared ($p_{0}$) of the most
important SD in the ground states converged using either the MOs or NOs as
orbital basis and DMRG or SHCI as algorithm.

Figure 4: Results obtained with naphtalene. (left) Absolute difference between
the amplitude in wavefunction A and the amplitude calculated by taking the
overlap of the corresponding SD with wavefunction B. Wavefunctions A and B
differ by the orbital basis in which they are expanded. Circles are for
wavefunction A in MOs and wavefunction B in NOs and squares for wavefunction A
in NOs and wavefunction B in MOs. The lines give the values of $\epsilon_{r}$
(and $\epsilon_{r}^{2}$), the norm of the difference between wavefunction A
and B. Blue is for DMRG and red for SHCI (the two different algorithms
employed to find the ground states). (right) Overlap matrix between the MOs
and NOs

.

We repeat the same steps for ground states converged using SHCI with no
perturbative step and a sweep epsilon of $10^{-4}$. The results are also given
in Table 1 and Figure 4 (in red). In this case
$|\braket{\psi_{0}^{\text{MO}}}{\psi_{0}^{\text{NO}}}|^{2}=0.9999991$. Here
again, $\epsilon_{r}$ is a very loose bound and the errors are much closer to
$\epsilon_{r}^{2}$.

As a second part, we shift our focus to larger acenes, namely anthracene
($n=3$), tetracene ($n=4$), pentacene ($n=5$) and hexacene ($n=6$). We first
find the ground state in the canonical MOs with spin-adapted DMRG. For
generality, we look for both the singlet and the triplet ground states. The
convergence of the energy with the bond dimension is shown in Table 2. The
selected ground state is highlighted in bold. The CSFs are sampled with a
threshold for sampling of $10^{-4}$. We transform them and keep only the SDs
above $10^{-4}$. The number $N_{\text{SD}}$ of SDs in the resulting
wavefunction together and its corresponding norm are shown in Table 3. They
form our reference states, $\ket{\tilde{\Psi}_{0}}$. In the same table, we
also give $p_{0}$(MO), the amplitude squared of the main determinant in the
$\ket{\tilde{\Psi}_{0}}$.

In a second step we obtain the NOs by diagonalizing the 1-particle density
matrix and we calculate the weight of the Aufbau filled SD is the NOs with
$\ket{\tilde{\Psi}_{0}}$. We also find a set of optimized orbitals by
maximizing the weight of a single SD with $\ket{\tilde{\Psi}_{0}}$. This
single SD is chosen to be the same as the largest amplitude one in the
reference state. We repeat the optimization for two different truncation
threshold, $N_{\text{SD}}^{\text{opt}}$, of the sum in Eq. 1: the first one
corresponds to the first 10% SDs in $\ket{\tilde{\Psi}_{0}}$ and the second
one to the first 50% SDs. Once the optimized orbitals are obtained, we
calculate the final value of $p_{0}$ between the single SD in the optimized
basis with the full $\ket{\tilde{\Psi}_{0}}$ (corresponding to $N_{\text{SD}}$
and norm of Table 3). This final value is reported in Table 3.

The discrepancies observed in the different overlaps are minor here. However,
it is interesting to see that in all cases, the optimized orbitals lead to the
largest overlap values. Moreover, and maybe more importantly, the
optimizations using 10% and 50% of SDs give the same result. This is an
encouraging sign that $p_{0}$ converges fast with $N_{\text{SD}}$. To
emphasize this point, we calculate $p_{0}$ between the SD in the optimized
basis OPT0.1 (obtained using 10% $N_{\text{SD}}$) and $\ket{\tilde{\Psi}_{0}}$
by progressively increasing $N_{\text{SD}}$ in the sum of Eq. 1. Figure 5
shows the very fast convergence of $p_{0}$(OPT0.1) with $N_{\text{SD}}$ for
all the molecules considered here. The relative error of the weight drops
below $10^{-5}$ with only 10% of the total number of SDs for all the systems
considered here.

The $p_{0}$(NO) values are unexpectedly lower than the corresponding
$p_{0}$(MO). However, as we discussed above and numerically witnessed in the
case of naphtalene, these overlaps can be tainted by a large error due to
divergence in the reference states obtained in one or the other basis and due
to truncating of the sum in Eq. 1. In fact, if we approximate the ground state
again, in the NO basis, with DMRG (same bond dimension as for the MO state),
the resulting $p_{0}$(NO) values, taken as the most important amplitude
squared in the statevector, are very close to the ones obtained with the
optimized orbitals.

In Figure 6, we compare the NOs and optimized orbitals (obtained with 10% of
the SDs) by showing their overlap matrix with the original MOs. We clearly see
that they are distinct orbital sets, the optimized ones remaining closer to
the MOs. However, it is interesting to remark that the two types of orbitals
also hold similarities, i.e., we observe a common structure in their overlap
matrix with the original MOs.

anthracene (14,14) | tetracene (18,18) | pentacene (22,22) | hexacene (26,26)
---|---|---|---
D=500 -536.049037 | D=500 -688.687930 | D=1000 -841.361452 | D=1500 -993.998918
D=1000 -536.049038 | D=1000 -688.687953 | D=1500 -841.361630 | D=2000 -993.999350
| D=1500 -688.687954 | D=2000 -841.361700 | D=2500 -993.999574
D=500 -535.975915 | D=1000 -688.622454 | D=1500 -841.320648 | D=2000 -993.959774
D=1000 -535.975915 | D=1500 -688.622460 | D=2000 -841.320829 | D=2500 -993.960074
Table 2: Convergence of the energy of different acenes, with active space (number of electrons,number of orbitals), obtained with DMRG at bond dimension $D$. The upper and lower rows are for singlet and triplet states, respectively. Molecule | S | Energy [Ha] | $N_{\text{SD}}$ | Norm |  $p_{0}$(MO) |  $p_{0}$(NO) |  $p_{0}$(OPT0.1) |  $p_{0}$(OPT0.5)
---|---|---|---|---|---|---|---|---
anthracene | 0 | -536.049037 | 7522 | 0.99995 | 0.81655 | 0.81481 | 0.81892 | 0.81892
tetracene | 0 | -688.687953 | 13029 | 0.99909 | 0.76942 | 0.76833 | 0.77301 | 0.77301
pentacene | 0 | -841.361700 | 50842 | 0.99506 | 0.67247 | 0.66881 | 0.67774 | 0.67774
hexacene | 0 | -993.999574 | 64098 | 0.99207 | 0.61155 | 0.60896 | 0.61800 | 0.61800
anthracene | 1 | -535.975916 | 4230 | 0.99959 | 0.78481 | 0.76678 | 0.79843 | 0.79843
tetracene | 1 | -688.622454 | 10208 | 0.99847 | 0.73751 | 0.72471 | 0.74703 | 0.74703
pentacene | 1 | -841.320814 | 29829 | 0.99308 | 0.67598 | 0.66598 | 0.68507 | 0.68506
hexacene | 1 | -993.960030 | 38755 | 0.98923 | 0.62550 | 0.61755 | 0.63316 | 0.63316
Table 3: Results for the four studies acenes in both singlet and triplet
states. The ground state is approximated with spin-adapted DMRG and leads the
reported energies. After sampling and transforming to the SD basis, the
statevectors contain $N_{\text{SD}}$ SDs and have the reported norm.
$p_{0}$(MO) is the amplitude squared of the most important determinant.
$p_{0}$(NO) is the weight of the Aufbau filled determinant in the NO basis
with the statevector. $p_{0}$(OPT$X$) is the weight of the most important
determinant in the reference statevector expressed in the optimized orbital
basis obtained by using the first $XN_{\text{SD}}$ SDs. Figure 5: Convergence
of $p_{0}$ with the number of SDs included in the reference statevector,
$\ket{\tilde{\Psi}_{0}}$. The full lines correspond to the error between
$p_{0}$ calculated using $N_{\text{SD}}$ with respect to the value obtained
from all sampled SDs. The error bound, $\epsilon_{r}$, is also shown in dashed
lines. The results are given for acenes of different sizes $n$ and for both
singlet and triplet states. Figure 6: Overlap matrix between the starting MOs
and either the NOs or the optimized orbitals (OPT) obtained with the first 10%
of the SDs in the reference state for each acenes in both singlet and triplet
states.

##  Appendix D Iron-Sulfur clusters

One of our focuses is on determining whether optimizing the orbital basis can
significantly enhance the weak overlaps observed in Ref. [6]. To maintain a
seamless progression in our research, we diligently replicate the findings of
Ref. [6] and build upon them. Here we focus on the clusters with a number of
iron atoms $N_{\text{Fe}}=2$ and $N_{\text{Fe}}=4$. Tab. 4 provides details
about the studied iron-sulfur clusters, supplying all the essential
information for our investigation.

In the following paragraph, we describe the procedure followed to obtain the
active space, which is the same as the one of Ref. [6]. We start by running a
high spin unrestricted density functional theory (DFT) calculation with PySCF
[51]. The specific spin numbers, as well as information about basis sets and
functionals, can be found in Tab. 4. For all systems we use the spin-free
exact two-component (sf-X2C) Hamiltonian.

molecule | oxidation | charge | high S | low S | AS | functional | basis set
---|---|---|---|---|---|---|---
[Fe2S2]-3 [52] | Fe${}^{\text{II}}$, Fe${}^{\text{III}}$ | -3 | 9/2 | 1/2 | (31e,20o) | BP86 | TZP-DKH
[Fe2S2]-2 [52] | Fe${}^{\text{III}}$, Fe${}^{\text{III}}$ | -2 | 5 | 0 | (30e,20o) | BP86 | TZP-DKH
[Fe4S4]-2 [52] | 2Fe${}^{\text{II}}$, 2Fe${}^{\text{III}}$ | -2 | 9 | 0 | (54e,36o) | BP86 | TZP-DKH
[Fe4S4] | 4Fe${}^{\text{III}}$ | 0 | 10 | 0 | (52e,36o) | BP86 | TZP-DKH
Table 4: Details about the studied iron-sulfur clusters

The spin averaged density matrix, $\rho=1/2(\rho_{\alpha}+\rho_{\beta})$ is
calculated in the basis of the resulting Kohn-Sham (KS) orbitals and
transformed to the Lowdin basis, $\rho^{\prime}=S^{1/2}\rho S^{1/2}$. The NOs
and corresponding occupation numbers are obtained by diagonalizing
$\rho^{\prime}$. The NOs are transformed to the non-orthogonal AO basis and
classified in a core, active, or virtual group of orbital according to their
occupation number, i.e., orbitals with occupation number above 1.95, between
0.05 and 1.95 and below 0.05 are classified as core, active and virtual
respectively. The core and active orbitals are then split localized with the
Pipek-Mezey algorithm. The active spaces are then constructed as in Ref. [6],
from the full valence space of the Fe 3d, S 3p, and bonding ligand orbitals
around each Fe atom. These active orbitals are localized one more time and
shown for [Fe2S2]-2 in Figure 7. The orbital pictures for all systems are
available in [50]. The resulting active space sizes are given in Tab. 4.

Figure 7: Pictures of the orbitals in the active space of [Fe2S2]-2.

We approximate the ground states of the four clusters with spin-adapted DMRG
at bond dimension $D=8000$. We employ a singlet embedding [34] scheme for
[Fe2S2]-3. The resulting energies, given in Table 5, are in close agreement
with Ref. [6]. We sample the CSFs with $t_{\text{CSF}}=10^{-4}$ for
$N_{\text{Fe}}=2$ and $t_{\text{CSF}}=5\cdot 10^{-4}$ for $N_{\text{Fe}}=4$.
We transform them and discard all SDs with absolute amplitude below $10^{-4}$.
The resulting number of SDs in the statevector, $\ket{\tilde{\Psi}_{0}}$, and
its norm are summarized in Table 5. In all cases, we sort these reference
statevectors by decreasing order in absolute amplitude. We optimize the
orbitals using the first 10% SDs and the final value we report in Table 5 is
$p_{0}(\text{OPT})=|\braket{J^{\prime}}{\tilde{\Psi}_{0}}|^{2}$ where
$\ket{J^{\prime}}$ is expressed in the optimized basis and
$\ket{\tilde{\Psi}_{0}}$ (expressed in the original MOs) contains the full set
of sampled SDs corresponding to $N_{\text{SD}}$ of Table 5. The orbital
optimization increases the weight by more than one order of magnitude in all
cases. Moreover, for all these clusters, optimizing the orbitals allows to
reach a higher $p_{0}$ than starting with a CSF. As a matter of comparison,
the CSF of largest amplitude in [Fe4S4]-2 comprises more than 48000 SDs.

Molecule | Energy [Ha] | Norm | $N_{\text{SD}}$ |  $p^{\text{CSF}}_{0}$(MO) |  $p_{0}$(MO) |  $p_{0}$(OPT) |  $p_{0}$(NO)
---|---|---|---|---|---|---|---
[Fe2S2]-2 | -5092.8710 | 9.85e-1 | 3121044 | 3.21e-2 | 5.35e-3 | 4.53e-2 | 3.36e-2
[Fe2S2]-3 | -5092.7187 | 9.94e-1 | 734408 | 3.19e-2 | 1.76e-2 | 1.80e-1 | 2.46e-2
[Fe4S4]-2 | -8432.8070 | 3.65e-1 | 2146432 | 2.94e-4 | 2.94e-5 | 2.67e-3 | -
[Fe4S4] | -8432.6267 | 3.65e-1 | 2127398 | 4.37e-3 | 3.97e-4 | 5.16e-3 | -
Table 5: Energy, norm characteristics and overlaps of the iron-sulfur
clusters. The ground state is approximated with spin-adapted DMRG ($D=8000$)
and leads to the reported energies. After sampling and transforming to the SD
basis, the statevectors contain $N_{\text{SD}}$ SDs and have the reported
norm. $p^{\text{CSF}}_{0}$(MO) and $p_{0}$(MO) are the amplitude squared of
the most important CSF and SD, respectively. $p_{0}$(OPT) is the overlap of
the most important SD in the reference statevector expressed in the optimized
orbital basis obtained by using the first 10% SDs. $p_{0}$(NO) is the
amplitude squared of the most important determinant in the state resulting
from running DMRG in the NO basis (see main text).

As for the acenes, we study the convergence of the weight with
$N_{\text{SD}}$. The results are shown in Figure 8. We see again that
$p_{0}$(OPT) converges very quickly with $N_{\text{SD}}$. We also study the
effect of using an approximate wavefunction to optimize the orbitals. To this
aim, we obtain additional states, $\ket{\tilde{\Psi}_{0}(D)}$ with spin-
adapted DMRG at bond dimensions $D\in\\{200,400,600,800,1000\\}$. We find an
optimized orbital basis, OPT$D$, using each of the
$\ket{\tilde{\Psi}_{0}(D)}$. In the top pannels of Figure 9 we show $\Delta
E=|E(8000)-E(D)|$, the difference between the energies of
$\ket{\tilde{\Psi}_{0}(8000)}$ and $\ket{\tilde{\Psi}_{0}(D)}$ as well as the
most important amplitude squared in $\ket{\tilde{\Psi}_{0}(D)}$. After the
optimization, we calculate $|\braket{J^{\prime}}{\tilde{\Psi}_{0}(8000)}|^{2}$
with $\ket{J^{\prime}}$ in OPT$D$ and show the resulting values in the bottom
pannels of Figure 9. There we also display $p_{0}$(OPT) of Table 5 for
reference. Interestingly, these results clearly show that statevectors
obtained at a smaller $D$, while inaccurate in energy, already allow to find
an orbital basis which substantially improves the overlap.

Moreover, we compare the NO and optimized bases. The NOs are obtained by
diagonalizing the 1-particle density matrix of the $D=8000$ state. In Figure
10, we show the overlap matrix between the starting MOs and the NOs as well as
between the MOs and the optimized orbitals. As we observed in the acene case,
the NOs and optimized orbitals are different and the latter are much closer to
the initial MOs.

Here, the direct calculation of $p_{0}$(NO) between the Aufbau filled SD in
the NO basis, and the MO reference state leads to values much lower than
$p_{0}$(MO). However, the density of non-diagonal elements in the MO-NO
overlap matrix is large. This translates to having more non-vanishing terms in
the sum of Eq. 1. We therefore expect the errors defined in Eqs. 4 and 5 to be
large, hence resulting in unreliable $p_{0}$ estimates. To allow for a more
fair comparison between the optimized and the NO basis, we run DMRG again in
the NO basis and get $p_{0}$(NO) as the amplitude squared of the most
important determinant in the resulting state. This is only done for the two
$N_{\text{Fe}}=2$ clusters, as we can afford a large enough bond dimension to
compensate for the loss of locality in the orbitals. As the energy and overlap
converge at $D=500$, we use this bond dimension to obtain our states in the NO
basis. The results are given in Table 5. In both cases we obtain better
weights in the optimized basis.

Finally, note that another way of verifying the accuracy of the calculated
$p_{0}$ is to rotate the MPS underlying $\ket{\tilde{\Psi}_{0}}$ to the
optimized basis and sample from it. However, this procedure was too costly for
the systems considered here as $D$ blows up due to the loss of locality in the
optimized orbitals. This explains our choice of demonstrating the convergence
of $p_{0}$ with $N_{\text{SD}}$ and $D$ instead.

Figure 8: Convergence of $p_{0}$(OPT) with the number of SDs included in the
reference statevector, $\ket{\tilde{\Psi}_{0}}$. The full line corresponds to
the error between $p_{0}$(OPT) calculated using $N_{\text{SD}}$ with respect
to the value of the weight obtained from all sampled SDs. The error bound
$\epsilon_{r}$ is shown in dashed line. Plots a, b, c and d are for [Fe2S2]-2,
[Fe2S2]-3, [Fe4S4]-2 and [Fe4S4], respectively. Figure 9: (Top) Convergence of
the energy and weight with the DMRG bond dimension, $D$. $\Delta
E=|E(8000)-E(D)|$ is the difference between the reference energy and the one
obtained at reduced bond dimension $D$. Circles represent the most important
absolute amplitude squared in $\ket{\Psi(D)}$. (Bottom) Circles correspond to
$|\braket{J^{\prime}}{\Psi(8000)}|^{2}$ where $\ket{J^{\prime}}$ is in the
orbital basis optimized using $\ket{\Psi(D)}$. The dashed line is
$p_{0}$(OPT), i.e., $|\braket{J^{\prime}}{\Psi(8000)}|^{2}$ where
$\ket{J^{\prime}}$ is in the orbital basis optimized using $\ket{\Psi(8000)}$.
Plots a, b, c and d are for [Fe2S2]-2, [Fe2S2]-3, [Fe4S4]-2 and [Fe4S4],
respectively. Figure 10: Overlap matrix between the starting MOs and either
the NOs or the optimized orbitals (OPT). Plots a, b, c and d are for
[Fe2S2]-2, [Fe2S2]-3, [Fe4S4]-2 and [Fe4S4], respectively.

##  Appendix E Cytochrome P450

Figure 11: Oxo model of the cytochrome P450 enzyme.

The cytochrome P450 enzyme holds significant importance in the pharmaceutical
industry. Recent studies showcased its complex electronic structure and
highlight the system as an ideal use case for the first generation of quantum
computers [41]. The numerical results presented in this paper also underscore
that a single determinant already exhibits a considerable overlap with the
ground state of P450. However, further enhancing this overlap can reduce the
runtime of quantum phase estimation, making the computation more feasible for
industrial applications. To evaluate our method’s performance on the P450
system, we utilize the O2-bound model (oxo, see Fig. 11) and active spaces
available in [53]. We approximate the ground state with spin-adapted DMRG at
bond dimension $D=1500$ for the various spins. In all cases we use a singlet-
embedding. For both the “X” (47e,43o) and “G” (63e,58o) active spaces, as well
as for the $S=1/2$ and $S=5/2$ states, our method yields an approximate 10%
improvement in $p_{0}$, see Tab. 6. We note that the overlaps reported in Fig.
9 of Ref. [41] are considerably higher as the authors report the highest
amplitude of a CSF, while we report the weight (amplitude squared) of a single
SD.

Active space | S | Energy [Ha] |  $p_{0}$(MO) |  $p_{0}$(OPT)
---|---|---|---|---
G | $5/2$ | -2756.7356 | 0.53 | 0.58
G | $1/2$ | -2756.8039 | 0.11 | 0.12
X | $5/2$ | -2756.9426 | 0.48 | 0.53
X | $1/2$ | -2757.0156 | 0.093 | 0.11
Table 6: Results for the cytochrome P450 enzyme. The active spaces G and X are
taken from Ref. [41] and correspond to (47e,43o) and (63e,58o), respectively.
The ground state of all systems are approximated with spin-adapted DMRG
($D=1500$).