Variational quantum eigensolver for approximate diagonalization of downfolded Hamiltonians using generalized unitary coupled cluster ansatz

As originally envisaged by Feynman [1], complex quantum system simulations should be performed using computational devices (quantum computers) operating based on principles of quantum mechanics and its stochastic character. The last few decades have been marked by the intensive development of quantum algorithms designed to tackle many-body problems characterized by complexities which cannot be effectively addressed even by employing computational power of near-exa-scale classical computer architectures. These difficulties originate in the necessity of systematic inclusion of collective many-body effects that quickly leads to unsurmountable computational barriers that preclude in many cases predictive simulations of strongly correlated systems in chemistry and material sciences. Quantum computing brings the hope that these problems can be overcome and reliable simulations of quantum systems can be performed in finite times. Another advantage of certain quantum computing algorithms (for example, quantum phase estimation (QPE) [2–9] algorithm) is the fact that they significantly reduce (although it does not entirely eliminate) the bias of the conventional approximate approaches, typically associated with the choice of the broadly understood zeroth-order approximations and levels of accuracy, and at least theoretically are capable of providing exact results at least in the stochastic sense.

The development of quantum algorithms for quantum chemistry have permeated several simulations including ground-state calculations [10, 11], excited-state modeling for low-lying valence excited states [12–17] and core-level transitions [18], phase estimation for relativistic and non-Born–Oppenheimer Hamiltonians [19, 20], imaginary time evolution and quantum Lanczos algorithms for thermal states [21–23], quantum filter diagonalization [24], quantum inverse iterations, [25], quantum power methods [26], and efficient preparation of correlated fermionic states [27, 28]. Nowadays, two major classes of quantum algorithms are being intensively developed with the aim of applying them to realistic systems on current noisy intermediate-scale quantum (NISQ) devices: (1) mentioned earlier family of QPE algorithms, and (2) broad class of VQE methods [13, 29–36]. These two algorithms are driven by different computational strategies and characterized by different resource requirements. While VQE draws heavily on the explicit form of the wave function ansatz $\vert {\Psi}\left( \overrightarrow {\theta }\right)\rangle$ with the polynomial number of variational parameters θ_i , the QPE utilizes unitary evolution operator e^−itH |Φ⟩ for judiciously chosen initial guess |Φ⟩ for the exact wave function.

The applicability of VQE and QPE algorithms is either limited by the number of variational variables and the numbers of necessary measurements or by the circuit depths, respectively, which was extensively discussed in the literature (see for example reference [37]). For the VQE formalism, the unitary coupled cluster (UCC) [38–46] expansion (as well as it generalized UCC variant (GUCC) [47]) is most frequently used form for the wave function ansatz, where the impact of higher-rank excitations and the role of the orbital choice [48, 49] have been explored. A non-orthogonal multi-reference variant of VQE logical ansatz (NOVQE) [35] has also been recently tested and verified for molecules characterized by the presence of strong correlation effects. On the other hand, the progress achieved in the development of QPE formalism should be mostly attributed to the development of novel Trotterization and Taylor expansion techniques [50–54] (or approximations of the exponential form of the unitary operator in general) and new representations of the second-quantized form of the many-body Hamiltonians stemming from the double decomposition of tensors describing pairwise interactions [37, 55, 56]. Using these techniques it was demonstrated that the overall number of terms entering the many-body Hamiltonian can be significantly reduced, which leads to a significant reduction in the gate depth of the corresponding quantum circuit.

In the last few years, several methodologies have been explored to reduced the dimensionality of many-body Hamiltonians used in quantum computing. This includes theory of transcorrelated Hamiltonians, canonical transformation theory, and subspace expansions [57–59]. In the case of the downfolding methods [60, 61] one can construct active-space representation of effective/downfolded Hamiltonians which integrate out high-energy Fermionic degrees of freedom while being capable of reproducing exact energy of quantum systems. This is accomplished by introducing an active space, based on energy, spatial, and/or time criterium, to define subsets of excitations either entirely within the active-space or involving some external orbital. The external excitations are used to define a similarity transformed effective Hamiltonian in which the corresponding correlation is captured and can be diagonalized to determine the remaining correlation. These formulations based on the double unitary CC (DUCC) ansatz have recently been used in the context of quantum computing for ground- and excited-state problems showing a promise in the recovering the exact energy of the full problem using limited quantum resources.

