An adaptive variational algorithm for exact molecular simulations on a quantum computer. — TEI-tagged view

Harper R. GrimsleypersName, Sophia E. EconomoupersName, Edwin BarnespersName, Nicholas J. MayhallpersName (2019date). “An adaptive variational algorithm for exact molecular simulations on a quantum computer..” Nature Communicationstitle·j.

Object of study — variational quantum algorithms for exact molecular electronic-structure simulation

Research question — how to avoid the approximation built into a pre-selected VQEterm·mth ansatzterm·key

Contribution — ADAPT-VQEterm·mth, which grows the ansatzterm·key one operator at a time as dictated by the molecule

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

An adaptive variational algorithm for exact molecular simulations on a quantum computer.

Abstract Quantum simulation of chemical systems is one of the most promising near-term applications of quantum computers. The variational quantum eigensolverterm·mth, a leading algorithm for molecular simulations on quantum hardware, has a serious limitation in that it typically relies on a pre-selected wavefunction ansatzterm·mth that results in approximate wavefunctionsterm·dom and energies. Here we present an arbitrarily accurate variational algorithm that, instead of fixing an ansatzterm·key upfront, grows it systematically one operator at a time in a way dictated by the molecule being simulated. This generates an ansatzterm·key with a small number of parameters, leading to shallow-depth circuits. We present numerical simulations, including for a prototypical strongly correlatedterm·key molecule, which show that our algorithm performs much better than a unitary coupled cluster approach, in terms of both circuit depthterm·key and chemical accuracyterm·key. Our results highlight the potential of our adaptive algorithm for exact simulations with present-day and near-term quantum hardware.

Introduction Anticipation that a useful quantum computer will be realized in the near future has motivated intense research into developing quantum algorithms which can potentially make progress on classically intractable computational problems. While many research areas expect to see transformative change with the development of such quantum devices, computational chemistry is poised to be among the first domains to significantly benefit from such new technologies. Due to the exponential growth in the size of the Hilbert space with increasing orbitals, a quantum computer with tens of qubits could potentially surpass classical algorithms 1 – 3 . Achieving such a capability depends not only on the quality of the qubits, but also critically on the efficiency of the algorithms.

The phase estimation algorithmterm·mth (PEAterm·mth) 4 was the first algorithm proposed for simulating electronic structure problems on a quantum computer 1 , 5 . PEAterm·mth provides a path for obtaining the exact ground state electronic energy for a molecule by evolving in time a quantum state with significant overlap with the ground state using the molecular Hamiltonian of interest. Due to the very long circuit depthsterm·key and complex quantum gates required by PEAterm·mth, the coherence times needed to simulate interesting electronic states would exceed the coherence times available on any existing or near-term quantum device. Improvements to PEAterm·mth still require significant resources and experimental demonstrations to date only involve a few qubits 6 – 8 .

In order to reduce the significant hardware demands required by PEAterm·mth and exploit the capabilities of noisy intermediate-scale quantumterm·dom (NISQterm·dom) devices 9 , the variational quantum eigensolverterm·mth (VQEterm·mth) algorithm was proposed and demonstrated using photonic qubits by Peruzzo et al. 10 . This was followed by several theoretical studies on VQEterm·mth 7 , 11 – 17 and demonstrations on other hardware such as superconducting qubits 7 , 16 , 18 and trapped ions 19 , 20 . Other approaches have been pursued as well, including methods for adiabatic quantum computation 21 and quantum machine learning 22 .

VQEterm·mth is a hybrid quantum-classical algorithmterm·mth, because the computational work is shared between classical and quantum hardware. VQEterm·mth starts with an assumption about the form of the target wavefunction. Based on this form, an ansatzterm·key with several tunable parameters is constructed, and a quantum circuit capable of producing this ansatzterm·key is designed. The ansatzterm·key parameters are variationally adjusted until they minimize the expectation value of the molecular Hamiltonian. Classical hardware is used to precompute all the Hamiltonian terms and to update the parameters during the circuit optimization. The quantum hardware is only used to prepare a state (defined by its current set of ansatzterm·key parameter values) and to perform measurements of the various interaction terms in the molecular Hamiltonian, Ĥ = ∑ i g i ô i . Because the individual operator terms, ô i , generally do not commute, the state preparation has to be repeated multiple times, until all the individual operators have been measured enough times to get sufficient statistics on their mean value. Details on all these steps can be found in ref. 12 .

