Stabilizer Approximation II: From H$_2$O To C$_6$H$_6$

TL;DR

This paper demonstrates that stabilizer states can effectively approximate the ground states of complex molecules like water and benzene, even under strong distortions, suggesting their usefulness in quantum chemical simulations.

Contribution

The study extends the stabilizer method to larger, more complex molecules, showing their effectiveness as reference states in strongly correlated systems.

Findings

Stabilizer states approximate ground states well at various bond lengths.

Effective for molecules with strong static correlation.

Potential as initial states in quantum chemical simulations.

Abstract

We apply the stabilizer method to the study of some complicated molecules, such as water and benzene. In the minimal STO-3G basis, the former requires 14 qubits, and the latter 72 qubits, which is very challenging. Quite remarkably, We are still able to find the best stabilizer states at all the bond lengths. Just as the previously studied H$_2$, LiH and BeH$_2$ molecules, here the stabilizer states also approximate the true ground states very well, especially when the molecules are strongly distorted. These results suggest stabilizer states could serve as natural reference states when the system involves strong static correlation. And in the language of quantum computing, one would expect stabilizer states to be natural initial states for chemical simulation.

Equations

$$H_0$$

$$\bar{Z}\bar{I}II,\bar{I}\bar{I}ZI,-\bar{I}\bar{Z}II,-\bar{I}\bar{I}IZ.$$

$$H_1$$

$$-\bar{Z}\bar{I}ZI,-\bar{Z}\bar{I}IZ,-\bar{I}\bar{Z}ZI.$$

$$\frac{|11;00⟩\pm |00;11⟩}{\sqrt{2}}.$$

$$|1^2;0^2⟩,|0^2;1^2⟩.$$

$$-\bar{Z}\bar{I}ZI,-\bar{I}\bar{Z}IZ,-\bar{X}\bar{X}YY,-\bar{I}\bar{I}ZZ.$$

$$\frac{|10;01⟩-|01;10⟩}{\sqrt{2}}.$$

$$|0^2{1}^5;0^2{1}^5⟩.$$

$$|1^4;0^4⟩,|0^4;1^4⟩.$$

$$-\bar{I}\bar{I}\bar{I}\bar{Z}IIIZ,-\bar{I}\bar{I}\bar{Z}\bar{I}IIZI,-\bar{I}\bar{Z}\bar{I}\bar{I}IZII,-\bar{Z}\bar{I}\bar{I}\bar{I}ZIII,$$

$$-\bar{I}\bar{I}\bar{X}\bar{X}IIYY,-\bar{I}\bar{X}\bar{X}\bar{I}IYYI,-\bar{X}\bar{X}\bar{I}\bar{I}YYII.$$

$$(+|1000;0111⟩+|1011;0100⟩-|0100;1011⟩-|0111;1000⟩$$

$$\text{leavevmode}+|0010;1101⟩-|0001;1110⟩-|1101;0010⟩+|1110;0001⟩)/2\sqrt{2}.$$

$$(+|1010;0101⟩-|0110;1001⟩-|1100;0011⟩-|0000;1111⟩$$

$$\text{leavevmode}-|1001;0110⟩+|0101;1010⟩-|0011;1100⟩-|1111;0000⟩)/2\sqrt{2}.$$

$$|0^{15}{1}^{21};0^{15}{1}^{21}⟩.$$

$$|1^{6};0^{6}\rangle,\quad|0^{6};1^{6}\rangle.\leavevmode\nobreak\$$

$$-\bar{Z}\bar{I}\bar{I}\bar{I}\bar{I}\bar{I}ZIIIII,$$

$$-\bar{I}\bar{Z}\bar{I}\bar{I}\bar{I}\bar{I}IZIIII,$$

$$-\bar{I}\bar{I}\bar{Z}\bar{I}\bar{I}\bar{I}IIZIII,$$

$$-\bar{I}\bar{I}\bar{I}\bar{Z}\bar{I}\bar{I}IIIZII,$$