It is worth discussing the differences between the DUCC formalism and transcorrelated Hamiltonians approaches and quantum subspace algorithms. While the quantum algorithms for transcorrelated approaches [57] were introduced to counteract the deficiencies of the Gaussian basis sets in restoring cusp-condition (and effectively use smaller basis sets), the quantum sub-space approximations [14, 59, 62] are based on the idea of extending a simple form of the VQE wave function represented in active space by probing configurations inside and outside of active space, which results in the form of selected configuration interaction formalism. The sub-space algorithms can be considered as 'bottom-up' approaches. In contrast to the sub-space methods, the DUCC formalism provides a rigorous many-body prescription on compressing the full Hamiltonian into an effective Hamiltonian in a way that ground-state energy can be recovered by diagonalizing effective Hamiltonian in active space. In practical applications, the DUCC formalism provides a basis for well-justified approximations. The DUCC formalism is a 'top-down' approach or a hybrid of Löwdin partitioning technique [63, 64] that takes full advantage of properties of exponential parametrization of the wave function. For example, if the active space is properly dissociating, the DUCC formalism leads to size-consistent energies (which is not the case for selected configuration interaction methods). Despite the seemingly similar structure, the DUCC formalism is also different from the canonical transformation theory (CTT) [65]. In the former type of approach, the appearance of the effective/downfolded Hamiltonians is a natural consequence of the double exponential ansatz for the exact wavefunction. Moreover, DUCC is a single reference formalism in its nature, which results in a different set of parameters used to describe the wave function compared to CTT formalism.

For further development of various downfolding methods a crucial role is played by the possibility of integrating them with various ansatzen that can provide optimal utilization of quantum VQE algorithm. For this reason, in this paper we propose a variant of the VQE algorithm that utilize generalized UCC ansatz to diagonalize downfolded Hamiltonian in active spaces. We will refer the resulting formalism to as the downfolding VQE algorithm with GUCC form of trial wave function or DVQE-GUCC for short. We also evaluate the performance of various classical VQE optimizers for systems characterized by various complexity levels of ground-state electronic wave functions. In particular, we compare DVQE-GUCC results for three benchmark systems (H₂O, C₂H₄, and N₂) with the exact results obtained by the full configuration interaction (FCI) method.

In reference [60], we introduced the unitary extension of the sub-system embedding sub-algebra CC approach (SES-CC) [66] which utilizes the double unitary CC expansion. The Hermitian form of CC downfolding is predicated on a special form of wave function expansion given by double unitary coupled cluster ansatz (DUCC) given by the expressions (for more details see reference [60])