Compared to PEAterm·mth, VQEterm·mth is much more suitable for NISQterm·dom devices, trading in the long circuit depthsterm·key for shorter state preparation circuits, at the expense of a much higher number of measurements. Although VQEterm·mth has been demonstrated to be more efficient and error-tolerant 7 , 12 , 16 , this comes with the compromise that the ansatzterm·key generally only allows one to obtain approximations to the ground state. Because the choice of ansatzterm·key determines the variational flexibility of the trial state, the quality of a VQEterm·mth simulation is only as good as the ansatzterm·key.

Several approaches have been explored with the goal of creating a compact ansatzterm·key which provides high accuracy with few parameters and shallow circuits. The first ansatzterm·key explored 10 was based on the unitary variant of coupled cluster theoryterm·mth truncated at single and double excitations (UCCSDterm·mth), inspired by early efforts in computational chemistry to improve coupled cluster theory 23 – 26 . In UCCSDterm·mth, trial states are generated by applying to a reference state a unitary operator in the form of an exponential of a sum of single and double fermion operators with their coefficients taken as free parameters. More recent proposals based on UCCSDterm·mth include the unitary Bogoliubov coupled cluster theory which takes a generalized Hartree–Fock (HF) state as the reference 27 and the k -UpCCGSD approach of Lee et al. 17 which uses k products of unitary paired generalized doubles excitations, along with the full set of generalized single excitations. The k -UpCCGSD approach builds on early work by Nakatsuji 28 – 31 and Nooijen 32 studying the use of generalized excitation terms in classical quantum chemistry algorithms, but prunes the expansive operator list by restricting the two-particle terms to only paired interactions, which provides a systematic way to converge to FCI without introducing higher excitation rank operators. Ryabinkin et al. 33 recently proposed a coupled cluster-like ansatzterm·key which is constructed directly in the qubit representation with the goal of achieving shallower circuits. While not directly a variation of the UCCSDterm·mth ansatzterm·key itself, ref. 16 developed an approach (termed the quantum subspace expansion) to extract not just the expectation value of Ĥ but all the matrix elements ⟨ I ∣ Ĥ ∣ J ⟩ in a small subspace consisting of single excitations from the trial state. This Hamiltonian matrix is then diagonalized on a classical computer, which reduces the impact of decoherence and gives access to excited states. Even further from the original UCC ansatzterm·key, Kandala et al. 18 have used an alternative ansatzterm·key for their VQEterm·mth experiments based on the native entangling gate in their superconducting qubit device, referred to as a “hardware-efficient ansatzterm·key”. This allows entanglement to be created directly from a device-wide unitary instead of through a more traditional gate decomposition of a fermionic operatorterm·key.

Despite these considerable improvements to the UCCSDterm·mth ansatzterm·key for VQEterm·mth, this remains an approximate approach that works best for systems that are not strongly correlatedterm·key. However, strongly correlatedterm·key systems are the hardest to simulate classically, and this is precisely the motivation for performing simulations using quantum computers. While an exact VQEterm·mth simulation could in principle be performed by adding higher rank excitations to the ansatzterm·key, this would be prohibitively expensive for both the classical subroutines and NISQterm·dom devices. To overcome these challenges, we need to avoid imposing an ad hoc ansatzterm·key and instead allow the system to determine its own compact, quasi-optimal ansatzterm·key.

In this paper, we achieve this by introducing a simple algorithm termed Adaptive Derivative-Assembled Pseudo-Trotter ansatzterm·key Variational Quantum Eigensolverterm·mth (ADAPT-VQEterm·mth). ADAPT-VQEterm·mth determines a quasi-optimal ansatzterm·key with the minimal number of operators for a desired level of accuracy. The key idea is to systematically grow the ansatzterm·key by adding fermionic operatorsterm·key one-at-a-time, such that the maximal amount of correlation energy is recovered at each step. This results in a wavefunction ansatzterm·mth that is discovered by the algorithm, and which cannot be predicted a priori from a traditional excitation-based scheme like UCCSDterm·mth. While intuitive, this approach can also be derived more rigorously as a particular optimization procedure for Full Configuration Interaction (FCI) VQEterm·mth and is more thoroughly discussed in Section 1 of the Supplement . We demonstrate the power of ADAPT-VQEterm·mth through numerical simulations of three molecules of increasing complexity: LiHterm·key, BeH 2 , and H 6 . In each case, we find vastly improved performance compared to UCCSDterm·mth, both in terms of the number of operators needed to form the trial states and in terms of chemical accuracyterm·key. Therefore, we believe that ADAPT-VQEterm·mth is an ideal hybrid algorithm for NISQterm·dom devices.