$$-\bar{I}\bar{I}\bar{I}\bar{I}\bar{Z}\bar{I}IIIIZI,$$

$$-\bar{I}\bar{I}\bar{I}\bar{I}\bar{I}\bar{Z}IIIIIZ,$$

$$-\bar{X}\bar{X}\bar{I}\bar{I}\bar{I}\bar{I}YYIIII,$$

$$-\bar{I}\bar{X}\bar{X}\bar{I}\bar{I}\bar{I}IYYIII,$$

$$-\bar{I}\bar{I}\bar{X}\bar{X}\bar{I}\bar{I}IIYYII,$$

$$-\bar{I}\bar{I}\bar{I}\bar{X}\bar{X}\bar{I}IIIYYI,$$

$$-\bar{I}\bar{I}\bar{I}\bar{I}\bar{X}\bar{X}IIIIYY.$$

$$|1^{18};0^{18}⟩,|0^{18};1^{18}⟩$$

$$-ZIIIIIZIIIII,$$

$$-IZIIIIIZIIII,$$

$$-IIZIIIIIZIII,$$

$$-IIIZIIIIIZII,$$

$$-IIIIZIIIIIZI,$$

$$-IIIIIZIIIIIZ,$$

$$+XXIIIIXXIIII,$$

$$+IXXIIIIXXIII,$$

$$+IIXXIIIIXXII,$$

$$+IIIXXIIIIXXI,$$

$$+IIIIXXIIIIXX.$$

Full text

1. I INTRODUCTION

2. II THE CALCULATION

3. II.1 H2

4. II.2 H2O

5. II.3 C6H6

6. II.3.1 Hückel case

7. II.3.2 Overall case

8. III DISCUSSION

I INTRODUCTION

Modern computational chemistry crucially depends on the Hartree-Fock (HF) approach Hartree ; Slater ; Fock . When the HF wavefunction provides a good approximation to the system, one could develop various perturbative methods to improve the accuracy to a rather high level. However, when the HF description deviates much from the true state, an accurate description of the system is usually difficult to achieve. For example, when a large molecule is distorted from the equilibrium configuration, many orbitals could become degenerate. The occupation of electrons in these degenerate orbitals will no longer be well described by a single HF state, but many HF states at the same time. To figure out these reference states in a huge Hilbert space will be a challenging task. Such a problem is recently illustrated in TE22 .

Something similar happens in the recent application of quantum computing to chemical simulations. When applying variational algorithms to study the ground state of molecules, one usually takes the HF states, or some other product states, as the initial states. Around the equilibrium point this leads to quite nice results, just as in the classical treatment. But when the molecules are distorted strongly, the results start to deviate HW22 . The reason is also similar to the classical case, though in a different language: the molecules could develop huge quantum entanglement in the dissociation procedure, which is not guaranteed in the shallow circuits commonly used.

In order to solve this problem, the key step is to prepare a nice reference/inital state with suitable entanglement. Such a state could be a nearly equal-weighted supposition of many product states, thus may be well approximated by stabilizers states Gottesman96 ; Gottesman97 . According to the Gottesman-Knill theorem, such states can be efficiently prepared with Clifford circuits Gottesman98 , if they are known in advance. However, if the required states are unknown, directly searching them using Clifford circuits could be very inefficient Shang2021 ; CAFQA ; CAFQA2 . An efficient way to find out the required stabilizer states, though heuristic, is to determine the stabilizers directly from the Pauli terms in the Hamiltonian. This is used long ago in the development of the toric code, in which the stabilizers encode the Hamiltonian and the exact ground state by construction Kitaev97 .

