Enhancing the electron pair approximation with measurements on trapped-ion quantum computers

The oo-upCCD AnsatzThe unitary pair coupled-cluster double (upCCD) ansatz is$$\vert {\Psi }_{{\rm{u}}{\rm{pCCD}}}\rangle ={e}^{T-{T}^{\dagger }}\left\vert {\rm{HF}}\right\rangle$$
(1)
in which T is the pair-double cluster operator, defined as$$T=\sum _{ia}{t}_{i}^{a}{a}_{a\alpha }^{\dagger }{a}_{a\beta }^{\dagger }{a}_{i\beta }{a}_{i\alpha }$$
(2)
in which i and a are indices for occupied and unoccupied orbitals in the HF state. \({a}_{p\alpha }^{\dagger }\) (\({a}_{p\beta }^{\dagger }\)) and apα (apβ) are the fermionic creation and annihilation operators in the p-th spin-up (α) or spin-down (β) orbital.The exponentiation of the pair-excitation operator can be efficiently implemented with the following circuit, given in Fig. 1. Once the circuit is defined, one needs to measure the energy expectation value \(\langle{\Psi }_{{\rm{upCCD}}}| H| {\Psi }_{{\rm{upCCD}}}\rangle\) for the second-quantized Hamiltonian H. In the full electronic Hamiltonian, H contains O(N4) terms, with N being the number of qubits. Yet, owing to pair symmetry, most of these terms do not contribute to the energy within the upCCD approximation. After removing these pair-breaking terms, one finds that only three measurements are needed (in the Pauli X, Y, and Z basis, respectively) to compute the energy. Moreover, the number of basis measurements is independent of problem size, and three circuits are all that is needed for any upCCD energy evaluation.Fig. 1: A quantum circuit that implements the Givens rotation.S is the phase gate. H is the Hadamard gate. \(R_{y}(\theta)\) is the single qubit Y-rotation gate by angle \(\theta\). Entangling gates are \(CNOT\) gates.One potential concern is that the upCCD ansatz defined in Equation (1) is not invariant to the choice of underlying orbitals. Previous studies28,30,32,33 on similar wave functions have found that it is necessary to optimize the orbitals along with the cluster amplitudes, especially for strongly correlated systems. The orbital optimized upCCD ansatz is$$\vert{\Psi }_{{\rm{oo}}-{\rm{upCCD}}}\rangle ={e}^{K}{e}^{T-{T}^{\dagger }}\left\vert {\rm{HF}}\right\rangle$$
(3)
in which there are two different sets of parameters: 1) circuit parameters in the cluster operator T, and, 2) orbital rotation parameters in the the orbital rotation operator K, which is defined as$$K=\sum _{p > q}\sum _{\sigma }{K}_{pq}({a}_{p\sigma }^{\dagger }{a}_{q\sigma }-{a}_{q\sigma }^{\dagger }{a}_{p\sigma })$$
(4)
Here K is an anti-Hermitian matrix and σ indexes the spin.The VQE perturbation theoryAs mentioned previously, the cluster operator in upCCD does not contain any broken-pair excitations. This is because in the cluster operator defined in Equation (2), the α and β electron creation and annihilation operators are restricted to occupy the same orbitals. As a result, terms such as \({a}_{a\alpha }^{\dagger }{a}_{b\beta }^{\dagger }{a}_{j\beta }{a}_{i\alpha }(i\,\ne\,j,a\,\ne\,b)\) are not included. By neglecting these terms, we observe a quantitative deviation from the exact energy. To remedy this, we would like to recover the energetic corrections from these terms. To account for the broken-pair contributions we use perturbation theory as follows. We first partition the molecular Hamiltonian as$$H={H}_{0}+\lambda V$$
(5)
in which H0 is a simple non-interacting Hamiltonian and V accounts for the electron interaction. λ is the scale of the interaction.We then write the energy and wave function in terms of λ$$\begin{array}{l}\left\vert \Psi \right\rangle =\left\vert {\Psi }^{(0)}\right\rangle +\lambda \left\vert {\Psi }^{(1)}\right\rangle +{\lambda }^{2}\left\vert {\Psi }^{(2)}\right\rangle +\cdots \\ E={E}^{(0)}+\lambda {E}^{(1)}+{\lambda }^{2}{E}^{(2)}+\cdots \,.\\ \end{array}$$
(6)
Plugging the wave function and energy into the time-independent Schrödinger’s equation, and matching terms up to second order, one could derive the working equation for VQE-PT2. The second-order energy correction becomes$${E}^{(2)}=\mathop{\sum }\limits_{P}^{{\rm{non-pair}}}{t}_{P}\left\langle {\Psi }^{(0)}| V| {\Psi }_{P}\right\rangle$$
(7)
in which P is a compounded index for a pair-breaking single or double excited configuration, and tP is computed by solving,$$\sum _{P}{G}_{PQ}{t}_{P}=-{Y}_{Q}$$
(8)
where we define$$\begin{array}{l}{G}_{PQ}=\left\langle {\Psi }_{Q}| {H}_{0}| {\Psi }_{P}\right\rangle -{E}^{(0)}\left\langle {\Psi }_{Q}| {\Psi }_{P}\right\rangle \\ {Y}_{P}=\left\langle {\Psi }_{P}| V| {\Psi }^{(0)}\right\rangle \end{array}$$
(9)
in which the G matrix and the \(\overrightarrow{Y}\) vector are measured on quantum computers. The readers are strongly encouraged to read the “Method” section to find an in-depth explanation of the VQE-PT2 working equations derivations. The method for efficient construction of these quantities, as well as a regularization method to address the numerical instability of these equations, are provided in the Supplementary Information.Numerical results on molecular dissociationsWe begin our numerical results with the dissociation of the N2 triple bond. In Fig. 2, we compare the energy predicted by restricted Hartree-Fock (RHF), oo-upCCD, oo-upCCD-PT2, and full configuration interaction (FCI). We observe that the oo-upCCD energy provides a significant improvement over RHF, but it remains far from the FCI energy. This is due to the missing broken-pair excitation contributions in the oo-upCCD wave function. After applying the perturbation correction, we find that it is now much closer to FCI. The nonparallelity error (NPE) defined as the difference between the largest and smallest error with respect to FCI along the potential surface, is reduced from 52 mH to 14 mH, which demonstrates the effectiveness of the perturbation correction.Fig. 2: Potential energy surface for the dissociation of the N2 molecule in the STO-3G basis computed with RHF, oo-upCCD, oo-upCCD-PT2, and FCI.VQE results are obtained from a noise-free simulator.Our second example is the symmetric dissociation of the Li2O molecule. Li2O is one of the secondary reaction products in lithium-air batteries, which is a potential candidate for next-generation lithium battery due to its high energy density. The results on an ideal simulator are shown in Fig. 3, comparing RHF, FCI, oo-upCCD, and oo-upCCD-PT2. As one could see, similar to the N2 dissociation, the oo-upCCD-PT2 significantly improves the accuracy of oo-upCCD, and reduces the NPE vs FCI from 69 mH to 24 mH. This is particularly true in the equilibrium geometry. When bonds are stretched, oo-upCCD-PT2 still exhibits some noticeable amount of errors compared to FCI. These remaining errors are due to the limitation of the chosen form of the second order perturbation theory, such as 1) higher order terms are needed 2) one needs to use a different zeroth-order Hamiltonian than the simple one-body one we use. 3) the diagonal approximation of the G matrix.Fig. 3: Potential energy surface for the symmetric dissociation of the Li2O molecule in the STO-3G basis computed with RHF, oo-upCCD, oo-upCCD-PT2, and FCI.VQE results are obtained from a noise-free simulator.Numerical results on chemical reactionsWe now move our attention to the chemical decomposition process of the CH2OH+ → HCO++H2 and the SN2 reaction CH3I + Br− → CH 3Br + I−. The simulated energy profile along the reaction path is shown in Fig. 4 for the CH2OH+ decomposition. After freezing the core orbitals, the remaining eleven spatial molecular orbitals are mapped to eleven qubits. Unsurprisingly, oo-upCCD without perturbative corrections is insufficient to produce quantitative accuracy for absolute energies, capturing roughly only 50% of correlation energy. However, due to error cancellation, the predicted reaction energy barrier (209 mH) is close to FCI (202 mH). The energy difference between reactants and products (ΔE) is predicted to be 50 mH by oo-upCCD, which overestimates the FCI prediction (38 mH) by 24%. Applying the PT2 correction significantly improves the results, and we are able to capture 88% of the correlation. The predicted reaction energy barrier and ΔE is 205 mH and 39 mH respectively, which reduces the original oo-upCCD error by 57% and 92%.Fig. 4: The simulated CH2OH+ → HCO++H2 reaction energy pathway using RHF, oo-upCCD, oo-upCCD-PT2, and FCI.VQE results are obtained from a noise-free simulator. The zero of the energy is set to the FCI reactant state.The simulated results of the SN2 reaction of CH3I + Br− are shown in Fig. 5. After freezing the core orbitals, we are left with 15 spatial orbitals—therefore, 15 qubits—to simulate. Working in this active space, we compare with classical complete active space configuration interaction (CASCI) calculations. Similar to the CH2OH+ decomposition, the perturbation correction brings the oo-upCCD energy much closer to the CASCI predictions. Besides absolute energies, relative energies are also improved. The predicted barrier height is 14 mH by oo-upCCD-PT2 vs 10 mH by CASCI, which reduces the error of oo-upCCD (19 mH) by 50%.Fig. 5: The simulated CH3I + Br− → CH 3Br + I− SN2 reaction reaction energy pathway using RHF, oo-upCCD, oo-upCCD-PT2, and FCI.VQE results are obtained from a noise-free simulator. The zero of the energy is set to the FCI product state.Motivated by the success of perturbation theory for oo-upCCD in quantum simulators, we performed these chemical reaction simulations on two generations of the IonQ’s quantum computers. We simulate reaction steps 0, 50, and 140 of the CH2OH+ decomposition process on the Aria QPU. These three points correspond to the reactants, transition state, and products, respectively. We first perform a circuit pruning process, in which gates whose parameters are below a chosen threshold are removed from the circuit. In our experience, the minor energy benefits derived from these small-parameter quantum operations are overshadowed by the introduction of system noise. In our study, we choose the threshold to be 0.04 radians. The pruned circuits have 10 CX gates and 46 single qubit gates. The experimental results are shown in Fig. 6. As expected, depolarizing hardware noise yields a systematic, positive bias to the total energy. This can be alleviated by applying a constant shift to the energy. In this case, we shift by a constant 364 mH, which is obtained by computing the average energy error between the simulator and the quantum computer for these three structures. As we have demonstrated in a previous study27, in the small error regime, the errors lead to a positive bias in the measured VQE energy, and the bias is consistent across different geometries. Due to this consistency, the predicted relative energies are indeed very accurate. After applying this shift, the experimental results roughly match the prediction of the noiseless simulator. This suggests that the hardware errors are consistent across reaction steps and the relative energy is not affected.Fig. 6: The simulated CH2OH+ → HCO++H2 reaction energy pathway on the IonQ Aria quantum computer.The zero of the energy is set to the CASCI reactant state. Data markers are staggered for legibility, and error bars correspond to a statistical error of ± 1 standard error.Although it is tempting to consider the shifted energies as a form of error mitigation, we note that the exact shift requires the availability of noiseless simulation results, which is generally not the case for circuits that are not classically simulatable. Therefore, one should think of the shift as merely a way to demonstrate that the relative energies—which are of more utility in practice—are still accurate despite the presence of noise.The total energy of oo-upCCD-PT2 contains two parts: the oo-upCCD energy (EVQE) and the energy correction (E(2)). Table 1 shows the energy contributions from these two parts comparing the statevector simulator and the Aria QPU, and we find that almost all the errors in energy are from the the oo-upCCD energy term. This could be understood in two ways. First, the magnitude of the oo-upCCD energy is much larger than the energy correction, and so with the same error rate, the former would result in larger absolute errors than the latter. Second, the perturbative energy correction also benefits from the error cancellation, as the errors in the numerator (\({Y}_{p}^{2}\)) and the denominator (Gp) may cancel with each other. Therefore, the energy correction appears significantly more resilient to errors than the oo-upCCD energy. In Fig. 7 we compare the energy contributions to E(2) from unpaired excitations with \({({Y}_{p}/{G}_{p})}^{2}\) above 1 mH, using both the statevector simulator and the Aria QPU. The results from Aria closely align with those of the simulator.Table 1 Energy contributions of oo-upCCD-PT2 for CH2OH+ between the ideal simulator and the Aria QPUFig. 7: Comparison of top energy contributions to \(E^{(2)}\) between simulator and quantum hardware.The plot shows the energy contribution to \(E^{(2)}\) from the top 12 unpaired excitations for reaction step 140 of the CH2OH+ decomposition, comparing results from the statevector simulator and the Aria QPU.In order to improve the absolute energy measurements of oo-upCCD, we take advantage of a simple error mitigation approach based on post-selection. As we have mentioned in section II A, three measurements are needed to measure the energy of the oo-upCCD ansatz: measuring all the qubits in Pauli X, Pauli Y, and Pauli Z basis. The measurements in the Pauli Z basis should only yield binary strings that contains the correct number of electron pairs. For example, consider a 4-qubit, 2-electron pair system: binary strings such as 0011, 1010 are valid, while binary strings like 0111 or 0001 are invalid. On noisy quantum systems, these invalid binary strings could occur due to bit flip errors.Therefore, in the case of Pauli Z-basis (computational basis) measurements, we retain only those that maintain particle number symmetry and discard the rest. The results are shown in Fig. 6. We find that doing so improves the energy measurements by about 200 mH. We then perform the uniform shift for the post-selected energies by 151 mH. As in the case of CH2OH+ decomposition, the shift is computed as the averaged energy error between the measured results on hardware and the noiseless simulator. After shifting, the results better match the results from the ideal simulator. This demonstrates that the post-selection not only improves absolute energies, but also relative energies.The experimental results of the CH3I + Br− SN2 reaction are shown in Fig. 8 and presented in Table 2. These results are obtained using the IonQ Forte QPU. The SN2 reaction of CH3I + Br− is particularly challenging for NISQ quantum systems as the energy differences across the reaction path are small. The energy difference between the transition state and the reactant is about 50 mH, which is a value on the order of the experimental uncertainty. As seen in Fig. 8, the raw energy demonstrates nonphysical behavior. That is, the product energy is even higher than the transition state and the reactant energies. Similar to the CH2OH+ decomposition, the error in the total energy is predominantly due to errors from the unperturbed VQE energy. Therefore, we apply the same Z-basis post-selection, and find that the error mitigated energies now yield the correct behavior, matching well with the predictions of the simulator and CASCI.Fig. 8: The simulated CH3I + Br− → CH 3Br + I− SN2 reaction energy pathway on the IonQ Forte quantum computer.The zero of the energy is set to the simulated product state. Data markers are staggered for legibility, and error bars correspond to a statistical error of ± 1 standard error.Table 2 Energy contributions of oo-upCCD-PT2 for the SN2 reaction between the ideal simulator and the Forte QPU