Results Specification of the adopted notation In order to define the approach, several definitions and notations need to be established. First, molecular orbital indices i and j denote occupied orbitals, a and b denote virtual orbitals, and p , q , r , and s denote arbitrary molecular orbitals. In coupled cluster theory, in particular CCSD, an expansion based on the HF state | ψ HF 〉 is created by using an exponential ansatzterm·key involving single and double excitation operators: 1 ψ CCSD = e T ^ 1 + T ^ 2 ψ HF , where the excitation operators are defined as: 2 T ^ 1 = ∑ i a t ^ i a = ∑ i a t i a â a † â i 3 T ^ 2 = ∑ i < j , a , b t ^ i j a b = ∑ i < j , a < b t i j a b â a † â b † â i â j .

For closed shell molecules near equilibrium, CCSD provides a robust ansatzterm·key for molecular simulations. Early efforts to combine size extensivity and variationality were pioneered by Bartlett, Kutzelnigg, and coworkers 23 – 25 . In this context, a unitary variant of coupled cluster theoryterm·mth (UCCSDterm·mth) was defined by replacing the excitation operators with an anti-Hermitian sum of excitation and de-excitation operators: 4 t ^ i j a b → t ^ i j a b - t ^ a b i j = τ ^ i j a b .

Because UCCSDterm·mth is based on a unitary operator, the adjoint is the inverse, and the expectation value of the UCCSDterm·mth wavefunction can be expanded using the Baker–Campbell–Hausdorff (BCH) formula to obtain a normalized Hamiltonian expectation value (Rayleigh quotient) for variational optimization. Unfortunately, the BCH expansion does not truncate at finite order, making UCCSDterm·mth computationally intractable on classical hardware. However, the unitary nature of UCCSDterm·mth is actually a benefit for quantum algorithms as it corresponds to a coherent time evolution, and this was the original motivation for using UCCSDterm·mth in VQEterm·mth 10 .

In addition to a unitary form, CCSD can also be generalized by including excitation operators which immediately annihilate the HF state. These would include excitations from occupied to occupied, virtual to virtual, etc. Generalized excitations or interactions of this form have been considered previously, and have been used in the context of VQEterm·mth recently by Lee and coworkers 17 . In this case the cluster operators are further generalized to remove the HF-based subspace restriction: τ ^ i j a b → τ ^ p q r s , where p , q , r , and s refer to any arbitrary orbital.

Although UCCSDterm·mth is perhaps a natural ansatzterm·key for VQEterm·mth, it cannot be implemented directly as written or as explored previously in the quantum chemistry context. Because the gate model of quantum computation is realistically bound to using gates acting on only a few qubits at a time, the UCCSDterm·mth operator must be broken up into a time-ordered sequence of few (one or two) particle operators. This is achieved by using a Trotter expansion of a matrix exponential 34 , 5 e A + B = lim n → ∞ e ( A ∕ n ) e ( B ∕ n ) n .

Because the generalized single and double excitation operators do not commute, the use of a truncated Trotter expansion represents an approximation to the underlying UCCSDterm·mth ansatzterm·key, and recent work has shown clearly that this does not strongly affect the results because the variational flexibility is sufficient to absorb this error 7 , and that even a single Trotter number ( n = 1) is sufficient to reproduce the results of UCCSDterm·mth. As a result, a unitary, generalized, Trotterized ansatzterm·key becomes: 6 ψ tUCC = ∏ s ∈ { p q } e t ^ s ∏ d ∈ { p q r s } e t ^ d ψ HF , where notation is introduced such that the generalized singles index, s , runs over all unique pairs of p , q and the doubles index, d , over unique combinations of p , q , r , s .

ADAPT-VQEterm·mth algorithm The above discussion described the Trotter expansion as an approximation to UCCGSD. However, as recognized previously 7 , 35 , if the parameters are optimized after the Trotterization, this is not so much an approximation to UCC as it is a wholly unique ansatzterm·key. In fact, the exact FCI solution could be obtained by simply going to an n th order Trotterized form of UCCSDterm·mth and allowing the different parameter replicas to vary independently. This is due to the fact that n -body interactions can be described as products of one- and two-body interactions. The exact (FCI) quantum state can thus be represented as an arbitrarily long product of one- and two-body operators, 7 ψ FCI = ∏ k ∞ ∏ p q e τ ^ p q ( k ) ∏ p q r s e τ ^ p q r s ( k ) ψ HF where τ ^ p q r s ( k ) is the k th instance, or “replica”, of the operators in t ^ p q r s - t ^ r s p q . It is important to note that this is not a Trotter approximation to any simple two-body ansatzterm·key, as each replica can assume different parameter values, e.g., τ p q r s ( k ) ≠ τ p q r s ( j ) .