$$\begin{equation}\vert {\Psi}\rangle ={\mathrm{e}}^{{\sigma }_{\text{ext}}}\enspace {\mathrm{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle ;,\end{equation} \tag{ 1 }$$

where σ_int and σ_ext are the general type anti-Hermitian operators

$$\begin{equation}{\sigma }_{\text{int}}^{{\dagger}}=-{\sigma }_{\text{int}},\end{equation} \tag{ 2 }$$

$$\begin{equation}{\sigma }_{\text{ext}}^{{\dagger}}=-{\sigma }_{\text{ext}}\end{equation} \tag{ 3 }$$

defined by amplitudes defining action within and outside of the pre-defined active space, respectively, i.e., the amplitudes defining the σ_ext operator must carry at least one inactive spin-orbital index whereas all amplitudes defining the σ_int operator carry active spinorbital indices only. In equation (1), |Φ⟩ designates properly chosen reference function (usually chosen as a Hartree–Fock (HF) Slater determinant). The exactness of the expansion (1) has been recently discussed in reference [61] where it was also shown that the standard UCC expansions can provide a basic approximation of the exact σ_int and σ_ext operators, i.e.,

$$\begin{equation}{\sigma }_{\text{int}}\simeq {T}_{\text{int}}-{T}_{\text{int}}^{{\dagger}},\end{equation} \tag{ 4 }$$

$$\begin{equation}{\sigma }_{\text{ext}}\simeq {T}_{\text{ext}}-{T}_{\text{ext}}^{{\dagger}},\end{equation} \tag{ 5 }$$

where T_int and T_ext are single-reference-type internal and external cluster amplitudes (in the sense defined above).

Introducing (1) into Schrödinger's equation, premultiplying both sides by ${\text{e}}^{-{\sigma }_{\text{int}}}\enspace {\text{e}}^{-{\sigma }_{\text{ext}}}$, and projecting the resulting equations onto subspace of excited determinants, with respect to the reference function |Φ⟩, defined by a projection operator Q and onto reference function (with corresponding projection operator P = |Φ⟩⟨Φ|, one obtains equations for cluster amplitudes and the corresponding energy, i.e.,

$$\begin{equation}Q{\text{e}}^{-{\sigma }_{\text{int}}}\enspace {\text{e}}^{-{\sigma }_{\text{ext}}}H{\text{e}}^{{\sigma }_{\text{ext}}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle =0,\end{equation} \tag{ 6 }$$

$$\begin{equation}\langle {\Phi}\vert {\mathrm{e}}^{-{\sigma }_{\text{int}}}\enspace {\mathrm{e}}^{-{\sigma }_{\text{ext}}}H{\mathrm{e}}^{{\sigma }_{\text{ext}}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle =E,\end{equation} \tag{ 7 }$$

Using the Campbell–Baker–Hausdorff formula, the above formulas can be put in explicitly connected form

$$\begin{equation}Q{\left({\mathrm{e}}^{-{\sigma }_{\text{int}}}\enspace {\mathrm{e}}^{-{\sigma }_{\text{ext}}}H\enspace {\text{e}}^{{\sigma }_{\text{ext}}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}\right)}_{\mathrm{C}}\vert {\Phi}\rangle =0,\end{equation} \tag{ 8 }$$

$$\begin{equation}\langle {\Phi}\vert {\left({\mathrm{e}}^{-{\sigma }_{\text{int}}}\enspace {\mathrm{e}}^{-{\sigma }_{\text{ext}}}H\enspace {\text{e}}^{{\sigma }_{\text{ext}}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}\right)}_{\mathrm{C}}\vert {\Phi}\rangle =E,\end{equation} \tag{ 9 }$$

where subscript 'C' stands for a connected part of a given operator expression. In contrast to standard CC formulations, expansions in above equations are non-terminating.

For the sake of simplicity, in the following analysis we assume the case of the exact limit (σ_int and σ_ext include all possible excitations). In reference [60], we showed that when σ_int contains all possible excitations/de-excitations within the active space, the energy of the system section 7 can be obtained by diagonalizing the DUCC effective Hamiltonian

$$\begin{equation}{\bar{H}}_{\text{eff}}^{\left(\text{DUCC}\right)}\enspace {\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle =E\enspace {\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle ,\end{equation} \tag{ 10 }$$

where ${\bar{H}}_{\text{eff}}^{\left(\text{DUCC}\right)}$ is defined as

$$\begin{equation}{\bar{H}}_{\text{eff}}^{\left(\text{DUCC}\right)}=\left(P+{Q}_{\text{int}}\right){\bar{H}}_{\text{ext}}^{\text{DUCC}}\left(P+{Q}_{\text{int}}\right)\end{equation} \tag{ 11 }$$

and

$$\begin{equation}{\bar{H}}_{\text{ext}}^{\text{DUCC}}={\text{e}}^{-{\sigma }_{\text{ext}}}H\enspace {\text{e}}^{{\sigma }_{\text{ext}}}.\end{equation} \tag{ 12 }$$

In the above equations, the ${\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle$ vector defines the corresponding eigenvector. To prove this property, it is sufficient to introduce the resolution of identity ${\text{e}}^{{\sigma }_{\text{int}}}\enspace {\text{e}}^{-{\sigma }_{\text{int}}}$ to the left of the ${\bar{H}}_{\text{ext}}^{\text{DUCC}}$ operator in

$$\begin{equation}\left(P+{Q}_{\text{int}}\right){\bar{H}}_{\text{ext}}^{\text{DUCC}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle =E\left(P+{Q}_{\text{int}}\right){\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle ,\end{equation} \tag{ 13 }$$

where we employed the fact that

$$\begin{equation}\left(P+{Q}_{\text{int}}\right){\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle ={\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle ,\end{equation} \tag{ 14 }$$

and to notice that ${\text{e}}^{-{\sigma }_{\text{int}}}{\bar{H}}_{\text{ext}}^{\text{DUCC}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}={\text{e}}^{-{\sigma }_{\text{int}}}\enspace {\text{e}}^{-{\sigma }_{\text{ext}}}H\enspace {\text{e}}^{{\sigma }_{\text{ext}}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}$. Next, using matrix representation of the σ_int operator in the CAS space, denoted as σ _int, this equation can be re-written as

$$\begin{equation}\left[{\mathrm{e}}^{{\boldsymbol{\sigma }}_{\text{int}}}\right]\left[\boldsymbol{y}\right]=0,\end{equation} \tag{ 15 }$$

where the first component of the [ y ] vector is equivalent to $\langle {\Phi}\vert {\text{e}}^{-{\sigma }_{\text{int}}}\enspace {\text{e}}^{-{\sigma }_{\text{ext}}}H\enspace {\text{e}}^{{\sigma }_{\text{ext}}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle -E$ while the remaining components correspond to projections of ${\text{e}}^{-{\sigma }_{\text{int}}}\enspace {\text{e}}^{-{\sigma }_{\text{ext}}}H\enspace {\text{e}}^{{\sigma }_{\text{ext}}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle$ onto excited configurations belonging to Q_int. The $\left[{\mathrm{e}}^{{\boldsymbol{\sigma }}_{\text{int}}}\right]$ matrix is also non-singular, which is a consequence of the formula

$$\begin{equation}\enspace \mathrm{det}\left({\mathrm{e}}^{{\boldsymbol{\sigma }}_{\text{int}}}\right)={\mathrm{e}}^{\mathrm{Tr}\left({\boldsymbol{\sigma }}_{\text{int}}\right)}=1\end{equation} \tag{ 16 }$$

and the anti-Hermitian character of the σ _int matrix, i.e., Tr( σ _int) = 0 (where real character of σ_int cluster amplitudes is assumed). Given the non-singular character of the $\left[{\mathrm{e}}^{{\boldsymbol{\sigma }}_{\text{int}}}\right]$ matrix (see also reference [60]), this proves the equivalence of these two representations.

The discussed variant of double unitary CC expansion extends properties of single reference CC sub-system embedding sub-algebra (SES-CC) [66] to the unitary CC case. In analogy to SES-CC case, the DUCC expansion based on ansatz (1) leads to the rigorous decoupling of fermionic degrees of freedom in the effective Hamiltonian. As discussed in earlier papers these degrees of freedom can be associated with energy [60], time [67], and spatial scales [60]. While the focus of this and previous work has been on the Hamiltonian and associated ground-state energy, the ground-state expectation value of any operator A (⟨A⟩) can be calculated if the ${\text{e}}^{{\sigma }_{\text{int}}}$ operator is known or can be effectively approximated. For this purpose one can use the formula

$$\begin{equation}\langle A\rangle =\langle {\Phi}\vert {\text{e}}^{-{\sigma }_{\text{int}}}{\bar{A}}_{\text{ext}}^{\text{DUCC}}\enspace {\text{e}}^{{\sigma }_{\text{int}}}\vert {\Phi}\rangle ,\end{equation} \tag{ 17 }$$

where ${\bar{A}}_{\text{ext}}^{\text{DUCC}}={\text{e}}^{-{\sigma }_{\text{ext}}}A\enspace {\text{e}}^{{\sigma }_{\text{ext}}}$.

In the following part of the paper we will study the approximate form of the ${\bar{H}}_{\text{eff}}^{\text{DUCC}}$ operator. In analogy to reference [60] we will use the approximate form of the ${\bar{H}}_{\text{eff}}^{\text{DUCC}}$, which is consistent with the second order energy expansion, i.e.,

$$\begin{equation}{\bar{H}}_{\text{eff}}^{\text{DUCC}}\simeq \left(P+{Q}_{\text{int}}\right)\left(H+\left[{H}_{\mathrm{N}},{\sigma }_{\text{ext}}\right]+\frac{1}{2!}\left[\left[{F}_{\mathrm{N}},{\sigma }_{\text{ext}}\right],{\sigma }_{\text{ext}}\right]\right)\left(P+{Q}_{\text{int}}\right),\end{equation} \tag{ 18 }$$

where H_N and F_N are the Hamiltonian and Fock operators in the normal product form and the cluster operator σ_ext is approximated as

$$\begin{equation}{\sigma }_{\text{ext}}\simeq {T}_{\text{ext,}\;1}+{T}_{\text{ext,}\;2}-{T}_{\text{ext,}\;1}^{{\dagger}}-{T}_{\text{ext,}\;2}^{{\dagger}},\end{equation} \tag{ 19 }$$

where T_ext,1 and T_ext,2 are external pars of the CCSD T₁ and T₂ operators. Moreover, ${\bar{H}}_{\text{eff}}^{\text{DUCC}}$ is approximated by one- and two-body effects as discussed in reference [60], since several algorithms for diagonalizing ${\bar{H}}_{\text{eff}}^{\text{DUCC}}$ are implemented with only one- and two-body elements of the Hamiltonian. In this work, we employ these approximations with an active space consisting of the occupied orbitals and lowest-lying virtual orbitals, where a large portion of the static correlation is contained. That means that correlation captured by the DUCC Hamiltonian is mainly dynamical, which is characterized by relatively small cluster amplitudes and is suitably captured by the underlying CCSD. As we have demonstrated in references [60, 68], and will once again show in this work, this series of approximations often leads to a recovery of more than 85% of the correlation energy. Errors can be improved upon by extending these approximations since the DUCC formalism is exact in the limit, or reshaping the correlation landscape as in the case when natural orbitals are used [68].

Although the quantum phase estimation algorithm [2–9] provide efficient means to obtain the exact full CI solution on a quantum computer, as far as presently understood, it requires large number of qubits, with quantum error correction, and long coherence times, which are (except for tiny model systems) beyond the reach of present quantum hardware. A promising alternative for near-term quantum devices is the variational quantum eigensolver (VQE) approach [13, 29–36], which combines classical variational energy minimization over normalized trial wave functions Ψ(ϑ) parameterized by a set of variables ϑ

$$\begin{equation}E={\mathrm{min}}_{\vartheta }\enspace \langle {\Psi}\left(\vartheta \right)\vert H\vert {\Psi}\left(\vartheta \right)\rangle \end{equation} \tag{ 20 }$$

with computation of the Hamiltonian matrix elements on a quantum device. The 'bare' second-quantized molecular Hamiltonian can be expressed as a sum of one- and two-body terms, formally

$$\begin{equation}H=\sum\limits _{X}{h}_{X}\hat{X}\end{equation} \tag{ 21 }$$

where $\hat{X}$ is a string of up to two creation and two annihilation operators. The DUCC effective Hamiltonian (18) can be approximated in the same form, when higher-body terms are neglected. In order to minimize the complexity of the quantum circuit, expectation values of each term in the summation (21) are evaluated separately and results added classically. The advantage of this approach consists in the fact that evaluation of these matrix elements for many trial wave functions, which would lead to exponential complexity on a classical computer, can be efficiently performed on a quantum computer with a moderate circuit complexity and coherence time requirements [69].

A popular wave function ansatz is the exponential CC form

$$\begin{equation}{\Psi}\left(t\right)={\mathrm{e}}^{T}\vert {\Phi}\rangle \end{equation} \tag{ 22 }$$

where at the singles and doubles truncation level the cluster operator has the form

$$\begin{equation}T={T}_{1}+{T}_{2}\end{equation} \tag{ 23 }$$

with

$$\begin{equation}{T}_{1}=\sum\limits _{ia}{t}_{i}^{a}{a}_{a}^{{\dagger}}{a}_{i}\end{equation} \tag{ 24 }$$

$$\begin{equation}{T}_{2}=\sum\limits _{i{< }j,a{< }b}{t}_{ij}^{ab}{a}_{a}^{{\dagger}}{a}_{b}^{{\dagger}}{a}_{j}{a}_{i}\end{equation} \tag{ 25 }$$

where we adopt the convention of i, j, ... running over occupied spinorbitals and a, b, ... running over the virtual ones. In the traditional CC method [70] this ansatz is not employed in the variational minimization (20), since that yields infinite diagrammatic expansions, which are inconvenient to evaluate and efficiently approximate on a classical computer. A projection technique is employed instead, which leads to closed explicit algebraic equations for the parameters ${t}_{i}^{a}$, ${t}_{ij}^{ab}$. For the same reason, the unitary coupled cluster ansatz (UCCSD)

$$\begin{equation}{\Psi}\left(t\right)={\mathrm{e}}^{T-{T}^{{\dagger}}}\vert {\Phi}\rangle \end{equation} \tag{ 26 }$$

was never widely employed in mainstream quantum chemistry, although considerable effort has been dedicated to development of computationally efficient approximations [71]. On the other hand, on a quantum computer the unitary ansatz is clearly advantageous as the wave function parametrization then straightforwardly and efficiently translates to a (unitary) quantum circuit. However, it still bears the main disadvantage of the single-reference CC method, that it is based on truncated excitations from a single reference determinant Φ, which hampers it application to strongly correlated (multiconfigurational) systems.

A possible generalization of the CCSD ansatz (22) consists in including both excitations and deexcitations into the cluster operator, which thus attains the GCCSD form

$$\begin{equation}{T}_{1}=\sum\limits _{pq}{t}_{q}^{p}{a}_{p}^{{\dagger}}{a}_{q}\end{equation} \tag{ 27 }$$

$$\begin{equation}{T}_{2}=\sum\limits _{p{< }q,r{< }s}{t}_{rs}^{pq}{a}_{p}^{{\dagger}}{a}_{q}^{{\dagger}}{a}_{s}{a}_{r},\end{equation} \tag{ 28 }$$

where the indices p, q, r, s run over all spinorbitals and the reference to a Fermi vaccum is thus eliminated from the cluster operator. Moreover, the 'reference' wave function Φ in (22) need not be a single Slater determinant, but any state that can be efficiently prepared on a quantum computer. In the early 2000s this ansatz was subject of an active debate in the literature [72–84], since Nooijen conjectured [72] that the ground state solution of a two-body Hamiltonian can be exactly expresed in this form. This conjecture was later disproved by Kutzelnigg and Mukherjee based on Lie-algebraic arguments [83, 84] (after excluding solutions in the infinite t limit where e^T is a projection operator and precisely defining what is allowed to be the reference wave function Φ). Nevertheless, several studies in the meantime have shown, that even if not exact, it yielded excellent numerical agreement compared to FCI [73, 78]. This has recently motivated its use in the VQE context [69], in the unitary modification (GUCCSD), when the GCCSD cluster operator (27) is combined with the unitary ansatz (26). Notice that since in the GCCSD both excitations and deexcitations are included, the antihermitian linear combination T − T^† is also expressible in the GCCSD form (27) where the unitarity manifests itself only in restrictions on the values of the amplitudes of mutually conjugated excitation operators, so that GUCCSD trial wave function has the same form as GCCSD, just with less free parameters, while the unitarity represents an advantage for the implementation in a quantum circuit. The mutual non-commutativity of the terms in the GCCSD cluster operator (27) has two implications. Firstly, in the qunatum implementation it has to be approximated by a Trotter expansion. However, previous numerical study has achieved excellent accuracy even with very low number of Trotter steps [69]. Secondly, it complicates the evaluation of the analytic gradient of the energy expectation value (20), since one cannot straightforwardly apply the chain rule. This problem has already been discussed and solved by Van Voorhis and Head–Gordon [73] using the Wilcox identity [85].

In our numerical study, we employed the DUCC effective Hamiltonian truncated to one- and two-body terms (18) and performed the variational minimization (20) using the BFGS method with analytic gradient. The GCCSD ansatz (27) has been implemented using an in-house code based on full-CI expansion of the excitation operators and numerical evaluation of the exponential (26) in a 'direct-CI' style. We have employed several choices of the reference state, either a single Slater determinant or matrix produxt state (MPS) wave function from a density matrix renormalization group (DMRG) calculations with several small values of bond dimension. Also the effect of different initial guess of amplitudes (zero, MP2, MPS) on the convergence rate has been investigated.

The C_2v -symmetric double dissociation of H₂O is a popular benchmark system to test the accuracy methods under different regimes of static/dynamical correlation, since the degree of static correlation is systematically increased while stretching both bonds. At equilibrium, simple truncation of the bare orbital space down to 7 or 9 orbitals disposes of over 85% of the correlation energy and results in errors of 0.185–0.211 Hartree, relative to the full configuration interaction (FCI) values (see table 1 and figure 1). As the bonds are stretched, the static correlation grows and the truncated orbital space recovers a larger percentage of the correlation. However, even when both bonds are stretched to three times the equilibrium length, over 25% of the correlation is still missing, amounting to 0.149–0.165 Hartree errors. When correlation from outside of the active space is incorporated using the DUCC downfolding procedure, 87%–93% of the correlation is consistently captured along the double-dissociation pathway. Errors for the 7 orbital active space are reduced to 0.019–0.063 Hartree, and are further reduced to 0.015–0.055 Hartree for the 9 orbital active space.

Zoom In Zoom Out Reset image size

In table 2 we performed a comparison of VQE with different types of wave function ansatz. In particular, we employed the standard unitary CCSD (UCCSD) and the generalized unitary CCSD one (GUCCSD), which contains also deexcitations and has thus more variational parameters. Except for the smallest orbital spaces, the number of GUCCSD parameters is significantly smaller than the number of FCI coefficients, though. In particular, in the largest 12-orbital space with 1 orbital frozen, the number of UCCSD, GUCCSD, and FCI variational parameters is ratio of about 1:10:100. Besides the total energies, we computed also overlaps of the trial wave function with the FCI one and the norm of difference between them, to assess the accuracy of the particular ansatz. As can be seen, the generalized ansatz yield energies much closer to the FCI one (sub-millihartree accuracy) than the regular UCCSD ansatz, where differences in the order of millihartrees appear. The superiority of the GUCCSD Ansatz is more pronounced at the stretched geometries, where the system has stronger multiconfigurational character. The overlaps and difference norms represent even stricter criteria than the energies only. For GUCCSD, overlaps differ from one at the 6th decimal place, while for UCCSD already at the 3rd decimal. Similarly, the norm of wave function difference from FCI is at least one order of magnitude better for GUCCSD. We thus consider the GUCCSD Ansatz advantageous in the VQE context, particularly for molecules far from the equilibrium geometry or in other multireference situations where no single Slater determinant dominates the wave function. Although GUCCSD is not equivalent to FCI, it is not biased toward any single Slater determinant, and its parameter number matches that of a two-body Hamiltonian, making it at least a plausible option. Our numerical data suggest that GUCCSD is able to yield energies as well as wave functions close to the FCI solution, certainly superior to the UCCSD ansatz. The price one has to pay is a larger number of variational parameters than in UCCSD and a slower convergence.

We investigated the performance of the DUCC downfolding procedure for the twisting of ethylene using two active spaces, 9 and 13 orbitals. As seen in table 3 and figure 2, the bare Hamiltonian recovers 6%–18% of the correlation along the potential energy surface for the 9 orbital active space and 10%–23% for the 13 orbital active space, with errors of 0.28–0.30 and 0.29–0.31 Hartree, respectively, when compared with the benchmark CCSDTQ calculations. When the correlations from outside of the active space are included through the DUCC downfolding procedure, 76%–86% of the correlation is recovered with either active space. In general, the DUCC Hamiltonians recover a more consistent portion of the correlation energy, with the largest errors occurring at the 90 degree twisted confirmation. The CCSDTQ curve shows the flattening of the potential energy surface at this point. The CCSD calculations that underlies the downfolding procedure do not produce the correct curvature, which propagates into the corresponding DUCC Hamiltonians. This demonstrates the importance of the source and quality of the amplitudes defining the downfolding procedure and the occasional need for amplitudes from higher-level methods. If the correlation effects inside and outside of the active space are well decoupled, the diagonalization of the DUCC Hamiltonian can overcome deficiencies in the underlying CCSD calculation, as we will see with N₂. It is important to mention that the active-space FCI calculations appear to have the correct curvature, but this is fortuitous as even CCSD will have similar curvature in the same active space.

Zoom In Zoom Out Reset image size

The bond breaking of N₂ is particularly difficult for conventional CC methods. Lower-order methods such as CCSD and even CCSDT have breakdowns that lead to nonphysical barriers (see figure 3), and methods such as CCSDTQ will run into difficulty converging at stretched bond lengths. The static and dynamical correlations are well decoupled along the potential energy surface when defined by an active space of just 10 orbitals. However, the dynamical correlation missing from FCI calculations in the active space is significant, especially near equilibrium. When those dynamical correlations are incorporated into the active space via the DUCC Hamiltonian, we recover energies comparable to the full conventional CC methods, especially at equilibrium. Errors at stretched bond lengths will be higher, but considering that the underlying approach to the downfolding technique is CCSD, this is rationalized. Once again, since the correlation effect in and out of the active space are well decoupled, the DUCC Hamiltonian results do not exhibit the breakdown that lower-order canonical CC methods, such as CCSD and CCSDT, have. With amplitudes from a more suited method and longer expansions defining the DUCC Hamiltonian will further improve upon the energies at the stretched bond lengths. Similarly to the results for water molecule, VQE yields more accurate results using the GUCCSD ansatz than with regular UCCSD. The former agree with FCI to sub-millihartree accuracy, while regular UCCSD exhibits errors in the millihartree range, particularly for the stretched bond lengths.

Zoom In Zoom Out Reset image size

In tables 4 and 5 we present VQE results for N₂ molecule with bare Hamiltonian in STO-3G basis to investigate its convergence for different reference kets and inital guess of GUCCSD amplitudes. These results correspond to computations with the frozen core approximation, i.e. correlate 10 electrons in 8 orbitals [complete active space, CAS(10,8)]. We varied the reference ket state between HF and three different MPS states and used initial amplitudes either zero, or from MP2, or from a CC analysis of the corresponding MPS wave function [86]. We would like to note that MPS wave functions with bounded bond dimensions can be efficiently prepared on a quantum register [87]. We employed MPSs with small bond dimensions equal to 8, 16, and 32. These wave functions were obtained by DMRG calculations [88] with a fixed number of sweeps equal to 6 and U(1) symmetries enforcing M_S = 0 spin projection. The initial ordering of orbitals was obtained by the Fiedler method [89] and all calculations were initialized with the CI-DEAS procedure [90, 91]. As the convergence criterium for the BFGS quasi-Newton energy minimization we employed maximum amplitude change 1 × 10⁻⁸ a.u. or maximum gradient component 1 × 10⁻⁵ a.u.

The first thing to note is that for both geometries, when HF is employed as a reference, the number of iterations before reaching convergence is relatively low compared to when the MPS8 or MPS16 are used. However, the final converged values have noticeable errors in energy compared to the FCI reference and relatively large errors in the norm of difference from FCI and/or overlap with FCI. On the other hand, for MPS32 the number of iterations till convergence is lower and accuracy higher, which can be explained by the fact that for such a small system the bond dimension 32 is already high enough for the reference being close to the exact solution. Nevertheless, the DMRG energies corresponding to the MPS32 guess states have error still higher than 1 mHartree.

Generally, MPS initial guess offers a better agreement with FCI, but this can comes at the cost of several more iterations before convergence in some cases. As the MPS gets larger, the accuracy generally increases toward FCI. However, we observed two 'pathological' cases where the method converged to a local minimum providing a solution with very high energy and very far from the ground state FCI wave function. This behavior seems sporadic, and it might disappear for other values of convergence thresholds, but we would like to understand this further in future studies. The choice of initial amplitudes played a little role in the convergence and achieved accuracy, except for the two aforementioned cases when the minimization did not converge to the correct ground state with one value of initial amplitudes, while it did for all the other ones employed.

In this paper, we have presented the quantum approach for approximate diagonalization of DUCC down-folded Hamiltonians, which combines the variational quantum eigensolver with the generalized unitary coupled cluster ansatz and is aimed at near-term noisy quantum devices. We have demonstrated its capability to capture a large part of dynamical electron correlation missing in the active space (about 80%–90% of total correlation) by numerical simulations of three benchmark systems of varying multireference character, namely N₂, H₂O, and C₂H₄. Last but not least, we have shown on the example of N₂, how the initial guesses in the form of matrix product states, which can be efficiently prepared on a quantum register, can improve the accuracy of the VQE-GUCC method for multireference systems.

The data that support the findings of this study are available upon reasonable request from the authors.

This work was supported by the 'Embedding Quantum Computing into Many-body Frameworks for Strongly Correlated Molecular and Materials Systems' project, which is funded by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, the Division of Chemical Sciences, Geosciences, and Biosciences. A portion of the calculations have been performed using the Molecular Science Computing Facility (MSCF) in the Environmental Molecular Sciences Laboratory (EMSL) at the Pacific Northwest National Laboratory (PNNL). PNNL is operated for the U.S. Department of Energy by the Battelle Memorial Institute under Contract DE-AC06-76RLO-1830. We also thank the Czech Ministry of Education for support of the Czech-US bilateral cooperation under the project no. LTAUSA17033 [88].