Abstract
In this paper, we discuss the utilization of variational quantum solver (VQE) and recently introduced generalized unitary coupled cluster (GUCC) formalism for the diagonalization of downfolded/effective Hamiltonians in active spaces. In addition to effective Hamiltonians defined by the downfolding of a subset of virtual orbitals we also consider their form defined by freezing core orbitals, which enables us to deal with larger systems. We also consider various solvers to identify solutions of the GUCC equations. We use N2, H2O, and C2H4, as benchmark systems to illustrate the performance of the combined framework.
Export citation and abstract BibTeX RIS
1. Introduction
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 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 (H2O, C2H4, and N2) with the exact results obtained by the full configuration interaction (FCI) method.
2. CC downfolding techniques
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])
where σint and σext are the general type anti-Hermitian operators
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.,
where Tint and Text 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 , 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.,
Using the Campbell–Baker–Hausdorff formula, the above formulas can be put in explicitly connected form
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
where is defined as
and
In the above equations, the vector defines the corresponding eigenvector. To prove this property, it is sufficient to introduce the resolution of identity to the left of the operator in
where we employed the fact that
and to notice that . Next, using matrix representation of the σint operator in the CAS space, denoted as σ int, this equation can be re-written as
where the first component of the [ y ] vector is equivalent to while the remaining components correspond to projections of onto excited configurations belonging to Qint. The matrix is also non-singular, which is a consequence of the formula
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 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 operator is known or can be effectively approximated. For this purpose one can use the formula
where .
In the following part of the paper we will study the approximate form of the operator. In analogy to reference [60] we will use the approximate form of the , which is consistent with the second order energy expansion, i.e.,
where HN and FN are the Hamiltonian and Fock operators in the normal product form and the cluster operator σext is approximated as
where Text,1 and Text,2 are external pars of the CCSD T1 and T2 operators. Moreover, is approximated by one- and two-body effects as discussed in reference [60], since several algorithms for diagonalizing 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].
3. VQE algorithm and generalized unitary CCSD ansatz
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 ϑ
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
where 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
where at the singles and doubles truncation level the cluster operator has the form
with
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 , . For the same reason, the unitary coupled cluster ansatz (UCCSD)
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
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 eT 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.
4. Results and discussion
The C2v -symmetric double dissociation of H2O 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.
Table 1. Energies obtained with various methods using both bare and DUCC Hamiltonians for the double dissociation of H2O
Hamiltonian: | HF | FCI | FCI | HF | HF | VQE | VQE | FCI | FCI | CCSD | FCI |
---|---|---|---|---|---|---|---|---|---|---|---|
Bare | Bare | Bare | DUCC | DUCC | DUCC | DUCC | DUCC | DUCC | Bare | Bare | |
Orbitals: | [7] | [9] | [7] | [9] | [7] | [9] | [7] | [9] | [24] | [24] | |
Eq | 0.2178 | 0.2113 (3) | 0.1846 (15) | 0.0253 (88) | 0.0429 (80) | 0.0189 (91) | 0.0148 (93) | 0.0189 (91) | 0.0148 (93) | 0.0037 (98) | −76.2419 |
1.5*Eq | 0.2700 | 0.2016 (25) | 0.1698 (37) | 0.0925 (66) | 0.1064 (61) | 0.0308 (89) | 0.0216 (92) | 0.0308 (89) | 0.0216 (92) | 0.0100 (96) | −76.0723 |
2*Eq | 0.3640 | 0.1778 (51) | 0.1603 (56) | 0.2268 (38) | 0.2265 (38) | 0.0468 (87) | 0.0372 (90) | 0.0468 (87) | 0.0372 (90) | 0.0220 (94) | −75.9517 |
2.5*Eq | 0.4767 | 0.1679 (65) | 0.1532 (68) | 0.3870 (19) | 0.3671 (23) | 0.0595 (88) | 0.0503 (89) | 0.0595 (88) | 0.0503 (89) | 0.0203 (96) | −75.9180 |
3*Eq | 0.5676 | 0.1654 (71) | 0.1491 (74) | 0.5106 (10) | 0.4781 (16) | 0.0633 (89) | 0.0545 (90) | 0.0633 (89) | 0.0545 (90) | 0.0108 (98) | −75.9119 |
aFCI energy values in the last column are in Hartree. All other quantities are errors relative to these energies. Values in parentheses are the percentages of correlation recovered.
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.
Table 2. Comparison of the GUCCSD and UCCSD VQE Ansatz for the double dissociation of H2O using DUCC Hamiltonian
Orbitals | E(VQE/UCCSD) | FCI overlap | FCI difference norm | E(VQE/GUCCSD) | FCI overlap | FCI difference norm | E(FCI) | |
---|---|---|---|---|---|---|---|---|
Eq | 7 | −76.222 874 | 0.999 984 | 0.005 53 | −76.222 952 | 0.999 999 | 0.000 53 | −76.222 954 |
1.5*Eq | 7 | −76.040 794 | 0.999 695 | 0.024 67 | −76.041 552 | 0.999 999 | 0.000 55 | −76.041 552 |
2*Eq | 7 | −75.904 402 | 0.999 704 | 0.024 30 | −75.904 863 | 0.999 998 | 0.001 44 | −75.904 866 |
2.5*Eq | 7 | −75.857 429 | 0.998 119 | 0.061 32 | −75.858 491 | 0.999 999 | 0.000 52 | −75.858 492 |
3*Eq | 7 | −75.606 818 | 0.999 939 | 0.010 98 | −75.606 867 | 0.999 999 | 0.000 34 | −75.606 867 |
Eq | 9 | −76.226 856 | 0.999 968 | 0.007 88 | −76.227 026 | 0.999 999 | 0.000 79 | −76.227 028 |
1.5*Eq | 9 | −76.050 040 | 0.999 716 | 0.023 83 | −76.050 791 | 0.999 998 | 0.001 41 | −76.050 796 |
2*Eq | 9 | −75.913 856 | 0.999 619 | 0.027 58 | −75.914 443 | 0.999 998 | 0.001 71 | −75.914 448 |
2.5*Eq | 9 | −75.866 377 | 0.997 809 | 0.066 19 | −75.867 658 | 0.999 999 | 0.001 05 | −75.867 660 |
3*Eq | 9 | −75.611 251 | 0.999 676 | 0.025 44 | −75.612 971 | 0.999 999 | 0.000 59 | −75.612 972 |
Eq | 12–1 | −76.226 562 | 0.999 882 | 0.015 37 | −76.227 370 | 0.999 997 | 0.002 28 | −76.227 394 |
1.5*Eq | 12–1 | −76.052 027 | 0.999 473 | 0.032 44 | −76.054 390 | 0.999 997 | 0.002 14 | −76.054 404 |
2*Eq | 12–1 | −75.922 196 | 0.998 821 | 0.048 55 | −75.926 965 | 0.999 998 | 0.001 98 | −75.926 976 |
2.5*Eq | 12–1 | −75.877 720 | 0.995 921 | 0.090 31 | −75.885 589 | 0.999 996 | 0.002 50 | −75.885 606 |
3*Eq | 12–1 | −75.864 577 | 0.990 644 | 0.136 79 | −75.876 669 | 0.999 997 | 0.002 42 | −75.876 684 |
aIn the 7-orbital space, FCI, GUCCSD and UCCSD have 441, 1638, and 140 variational parameters, respectively. In the 9-orbital space, FCI, GUCCSD and UCCSD have 15 876, 4572, and 560 variational parameters, respectively. In the 12-orbital space excluding 1 core orbital, FCI, GUCCSD and UCCSD have 108900, 10 340, and 1092 variational parameters, respectively. Energies are given in Eh , the FCI overlap is defined as |⟨ΨVQE|ΨFCI⟩| and the FCI difference norm is computed as , where ΨVQE has been phase-aligned to have a positive overlap with ΨFCI.
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 N2. 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.
Table 3. Energies obtained with various methods using both bare and DUCC Hamiltonians for the twisting of ethylene.
Hamiltonian: | HF | FCI | FCI | FCI | FCI | FCI | VQE | VQE | FCI | FCI | FCI | CCSD | CCSD | CCSDTQ |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Bare | Bare | Bare | Bare | Bare | Bare | DUCC | DUCC | DUCC | DUCC | DUCC | Bare | Bare | Bare | |
Orbitals | [9] | [9*] | [15] | [15*] | [13*] | [9*] | [13*] | [9*] | [15*] | [13*] | [48*] | [48] | [48*] | |
0 | −78.0397 | −78.0593 | −78.0593 | −78.0851 | −78.0850 | −78.0738 | −78.3063 | −78.3099 | −78.3063 | −78.3109 | −78.3099 | −78.3452 | −78.3499 | −78.3563 |
10 | −78.0375 | −78.0569 | −78.0569 | −78.0830 | −78.0829 | −78.0716 | −78.3042 | −78.3078 | −78.3042 | −78.3088 | −78.3078 | −78.3432 | −78.3479 | −78.3543 |
20 | −78.0307 | −78.0500 | −78.0500 | −78.0767 | −78.0766 | −78.0653 | −78.2978 | −78.3016 | −78.2978 | −78.3026 | −78.3016 | −78.3372 | −78.3419 | −78.3485 |
30 | −78.0196 | −78.0391 | −78.0390 | −78.0662 | −78.0662 | −78.0549 | −78.2874 | −78.2913 | −78.2874 | −78.2923 | −78.2913 | −78.3272 | −78.3319 | −78.3388 |
40 | −78.0040 | −78.0244 | −78.0244 | −78.0519 | −78.0518 | −78.0406 | −78.2729 | −78.2770 | −78.2729 | −78.2780 | −78.2770 | −78.3134 | −78.3181 | −78.3254 |
50 | −77.9841 | −78.0063 | −78.0063 | −78.0340 | −78.0339 | −78.0227 | −78.2545 | −78.2589 | −78.2545 | −78.2598 | −78.2589 | −78.2959 | −78.3005 | −78.3085 |
60 | −77.9599 | −77.9854 | −77.9854 | −78.0130 | −78.0129 | −78.0017 | −78.2324 | −78.2371 | −78.2324 | −78.2381 | −78.2371 | −78.2750 | −78.2796 | −78.2888 |
70 | −77.9317 | −77.9629 | −77.9629 | −77.9902 | −77.9903 | −77.9791 | −78.2071 | −78.2123 | −78.2071 | −78.2133 | −78.2123 | −78.2512 | −78.2559 | −78.2675 |
80 | −77.8996 | −77.9421 | −77.9421 | −77.9688 | −77.9687 | −77.9579 | −78.1794 | −78.1852 | −78.1794 | −78.1865 | −78.1852 | −78.2259 | −78.2305 | −78.2497 |
90 | −77.8640 | −77.9353 | −77.9308 | −77.9599 | −77.9599 | −77.9493 | −78.1515 | −78.1584 | −78.1515 | −78.1600 | −78.1584 | −78.2020 | −78.2066 | −78.2387 |
aNumbers with an asterisk indicate that the atomic core orbitals were frozen in the corresponding calculation.
Download figure:
Standard image High-resolution imageThe bond breaking of N2 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.
Download figure:
Standard image High-resolution imageIn tables 4 and 5 we present VQE results for N2 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 MS = 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−8 a.u. or maximum gradient component 1 × 10−5 a.u.
Table 4. Energies obtained with the VQE method using a bare Hamiltonian in the STO-3G basis set for N2 and CAS(10, 8) at r = 1 Å.
Ket | Initial amplitudes | GCCSD-FCI overlap | Iterations | Final energy | GCCSD-FCI Δ norm |
---|---|---|---|---|---|
FCI (reference) | −107.548 966 77 | ||||
HF | 0 | 0.999 985 11 | 208 | −107.548 889 70 | 0.00545666 |
HF | MP2 | 0.999 985 11 | 217 | −107.548 889 70 | 0.00545533 |
HF | MPS8 | 0.999 985 12 | 211 | −107.548 889 75 | 0.00545399 |
HF | MPS16 | 0.999 985 12 | 214 | −107.548 889 74 | 0.00545419 |
HF | MPS32 | 0.999 985 12 | 219 | −107.548 889 66 | 0.00545513 |
MPS8 | 0 | 0.999 992 10 | 857 | −107.548 906 41 | 0.003 973 16 |
MPS8 | MP2 | 0.999 992 30 | 831 | −107.548 907 60 | 0.003 922 53 |
MPS8 | MPS8 | 0.999 992 82 | 831 | −107.548 910 98 | 0.003 788 11 |
MPS16 | 0 | 0.999 996 10 | 814 | −107.548 933 61 | 0.002 789 31 |
MPS16 | MP2 | 0.999 994 57 | 711 | −107.548 922 60 | 0.003 293 18 |
MPS16 | MPS16 | 0.999 995 69 | 814 | −107.548 930 83 | 0.002 935 01 |
MPS32 | 0 | 0.996 835 91 | 256 | −107.526 601 87 | 0.079 549 80 |
MPS32 | MP2 | 0.893 801 38 | 504 | −106.715 349 43 | 0.460 865 73 |
MPS32 | MPS32 | 0.999 993 11 | 1150 | −107.548 913 19 | 0.003 711 58 |
Table 5. Energies obtained with the VQE method using a bare Hamiltonian in the STO-3G basis set for N2 and CAS(10, 8) at r = 2.5Å
Ket | Initial amplitudes | GCCSD-FCI overlap | Iterations | Final energy | GCCSD-FCI Δ norm |
---|---|---|---|---|---|
FCI (reference) | −107.440 409 82 | ||||
HF | 0 | 0.965 976 82 | 315 | −107.440 204 46 | 0.260 856 97 |
HF | MP2 | 0.965 905 47 | 271 | −107.440 200 84 | 0.261 130 31 |
HF | MPS8 | 0.965 708 57 | 351 | −107.440 206 78 | 0.261 883 27 |
HF | MPS16 | 0.038 534 23 | 313 | −105.685 445 05 | 1.386 698 06 |
HF | MPS32 | 0.967 042 49 | 266 | −107.440 200 92 | 0.256 739 18 |
MPS8 | 0 | 0.999 799 56 | 507 | −107.440 357 31 | 0.020 021 57 |
MPS8 | MP2 | 0.988 668 06 | 113 | −107.434 149 46 | 0.150 545 21 |
MPS8 | MPS8 | 0.999 956 91 | 486 | −107.440 347 73 | 0.009 283 31 |
MPS16 | 0 | 0.999 737 94 | 1191 | −107.440 367 18 | 0.022 893 64 |
MPS16 | MP2 | 0.999 495 04 | 1023 | −107.440 361 40 | 0.031 779 10 |
MPS16 | MPS16 | 0.999 324 60 | 1139 | −107.440 361 68 | 0.036 753 16 |
MPS32 | 0 | 0.999 998 98 | 31 | −107.440 407 35 | 0.001 424 49 |
MPS32 | MP2 | 0.999 997 94 | 115 | −107.440 406 25 | 0.002 025 58 |
MPS32 | MPS32 | 0.999 999 02 | 122 | −107.440 408 06 | 0.001 397 29 |
aThe results are partially converged (unable to achieve full convergence under calculation limits).
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.
5. Conclusions
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 N2, H2O, and C2H4. Last but not least, we have shown on the example of N2, 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.
Table 6. Energies obtained with various methods using both bare and DUCC Hamiltonians for the dissociation of N2 from 0.8 to 3.5 Å.
Hamiltonian: | HF | FCI | HF | FCI | VQE(GUCCSD) | VQE(UCCSD) | CCSD | CCSDT | CCSDTQ |
---|---|---|---|---|---|---|---|---|---|
Bare | Bare | DUCC | DUCC | DUCC | DUCC | Bare | Bare | Bare | |
Orbitals: | [10] | [10] | [10] | [10] | [10] | [60] | [60] | [60] | |
0.8 | −108.4988 | −108.5177 | −108.8202 | −108.8372 | −108.8372 | −108.8370 | −108.8556 | −108.8663 | −108.8667 |
0.9 | −108.8399 | −108.8719 | −109.1621 | −109.1892 | −109.1891 | 109.1889 | −109.2081 | −109.2209 | −109.2215 |
1.0 | −108.9680 | −109.0165 | −109.2916 | −109.3305 | −109.3304 | −109.3301 | −109.3499 | −109.3653 | −109.3664 |
1.1 | −108.9830 | −109.0504 | −109.3089 | −109.3611 | −109.3611 | −109.3605 | −109.3809 | −109.3996 | −109.4012 |
1.2 | −108.9396 | −109.0287 | −109.2682 | −109.3361 | −109.3360 | −109.3353 | −109.3561 | −109.3786 | −109.3810 |
1.3 | −108.8684 | −108.9846 | −109.1961 | −109.2860 | −109.2858 | −109.2849 | −109.3056 | −109.3325 | −109.3362 |
1.4 | −108.7864 | −108.9383 | −109.1043 | −109.2272 | −109.2270 | −109.2257 | −109.2461 | −109.2780 | −109.2834 |
1.6 | −108.6221 | −108.8572 | −108.9106 | −109.1146 | −109.1144 | −109.1119 | −109.1305 | −109.1757 | −109.1848 |
1.8 | −108.4775 | −108.7952 | −108.7445 | −109.0279 | −109.0277 | −109.0244 | −109.0390 | −109.1099 | −109.1162 |
2.0 | −108.3575 | −108.7510 | −108.6026 | −108.9710 | −108.9708 | −108.9685 | −108.9814 | −109.1128 | −109.0826 |
2.2 | −108.2602 | −108.7500 | −108.4710 | −108.9369 | −108.9368 | −108.9351 | −108.9687 | −109.1588 | |
2.4 | −108.1818 | −108.7465 | −108.3476 | −108.9158 | −108.9157 | −108.9128 | −108.9927 | −109.1899 | |
2.6 | −108.1186 | −108.7453 | −108.2493 | −108.9047 | −108.9046 | −108.9002 | −109.0171 | −109.2066 | |
2.8 | −108.0675 | −108.7446 | −108.1762 | −108.9001 | −108.9000 | −108.8943 | −109.0329 | −109.2156 | |
3.0 | −108.0259 | −108.7435 | −108.1218 | −108.8987 | −108.8986 | −108.8921 | −109.0427 | −109.2204 | |
3.25 | −107.9842 | −108.7438 | −108.0725 | −108.8991 | −108.8990 | −108.8912 | −109.0501 | −109.2236 | |
3.5 | −107.9511 | −108.7435 | −108.0374 | −108.9003 | −108.9002 | −108.8923 | −109.0546 | −109.2253 |
aIn the employed orbital space, the number of variational parameters is 7020 in the GUCCSD Ansatz, 609 in UCCSD, and 14 400 in FCI.
Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.
Acknowledgment
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].