The main goal in this paper is to approximate FCI with arbitrary accuracy using a maximally compact sequence of unitary operators. The basic outline of the algorithm is drawn schematically in Fig. 1 and is as follows: On classical hardware, compute one- and two-electron integrals, and transform the fermionic Hamiltonian into a qubit representation using an appropriate transformation: Jordan–Wigner, Bravyi–Kitaev, etc. This is a standard step in regular VQEterm·mth. Define an “Operator Pool”. This is simply a collection of operator definitions which will be used to construct the ansatzterm·key. For the examples presented in the next section, we consider the set of all unique spin-complemented one- and two-body operators, but one might imagine adding a few three-body or four-body terms as well. Initialize qubits to an appropriate reference state, ideally one with the correct number of electrons. The HF state would be a sensible choice here. Initialize the ansatzterm·key to the identity operator. On a quantum computer, prepare a trial state with the current ansatzterm·key. If multiple quantum computers are available, perform this step on all devices simultaneously. Measure the commutator of the Hamiltonian with each operator in the pool to get the gradient. Repeating this multiple times and averaging gives the gradient of the expectation value of the Hamiltonian with respect to the coefficient of each operator. This can be done in parallel. If the norm of the gradient vector is smaller than some threshold, ε , exit. Identify the operator with the largest gradient and add this single operator to the left end of the ansatzterm·key, with a new variational parameter. Note that this does not “drain” the pool in the sense that choosing an operator does not remove it from the pool so it can be used again later. Perform a VQEterm·mth experiment to re-optimize all parameters in the ansatzterm·key. Go to step 4. Fig. 1 Schematic depiction of the ADAPT-VQEterm·mth algorithm described presented. Since step 1 occurs on classical hardware, it is not included in the illustration. θ → ( n ) is the list of ansatzterm·key parameters at the n th iteration. The number of parameters, l e n ( θ → ( n ) ) , is equal to the number of operators in the ansatzterm·key. “Operator Pool” refers to the collection of operators which are used to grow the ansatzterm·key one-at-a-time. Each τ p q represents a generalized single or double excitation, and these operators are then spin-complemented. The orbital indices refer to spatial orbitals, and the overbar indicates β spin. Orbital indices without overbars have α spin. Note that growing the ansatzterm·key does not drain the pool, and so operators can show up multiple times if selected by the algorithm As described above and illustrated in Fig. 1 , each iteration starts as a series of uncoupled experiments to obtain the parameter gradients via measurements of operator commutators (the gradient expression in step 5 is derived in section IB of the Supplemental Information ). The purpose of these gradient measurements is to determine the best operator with which to grow the ansatzterm·key, as the operator with the largest gradient is likely to recover the most correlation energy in the subsequent VQEterm·mth minimization. This process is continued iteratively, until a convergence threshold is met. In the classical numerical examples presented below, we chose to consider the L 2 norm of the gradient vector to determine convergence. This is just one possibility, and alternative convergence indicators could be used instead in step 6. At convergence, the ADAPT-VQEterm·mth algorithm obtains the following ansatzterm·key: 8 ψ ADAPT ( ε ) = e τ ^ N e τ ^ N - 1 ⋯ e τ ^ 2 e τ ^ 1 ψ HF where the identity of each τ ^ i is determined by the algorithm.

The re-optimization subroutine in step 8 can be implemented on either a classical or quantum processor using any of the gradient- or non-gradient-based optimization routines that have been proposed or demonstrated for VQEterm·mth 12 , 15 , 16 , 18 . Note that this subroutine is distinct from the gradient computed in step 5 of the algorithm. Additional possible modifications to the algorithm are mentioned in the “Discussion” section.

The evaluation of all the gradient terms could in principle be achieved in a NISQ-friendly, highly parallel manner with a large number of uncoupled quantum computers all tasked with preparing the same state and measuring a different operator. This is the same potential for parallelization that the underlying VQEterm·mth subroutine has. Just as with the original motivation for VQEterm·mth, ADAPT-VQEterm·mth decreases the circuit depthterm·key at the expense of a larger number of measurements. In our case a sequence of VQEterm·mth experiments is performed, with the most resource-demanding experimental steps happening at the end. This constitutes a rather large prefactor which would scale with the size of the system, but the crucial advantage is controllability over the ansatzterm·key accuracy (in principle approaching FCI). Because the number of non-zero parameters equals the number of iterations, in order to discover an ansatzterm·key for a large system, an equal number of VQEterm·mth re-optimizations will need to be performed. One strategy to minimize this prefactor could simply be to add a few operators at a time.