Recently we applied this approach to the chemical Hamiltonians, and found the required stabilizer states very quickly SA-I for the small molecules such as H2, LiH and BeH2. In this paper we try to extend this approach to more complicated molecules, to see if the stabilizer states could still be efficiently obtained. According to the study in CAFQA , the Clifford-circuit search is already difficult for the water molecule, and some ambiguous results have been reported. As for benzene, a complete description is beyond the ability of current quantum hardware and simulators, so the frozen-core approximation is usually taken TE22 . Therefore, we use these two cases to test the capability and efficiency of the stabilizer method.

II THE CALCULATION

II.1 H2

First we repeat the analyses for the hydrogen molecule SA-I to illustrate the procedure, but choose the Jordan-Wigner (JW) transformation instead of the parity one Parity . At equilibrium point $d=0.74$Å, the qubit-form Hamiltonian reads:

[TABLE]

We use $X,Y,Z,I$ for the gates acting on the spin-down orbitals, and $\bar{X},\bar{Y},\bar{Z},\bar{I}$ for the gates on the spin-up ones. To approximate the ground state, one must choose the stabilizers in such a way that the corresponding energy is as low as possible. The numerical coefficients of the Pauli terms above suggest the following choice:

[TABLE]

The corresponding state is the Hatree-Fock state $|01;01\rangle$, and the corresponding energy can be obtained by inserting the stabilizing condition back into the Hamiltonian.

When the molecule is stretched, the coefficients in the Hamiltonian change. For example, at $d=2.8$Å, the Hamiltonian becomes:

[TABLE]

An intuitive choice for the stabilizers is

[TABLE]

When these are selected, the Hamiltonian could be completely simplified. So the stabilizer states are doubly degenerate. To fix the states explicitly, one may choose one more stabilizer as either $\bar{X}\bar{X}YY$ or $-\bar{X}\bar{X}YY$. Then the correspond states read

[TABLE]

These are variants of the 4-qubit Greenberger-Horne-Zeilinger (GHZ) state. We will write $n$ repeating bits as $1^{n}$ or $0^{n}$ for short. Actually the degeneracy destroys the entanglement, and one could represents them as product states. This can be done by choosing the last stabilizer as $\pm\bar{I}\bar{I}IZ$, and get:

[TABLE]

Now these states are of spin 1, and different resonating forms of them give (5). The appearance of such high-spin states should be expected, according to the Hund’s rule Hund .

A careful analysis of the above Hamiltonian reveals that there is another degenerate state, with the stabilizers:

[TABLE]

And the state reads:

[TABLE]

The reason why such a state could achieve an energy as low as (6) could be explained by the interference effect of the superposition, or resonating Pauling . As the resonating is between two different spin sectors, we call it a “spin-resonating state". This also explains why the other combination, with a plus sign inbetween, has a higher energy. Notice that for all the states (5,6,8), the two spatial orbitals are half-filled, in agreement with the asymptotic behavior of full configuration interaction (FCI) results HJO .

Repeating this procedure for all the bond lengths, we obtain the approximate potential energy curve. The results are plotted in Fig. 1, together with those from HF and FCI. As emphasized in CAFQA and SA-I , the stabilizer method behaves correctly in the dissociation limit, and thus recover almost the whole correlation energy missed in the HF approximation.

II.2 H2O

For H2O there are 7 atomic orbitals involved, and we need 14 qubits to describe them. After JW transformation the Hamiltonian contains about 3000 Pauli terms. More than half of these terms have coefficients smaller than $10^{-6}$ Hartree, and are then neglected in previous studies. We focus on the dissociation with the H-O-H bond angle fixed at equilibrium value $105^{\circ}$.

Around the equilibrium point, the stabilizer method reproduces the HF state, in which 5 spatial orbitals are fully occupied, and the other 2 are empty. Explicitly, the state is

[TABLE]

When the O-H bond is stretched, the lowest 3 spatial orbitals are unaffected. For simplicity, We will omit them in the following expressions. The higher 4 orbitals start to change when the O-H distance goes beyond about 1.5Å. Same as the hydrogen molecule case, they become half-filled eventually. This is also in agreement with the FCI results HJO . Again we have different spin states. First we will immediately get the high spin states:

[TABLE]

Further analysis of the Hamiltonian shows that there are also spin-resonating states, determined by the stabilizers:

[TABLE]

Now we have 7 stabilizers for 8 qubits, which means a doubly degenerate state. The degeneracy could be lifted by adding an extra stabilizer, which we choose to be $\pm\bar{I}\bar{I}\bar{I}\bar{I}ZZZZ$. With $-\bar{I}\bar{I}\bar{I}\bar{I}ZZZZ$ included, we obtain an odd-spin state:

[TABLE]

If choosing $\bar{I}\bar{I}\bar{I}\bar{I}ZZZZ$ instead, we get the following even-spin state:

[TABLE]

One should notice that these states do not have conserved numbers of spin-up and spin-down electrons. However, since (13) is degenerate with the high-spin ones (10), we could make a new combination of them to recover the conservation of spin-up and spin-down electrons. The resulting state will have total spin 0. In practice, these low-spin stabilizer states (12,13) are not easy to identify, in contrast to the high-spin ones (10). This is similar to the situation in the Clifford-circuit approach, in which the lowest singlet state for H2O is never found actually CAFQA .

The corresponding dissociation curve is plotted in Fig. 2, together with the HF and FCI results. The correct dissociation limit is clearly recovered again. Moreover, the stabilizer method gives a single degenerate curve, in contrast to those from the Clifford-circuit search CAFQA . In HW22 it is pointed out that variational algorithms with initial product states fail to give accurate results for stretched molecules, including H2O. Fig. 2 shows that stabilizer states may help improving the accuracy in such cases, if taken as initial states.

II.3 C6H6

Now we turn to a challenging system, the benzene molecule. It involves 42 electrons and 36 spatial orbitals, so we need 72 qubits to describe them. And the Hamiltonian contains about 2.5 million Pauli terms. This makes it very difficult to treat the whole system directly. Previous study utilizes the frozen-core approximation, and chooses active spaces with 8-16 qubits TE22 . This is in some sense comparable to the traditional Hückel approximation Huckel . Here we try to go beyond such an approximation, utilizing the stabilizer method.

II.3.1 Hückel case

First we consider the situation when only the C-C bonds are stretched, keeping the hexagonal symmetry. This is recently studied in detail in TE22 . Due its similarity with the Hückel picture Huckel , we call it Hückel case.

At equilibrium, we expect 21 of the total 36 spatial orbitals would be fully occupied, and get the HF state:

[TABLE]

When the C-C bonds are stretched, the 6 active electrons in the delocalized $\pi$ orbitals will first be excited. The other 36 electrons may stay paired up in the remaining orbitals. Thus we should get a state with 18 fully-occupied orbitals, 6 half-filled ones, and 12 empty ones. This is indeed the case. For the half-filled orbitals, the 6 active electrons could take the same spin orientation and give rise to the high-spin state:

[TABLE]

Alternatively, they could hop between orbitals with different spin orientation, leading to the spin-resonating state. This is described by the following stabilizers:

[TABLE]

The corresponding state is a little complicated, and will not be shown explicitly here. It is doubly degenerate, and also degenerate with the high-spin states (15).

The appearance of 6 unpaired electrons has been numerically indicated in recent study TE22 . In our calculation, we find another stabilizer state with even lower energy, at least in the asymptotic region. This state has 12 fully-occupied orbitals, 18 half-filled ones, and 6 empty ones. We will investigate the details of this state later. Therefore, the frozen-core approximation fails to give the right physical picture here. This is not surprising: the resulting CH radicals are not stable, and contain active electrons. The comparison of the energy profiles calculated from different methods are shown in Fig. 3. Since FCI is not available anymore, we provide the results calculated with the Coupled-Cluster Singles-and-Doubles (CCSD) method. As a whole, our results from the stabilizer approximation resemble and improve the “full diagonalization of the active-space Hamiltonian" results in TE22 .