Determining resource requirements for adaptive procedures is rather difficult. The classical resources are not expected to be significant in the foreseeable future. However, as quantum technology progresses toward deeper circuits, the parameter manipulation and updating on a classical computer could become costly. However, we expect the dependence between parameters at the beginning and end of the ADAPT-VQEterm·mth circuit to decay with circuit depthterm·key, such that one could imagine freezing the early parameters after a certain number of iterations. This would possibly establish an approach for FCI with only a polynomial number of variables, completely avoiding any exponential cost for the classical hardware.

Molecular dissociation simulation results In this section, we explore the convergence properties of the ADAPT-VQEterm·mth algorithm with a few small molecular systems, LiHterm·key, BeH 2 , and linear H 6 . The former two molecules have been simulated using quantum hardware 18 , 20 . H 6 is included as a prototypical strongly correlatedterm·key molecule, which allows us to test the ADAPT-VQEterm·mth approach for systems which are not well described with unitary coupled cluster.

In order to perform the simulations, an in-house code was written, using Psi4 36 , 37 for the integral calculation (via the OpenFermion-Psi4 38 interface) and OpenFermion was used for the Jordan–Wigner operator transformation. All calculations used the Broyden–Fletcher–Goldfarb–Shannon (BFGS) minimization implemented within Scipy 39 . To classically simulate the re-optimization subroutine in step 8 of ADAPT-VQEterm·mth, we could use a standard numerical gradient method. However, in order to improve the efficiency and allow precise gradients for tight convergence, we derived and implemented an efficient analytic gradient function, which is detailed in Section D of the Supplement . By reusing intermediates between individual parameter gradients, this algorithm obtains the full gradient vector for all parameters at a cost which is only roughly 2× that of the base energy evaluation. A table with timing data is included in Table 1 of the Supplement . Table 1 Average errors across the PES scan for the different methods assessed. Units in kcal/mol UCCSDterm·mth ADAPT( ε 1 ) ADAPT( ε 2 ) ADAPT( ε 3 ) LiHterm·key 0.0480 0.3000 0.0058 0.0002 BeH 2 2.2384 0.8023 0.0907 0.0041 H 6 3.7387 4.5297 0.3023 0.0047 As discussed in the previous section, the ADAPT ansatzterm·key uses a convergence threshold to determine when the calculation should terminate (step 6). Here we use the norm of the gradient vector and compare it to threshold ε m , which we define as 9 ε m = 1 0 - m .

For example, an ADAPT-VQEterm·mth calculation where the norm of the operator pool gradient is converged to less than 0.001 would be denoted as ADAPT( ε 3 ). In what follows we present numerical results for bond-dissociation curves for LiHterm·key, BeH 2 , and H 6 for three different choices of the threshold ( m = 1, 2, 3). We also investigate alternate protocols for the ansatzterm·key growth and demonstrate the superiority of the ADAPT ansatzterm·key.

Here, we study the LiHterm·key bond dissociation computed using several methods, including FCI, UCCSDterm·mth (un-Trotterized), HF, ADAPT( ε 1 ), ADAPT( ε 2 ), and ADAPT( ε 3 ), all with the STO-3G basis set. In this basis set, LiHterm·key has 6 spatial orbitals and a Hilbert space of dimension 4096. By starting with the HF state with two α (spin-up) and two β (spin-down) electrons and using only number conserving operators, the relevant subspace to explore has a dimension of 6 2 ⋅ 6 2 = 225 . In this basis, the occupied orbitals are {1, 2}, and the virtual orbitals are {3, 4, 5, 6}.

The bond dissociation curves are shown in Fig. 2a , where all the curves, with the exception of HF, cannot really be distinguished on this scale. However, as shown in Fig. 2b , when the FCI energy is subtracted and the scale is adjusted, significant differences become evident. Shading is used to indicate chemical accuracyterm·key, which is achieved in all cases other than HF. LiHterm·key has only a single pair of electrons (a σ bond) breaking along the dissociation coordinate, and UCCSDterm·mth exhibits chemical accuracyterm·key throughout the curve. While ADAPT( ε 1 ) is not as accurate as UCCSDterm·mth, ADAPT( ε 2 ) is comparable to UCCSDterm·mth at short bond distances and comfortably outperforms it at longer distances. This is also evident in Table 1 , where the average error across the potential energy surface (PES) is shown. Remarkably, ADAPT( ε 3 ) outperforms UCCSDterm·mth throughout the whole curve by at least an order of magnitude and in some cases up to four orders of magnitude. Fig. 2 Dissociation curves for LiHterm·key, BeH 2 , and H 6 . Potential energy as a function of nuclear coordinate, Hartree units ( a , d , g ). Absolute energy differences from FCI, kcal/mol units ( b , e , h ). Shaded blue region represents area within “chemical accuracyterm·key” as 1 kcal/mol. Number of variational operators in associated ansatzterm·key ( c , f , i ). Notation: ε indicates gradient norm threshold used such that ε m = 10 − m . In all curves, the FCI curve lies directly underneath the ADAPT( ε 3 ) curve, and so is not visible Even more impressive is how few parameters are needed to achieve this level of accuracy. As shown in Fig. 2c , in all three cases and for all bond distances, ADAPT is much more compact than UCCSDterm·mth. UCCSDterm·mth has 92 parameters, which can be reduced to 64 by combining spin-complements. In all three ADAPT calculations, fewer than half of the parameters are needed compared to UCCSDterm·mth. Although UCCSDterm·mth is noticeably more accurate than the simplest ADAPT calculation with a gradient norm threshold of 0.1, the ADAPT( ε 1 ) ansatzterm·key is incredibly compact, consisting of fewer than 10 parameters across the curve. For example, the ADAPT( ε 1 ) ansatzterm·key for LiHterm·key at bond distance 2.39 Å is 10 ψ ADAPT ( ε 1 ) = e τ ^ 12 16 e τ ^ 2 2 ¯ 5 5 ¯ e τ ^ 2 2 ¯ 4 4 ¯ e τ ^ 12 13 e τ ^ 1 1 ¯ 3 3 ¯ e τ ^ 2 2 ¯ 3 3 ¯ e τ ^ 2 2 ¯ 3 6 ¯ e τ ^ 2 2 ¯ 6 6 ¯ ψ HF , which includes a mixture of both double excitations and correlated single excitations n ^ j â a † â i . The indices denote spatial orbitals, overbar on an index denotes β spin, and spin-complemented interactions are implied. For example τ ^ 01 06 is really τ ^ 01 06 + τ ^ 0 ¯ 1 ¯ 0 ¯ 6 ¯ . An interesting feature of the ansatzterm·key returned by ADAPT-VQEterm·mth, Eq. ( 10 ), is that the HOMO–LUMO double excitation e τ ^ 2 2 ¯ 3 3 ¯ is not the first operator, but instead the third. This is different from what one might expect if classical MP2 or CCSD amplitudes were used to order the ansatzterm·key. The reason is that in choosing the next operator no state energy information is used, for instance in the form of a denominator penalizing high energy terms. Interestingly, at convergence it is not the HOMO–LUMO term or the first operator with the largest amplitude, but rather the second operator, e τ ^ 2 2 ¯ 3 6 ¯ .

In Fig. 2d–f , the dissociation curves for BeH 2 are shown. In the STO-3G basis, BeH 2 has 7 spatial orbitals, for a total Hilbert space dimension of 16,384, and a neutral molecule subspace of dimension 7 3 ⋅ 7 3 = 1225 . Unlike with LiHterm·key, UCCSDterm·mth does not provide chemically accurate results across the full PES. UCCSDterm·mth and ADAPT( ε 1 ) are comparable at smaller bond distances. Beyond ~3 Å, they both go above 1 kcal/mol in absolute error. However, still with a small fraction of the number of parameters in UCCSDterm·mth, both ADAPT( ε 2 ) and ADAPT( ε 3 ) provide nearly exact results, with average deviations from FCI listed in Table 1 .

Now we move our focus to the H 6 data. At bond-breaking, the previous two molecules involved strong correlation between only two and four electrons, respectively. In order to evaluate the ability of ADAPT-VQEterm·mth to converge to FCI in the presence of much stronger correlations, we have computed the simultaneous stretching of H 6 , with the results presented in Fig. 2g–i .

The complexity of this strongly correlatedterm·key system is reflected in two obvious ways: (1) the failure of UCCSDterm·mth to achieve chemical accuracyterm·key across the curve in Fig. 2h , and (2) the increased number of parameters selected in the ADAPT calculations in Fig. 2i . Despite being strongly correlatedterm·key, such that higher excitation rank operators should be needed, both ADAPT( ε 2 ) and ADAPT( ε 3 ) provide accurate results with only one- and two-body operators. Moreover, in the case of ADAPT( ε 2 ) this is achieved with fewer operators than UCCSDterm·mth for most bond distances. ADAPT( ε 3 ) also uses fewer parameters than UCCSDterm·mth up to the distance where UCCSDterm·mth fails to reach chemical accuracyterm·key.