We also notice that the stabilizer curves do not approach the dissociation limit accurately. This is reasonable, since the HF/stabilizer description of the CH radical at fixed bond length is not accurate enough.

II.3.2 Overall case

To probe the state with 18 active electrons, we dissociate the benzene completely into individual atoms. This is done by stretching the C-C bonds and C-H bonds simultaneously, keeping their ratios fixed.

Again, the 18 electrons could occupy spin orbitals of the same orientation and lead to the high-spin states:

[TABLE]

Among them, the 6 electrons from the delocalized $\pi$ orbital could further hop between orbitals with different spin orientations, and gives the spin-resonating state (16). As for the remaining 12 unpaired electrons, it would be convenient to consider them from the atomic orbital point of view. Then one could guess that 6 electrons would occupy the $1s$ orbitals of the H atoms, and the other 6 occupy the $2p$ orbitals of the C atoms. Since there are 6 $2p_{x}$ orbitals and 6 $2p_{y}$ orbitals, the 6 electrons may hop between them. Indeed, for these 6 electrons we get the spatial-resonating state:

[TABLE]

The explicit form of the stabilizers differ sightly from those of (16). All the 12 electrons still remain in the same spin sector, as required by Hund’s rule. One may expect some fully-entangled state of the total 18 electrons, but we do not find such states yet.

It turns out the high-spin states (17), the one with spin resonating(16), and the one with both spin resonating(16) and spatial resonating (18) all have nearly the same energy in the asymptotic region. In Fig. 4 we plot the dissociation curves calculated from different methods. For the results from the stabilizer method, all the mentioned states fall on a single curve. remarkably, the stabilizer method nearly reproduces the asymptotic result. The deviation from the asymptotic line is about $0.12$ Hartree when the C-C distance is $5.0$Å. We are unable to obtain the CCSD results in the large distance region due to appearance of numerical singularities.

III DISCUSSION

We have demonstrated the efficiency and efficacy of the stabilizer method with two examples, the water and benzene molecules. Even for the challenging benzene molecule, we are able to find the proper stabilizers quickly for all the molecular configuration. In contrast to the HF method, the stabilizer states always stay close to the true ground states. Therefore, they could serve as proper reference/initial states for further refined calculations, either for classical methods or quantum simulation. In particular, it would be interesting to check if they qualify as suitable initial states for the quantum phase estimation algorithm, to guarantee the exponential quantum advantage in quantum chemistry GKC2208 .

To confirm the results obtained in this paper, we need to apply the method to even more complicated systems, for example, to molecules containing transition metal atoms. Recently the chromium dimer has already been studied with Clifford-circuit search CAFQA . It would be interesting to investigate such system directly with stabilizers.

States from the JW transformation have a direct interpretation on the chemistry side. In studying the water molecule we find that the stabilizer states do not preserve electrons with individual spin orientations. Then what would these states present themselves with the parity and Bravyi-Kitaev transformation BK ? Some preliminary study has been performed recently with Clifford circuits CAFQA2 . The stabilizer analysis may provide alternative answers.

Many combinational optimization problems can be reduced to a groundstate-finding problem for a classical Hamiltonian. In CAFQA2 it is concluded that Clifford circuits may not help solving such problems. However, from the stabilizer point of view, these classical Hamiltonians should be even easier to solve. Moreover, stabilizer states could provide the exact solution to such classical problems, not just an approximation.

We leave all these for future study.

Framework and Resource

The Hartree-Fock, Coupled-Cluster Singles-and-Doubles, and Full Configuration Interaction calculation are performed using Python-based Simulations of Chemistry Framework (PySCF) PySCF . We develop original submodules to implement the Jordan-Wigner transformation and the stabilizer-finding process.

All the calculation has been done using a desktop computer with a 8-core CPU and 64 GB RAM. With the stabilizer method, a single-point calculation of the benzene molecule can be finished within 10 minutes.