Because the algorithm is adaptive, during the course of a chemical event (bond breaking, isomerization, etc.) the number of parameters can change abruptly, leading to discontinuous potential energy curves. Two notable examples of this can be seen in Fig. 2h , first at R(H-H) = 1.8 Å where ADAPT( ε 1 ) experiences a large jump in energy, and second at 2.5 Å where ADAPT( ε 2 ) increases in energy. Figure 2i shows that these energy jumps correspond to sudden drops in parameter counts.

The cause of the discontinuities in the H 6 data can be explained from the convergence data provided in the Supplement (see Supplement Fig. 1 ). For larger bond lengths, as additional operators are added to the ansatzterm·key, the energy flattens out before dropping substantially again. If the convergence criterion is too lenient, then the ADAPT-VQEterm·mth optimization will abort at such “false gradient troughs”. In the ADAPT( ε 2 ) data of Fig. 2h, i , the jump in energy error and drop in parameter number, respectively, are caused by the 2.5 Å optimization aborting at a false gradient trough, while the optimizations at other bond lengths do not. Of course, if a tighter threshold is used (such as 0.001), the ADAPT-VQEterm·mth algorithm does not prematurely abort, and ultimately yields high-accuracy results, even for this strongly correlatedterm·key system. More sophisticated convergence checks in step 6 might avoid these situations and will be one focus of future work.

Dependence of convergence on operator ordering To demonstrate the importance of the gradient-based operator ordering chosen by ADAPT-VQEterm·mth, we compare it to a few alternate procedures for growing the ansatzterm·key: (a) Random (ijab): Randomly select from a pool of τ i j a b , where the indices are restricted to those which do not annihilate the HF reference state. (b) Random (pqrs): Randomly select from a pool of τ p q r s , where the indices are not restricted. (c) Lexical (ijab): Select from an ordered pool of τ i j a b , where the indices are restricted to those which do not annihilate the HF reference state. (d) Lexical (pqrs): Select from an ordered pool of τ p q r s , where the indices are not restricted.

In Fig. 3 , we show the convergence of each of these orderings and compare them to ADAPT using BeH 2 as a typical example. What stands out is that the ADAPT ansatzterm·key converges dramatically faster than the other four cases considered. While the two random-growth ansätze converge relatively similar to each other regardless of whether restricted indices are used or not, the lexically ordered ansatzterm·key shows a clear distinction between the restricted index (singles and doubles) and un-restricted index (generalized singles and doubles) ordering. This is due to the fact that the first operators in the ansatzterm·key involve creation operators on the occupied orbitals, and these do not contribute until the wavefunction has become entangled. The un-Trotterized UCCSDterm·mth result is also marked for reference. Overall, the data in Fig. 3 demonstrate that an iterative gradient minimization algorithm yields a highly compact ansatzterm·key for a given state. Fig. 3 Convergence of the energy as a function of parameter count for BeH 2 with a Be–H distance of 2.39 Å. The black dot indicates the UCCSDterm·mth point Discussion An obvious metric for evaluating the performance of any simulation algorithm can be simply described as some accuracy measure vs. some cost measure. While the accuracy measure in a simulation is often easy to define, the cost measure is more nuanced. For variational quantum simulations, there are two factors which largely determine the overall cost: circuit depthterm·key and number of measurements (or shot count). Shot count is important as it determines the time to solution. It is possible that due to the sheer number of measurements, a particular quantum simulation becomes intractable. However, for NISQterm·dom devices in which coherence times (and thus number of gates) are limited, circuit depthterm·key is usually the most critical cost metric, as it determines whether or not a simulation can occur at all. By taking circuit depthterm·key as the most important cost metric to address, the original VQEterm·mth has been successful by minimizing circuit depthterm·key at the cost of increased number of measurements. Similar to the original VQEterm·mth, our new ADAPT-VQEterm·mth algorithm seeks to further minimize the circuit depthterm·key with an increased number of measurements.

In this direction, the data clearly demonstrates that ADAPT-VQEterm·mth succeeds in creating a more compact and accurate wavefunction ansatzterm·mth than UCCSDterm·mth. The algorithm achieves this by systematically identifying the optimal set and ordering of operators to use in the wavefunction ansatzterm·mth for a given problem. The efficiency of ADAPT-VQEterm·mth makes it very promising for quantum chemistry simulations on NISQterm·dom devices, where circuit depthterm·key limitations remain a significant challenge.

In terms of shot-count, ADAPT-VQEterm·mth will likely have an increased number of measurements compared to UCCSD-based VQEterm·mth due to the necessary gradient measurements. However, this is perhaps an easier problem to address (compared to circuit depthterm·key) as the individual runs can in principle occur simultaneously if several devices exist. Further, the shot count also depends on the number of iterations required for the classical optimization of the ansatzterm·key parameters. For strongly correlatedterm·key systems where perturbation theory fails, the existing approach of using classical MP2 amplitudes to initialize the UCCSDterm·mth parameters 12 is not likely to provide much improvement in the UCCSD-based VQEterm·mth. Alternatively, each iteration of ADAPT-VQEterm·mth only adds a single new parameter, with the previously optimized parameters already being initialized to rather sensible values. This might ultimately decrease the number of iterations needed for the VQEterm·mth subroutine in ADAPT-VQEterm·mth, thus decreasing the shot count (although this is not likely to fully compensate for the large number of measurements for the gradient). As hardware capabilities continue to increase, in terms of both the size and number of quantum processors available, ADAPT-VQEterm·mth will offer an ideal quantum-parallel approach to performing nontrivial quantum chemistry simulations. We therefore expect this algorithm to have a strong impact on these efforts in the near term.

As the name suggests, ADAPT-VQEterm·mth could be classified as one member of a family of adaptive-basis strategies that has seen success in constructing compact many-electron wavefunctions 40 – 46 and single-electron wavefunctions 47 – 53 , or as a relative of methods using sequential transformations which have been explored in the context of multireference coupled-cluster theory 54 , 55 . Of these, the ADAPT ansatzterm·key is perhaps most closely related to the @-CC method of Lyakh and Bartlett 47 , in which a compact set of cluster operators is iteratively determined to describe the state of interest on a classical computer. Our approach is distinct in that it is not only designed for a quantum computer implementation, but also defined for a different wavefunction form (product of unitary operators vs. coupled cluster) and a different importance metric (operator gradient of the many-electron state vs. a single electron-defined importance function, see ref. 47 ) for determining new parameters.

An important aspect of ADAPT-VQEterm·mth is that several steps of the algorithm can be implemented in multiple ways, lending it still greater versatility across a wide landscape of problems and suggesting that it should perhaps be thought of as a class of algorithms rather than a specific one. In the “Results” section, we already discussed a few algorithmic options, including different ways to perform the gradient-based parameter update and to determine convergence. We also mentioned the possibility of freezing early parameters at later stages of the algorithm in order to speed up the re-optimization steps. Below, we discuss a few more modifications to explore.

Although the ADAPT-VQEterm·mth algorithm is notably not a perturbative approach, it still has a perturbative flavor in that the suitability of the next iteration’s best operator only involves the interaction of that operator with the Hamiltonian. As such, the algorithm may not be able to recognize the best quadruple excitation (for example) during one update. That being said, the physics described by quadruple excitations is ultimately captured after multiple iterations through the product of at least two two-body interactions. The consequence of this is that convergence will likely not be as fast for strongly correlatedterm·key systems because the algorithm can only “see” two body operators at a time. Because only local knowledge of the FCI energy landscape is used to update the ADAPT-VQEterm·mth ansatzterm·key construction, the “true optimally compact ansatzterm·key” is not guaranteed. As a result, flat energy landscapes (associated with “false gradient troughs”) are possible. Further classical simulations and device implementations are needed to provide better insight into the numerical behavior.

Fortunately, however, multiple strategies can be pursued to address any possible slow convergence issues. One possible approach would be to add a selection of three- or four-body interactions into the operator pool, such that these could be inserted when needed. Alternatively, one might imagine trying to update the ansatzterm·key with two (or more) operators in each iteration, such that the best set of operators is added. The operator pool would still consist of only one- and two-body interactions, but higher-body interactions could be incorporated through products of operators. Even further, one might imagine computing the second derivative and using Hessian matrix elements to identify cooperative effects between operators in the pool. We will explore each of these approaches in future work, with the aim of determining the fastest converging algorithm in different chemical scenarios.

In this paper, we presented ADAPT-VQEterm·mth, a novel variational hybrid quantum-classical algorithmterm·mth designed to achieve exact results at convergence. Unlike typical ansätze, which tends to be ad hoc, our approach is based on an ansatzterm·key that is determined by the system being simulated, and it features a well-defined, built-in convergence criterion. Moreover, the parameter count, and thus the gate depth, is kept to a minimum. A detailed description of the algorithm is given, and numerical examples are provided to demonstrate the performance of the ADAPT method with both weakly and strongly correlatedterm·key systems. Based on these results, we find the ADAPT-VQEterm·mth algorithm to be an operator- and parameter-efficient method capable of high accuracy, with controllable errors, that routinely outperforms UCCSDterm·mth. Its compatibility with classical routines for compiling state preparation circuits and quantum-parallelism should make ADAPT-VQEterm·mth extremely useful for simulations of molecules on both currently available and future quantum computers.

Supplementary Information Supplementary Information