Hardware Efficient Quantum Algorithms for Vibrational Structure Calculations
Pauline J. Ollitrault Alberto Baiardi Markus Reiher Ivano Tavernelli
quant-ph






    Daniil M. Lukin,^†1  Alexander D. White,^†1 Rahul Trivedi,^1 Melissa A. Guidry,^1 Naoya Morioka,^2 Charles Babin,^2 Öney O. Soykal,^3 Jawad Ul-Hassan,^4 Nguyen Tien Son,^4 Takeshi
Ohshima,^5 Praful K. Vasireddy,^6 Mamdouh H.  Nasr,^6 Shuo Sun,^1 Jean-Philippe W. MacLean,^1 Constantin Dory,^1 Emilio A. Nanni,^6 Jörg Wrachtrup,^2 Florian Kaiser,^2 and Jelena Vučković^1
    

=====================================================================================================================================================================================================================================================================================================================================================================================




We introduce a framework for the calculation of ground and excited state energies of bosonic systems suitable for near-term quantum devices and apply it to molecular vibrational anharmonic Hamiltonians.
Our method supports generic reference modal bases and Hamiltonian representations, including the ones that are routinely used in classical vibrational structure calculations.
We test different parametrizations of the vibrational wave function, which can be encoded in quantum hardware, based either on heuristic circuits or on the bosonic Unitary Coupled Cluster Ansatz.
In particular, we define a novel compact heuristic circuit and demonstrate that it provides the best compromise in terms of circuit depth, optimization costs, and accuracy.
We evaluate the requirements, number of qubits and circuit depth, for the calculation of vibrational energies on quantum hardware and compare them with state-of-the-art classical vibrational structure algorithms for molecules with up to seven atoms.




§ INTRODUCTION


Within the Born-Oppenheimer approximation, a molecular wave function is factorized as a product of an electronic part, which is the solution of the electronic Schrödinger equation, and a vibro-rotational one, which is the solution of the nuclear Schrödinger equation in the potential energy surface (PES) generated by sampling the eigenvalues of the electronic Schrödinger equation for different geometries.

The nuclear Schrödinger equation is usually solved in two steps, in analogy with its electronic counterpart. 
A single-particle basis (the basis functions are called, in this case, modals) is obtained either by the harmonic approximation applied to the PES or from a vibrational self-consistent field (VSCF) <cit.> calculation.
Vibrational anharmonic correlations are added a-posteriori with perturbative<cit.> or variational approaches.
The latter include Vibrational Configuration Interaction (VCI) <cit.> and Vibrational Coupled Cluster (VCC) <cit.> for highly-accurate anharmonic energies. 
Unlike perturbation theories, the accuracy of VCI and VCC can be systematically improved, but their applicability is limited to small molecules with up to about 10 atoms due to their unfavorable scaling with system size.
This unfavorable scaling can be tamed by pruning the VCI basis limiting, for instance, the maximum degree of excitation, or with precontraction algorithms.<cit.>
Such simplifications make calculations feasible for systems with up to 15-20 atoms.
Alternatively, the computational cost of VCI can be reduced with non-linear wave function parametrizations. 
This is the case, for example, of the vibrational Density Matrix Renormalization Group (vDMRG) <cit.> which encodes the wave function as a Matrix Product State <cit.>, or of VCC <cit.>.

The emerging development of quantum computers has refreshed the prospect of computing energies of large molecules by leveraging the exponentially-large multi-qubit Hilbert space.
However, current quantum computers based on superconducting qubits technology have limited coherence times (≈100 μs) and sizable gate error rates (≈2×10^-2 for a two-qubit gate), restricting the possible number of operations that can be executed to evolve a quantum state.
Under these limitations, hybrid quantum-classical algorithms are the most promising route to calculate molecular energies on quantum hardware.
In particular, the ground state energy of a general Hamiltonian can be obtained with quantum circuits of relatively low depth (i.e., a small number of operations) with the Variational Quantum Eigensolver (VQE) <cit.>.
The VQE has already been applied in hardware calculation of the electronic ground state of small molecules <cit.> and can also be extended to excited states <cit.>.

Despite the extensive work on the applications of VQE to the solution of the electronic Schrödinger equation, the extension to vibrational has not yet been fully investigated.
Molecular vibrations are described by Bose-Einstein statistics and, therefore, the modal basis must be mapped to the qubits by preserving such symmetry.
Moreover, any many-body expansion of a L-mode PES, contains, in principle, up to L-body coupling terms <cit.>. 
The potential energy operator of the nuclear Schrödinger equation is therefore much more complex than the pairwise Coulomb interaction of the electronic Schrödinger equation.
McArdle et al. <cit.> adapted the VQE to find the ground state of vibrational Hamiltonians of small molecules on a universal quantum computer based on the unitary extension of the VCC theory (UVCC).
In particular, they represent vibrational levels with the so-called compact mapping (the problem of mapping bosonic states to qubits has been discussed in details in a recent work <cit.>).
Two key limitations hinder the application of the theory presented in Ref.  to complex vibrational Hamiltonians.
First, the algorithm can be applied only to vibrational ground states and, therefore, does not allow to access vibrational excitation energies that are key for vibrational spectroscopy. 
Second, it approximates the PES as power series of Cartesian-based normal modes, which relies on harmonic-oscillator eigenfunctions as basis functions.
This is an important limitation in the case of strongly anharmonic molecules, whose PES is represented by highly non-compact Taylor expansions.
For these systems, the VSCF modals will lead instead to more compact VCC and VCI expansions.

In the present paper, we design a general framework for the calculation of vibrational structures with a quantum algorithm.
We introduce the qubit encoding of the vibrational levels based on the generalized second quantization representation <cit.> of the nuclear Schrödinger equation that enables the potential to be expressed as a general n-mode expansion <cit.>.
We then discuss the parametrization of the vibrational wave function introducing a new quantum circuit Ansatz as a compact approximation of UVCC.
We emphasize that, although UCC has mostly been applied in electronic-structure calculations <cit.> and its extension to vibrational quantum computation has hardly been explored <cit.>, its use for the solution of the vibronic problem is very promising.
We discuss how vibrational excited states can be targeted with equation-of-motion (EOM)-based UVCC algorithms <cit.>.
Finally, we also discuss the scaling of the UVCC resources in terms of qubits and gate counts as a function of the molecular size and show the resources necessary to compute the vibrational structure of molecules with up to five atoms.
The proposed framework offers us the possibility to estimate the hardware requirements that will allow to reach quantum advantage over classical vibrational-structure calculations using near-term quantum computers.



§ THEORY




 §.§ Second quantization theories for molecular vibrations


The real-space representation of the Watson Hamiltonian for the L modes of a molecular sytem can be written as


    ℋ_vib(Q_1, …, Q_L) = - 1/2∑_l=1^L∂^2/∂ Q_l^2
        + V(Q_1, …, Q_L)


where Q_l are the harmonic mass-weighted normal coordinates and the Coriolis couplings <cit.> have been neglected.
ℋ_vib must be mapped to an operator that acts on the states of a given set of N_q qubits in order to calculate its eigenfunctions on quantum hardware.
In electronic structure calculations, the mapping is achieved by expressing the non-relativistic electronic Hamiltonian in second quantization, i.e. by projecting it onto the complete set of antisymmetrized occupation number vectors (ONV) generated by a given (finite) set of orbitals.
To encode the vibrational Hamiltonian of Eq. (<ref>) in terms of the second quantization operators defined in Eq. (<ref>), we expand the potential V(Q_1, …, Q_L) with the n-body expansion <cit.>, as follows:


    V(Q_1, …, Q_L) = V_0 + ∑_l=1^L V^[l](Q_l) + ∑_l<m^L V^[l,m](Q_l, Q_m)
        + ∑_l<m<n^L V^[l,m,n](Q_l, Q_m, Q_n) + …


where V_0 is the electronic energy of the reference geometry, the one-mode term V^[l](Q_l) represents the variation of the PES upon change of the l-th normal coordinate from the equilibrium position.
Similarly, the two-body potential V^[l,m](Q_l, Q_m) represents the change in the exact PES upon a simultaneous displacement along the l-th and m-th coordinates <cit.>. 
The exact representation of a PES for an L-mode system requires an L-body expansion.
Often, including terms up to three-body in the L-body expansion is sufficient to obtain an accuracy of about 1 cm^-1.

A representation of Eq. (<ref>) that is suitable to encode on a quantum computer can be obtained with the so-called canonical quantization <cit.> that maps the l-th normal coordinate Q_l and its conjugate momentum P_l to a pair of bosonic creation and annihilation operators (a_l^+ and a_l) defined as


    Q_l    = 1/√(2)( a_l^+ + a_l ) 
    
       P_l    = i/√(2)( a_l^+ - a_l )   ,


where the a_l^+/a_l operators are defined as


    a_l^+ | n_1 ⋯ n_l ⋯ n_L ⟩   = √(n_l+1) | n_1 ⋯ n_l+1 ⋯ n_L⟩
    
       a_l   | n_1 ⋯ n_l ⋯ n_L ⟩   = √(n_l) | n_1 ⋯ n_l-1 ⋯ n_L ⟩.


Each index of the ONV | n_1 ⋯ n_L ⟩ is associated to a mode and n_l is the degree of excitation of the l-th mode.
Different vibrational-structure methods have been derived based on this canonical representation,<cit.> including VCC,<cit.>  although it is not flexible enough to target strongly anharmonic systems. 
In this formalism the PES V(Q_1,…,Q_L) is expressed as a power series to encode it in a second quantization format based on Eq. (<ref>).
In addition, the operators of Eq. (<ref>) imply that the reference basis set for every mode l are the harmonic oscillator eigenfunctions.
However, such a basis does not lead to a compact representation of vibrational wave functions for strongly anharmonic systems for which modals obtained, for instance, from VSCF <cit.> are better suited.

A more flexible second quantization form is the so-called n-mode representation introduced by Christiansen <cit.>.
Instead of labelling each basis function with a single integer, as in Eq. (<ref>), we expand each mode l into a basis of N_l modals (labelled as i_1 ⋯ i_N_l) which generates an ONV basis for that mode.
Let us consider the following, general VCI expansion


    |Ψ⟩ = ∑_k_1=1^N_1⋯∑_k_L=1^N_L C_k_1,…,k_Lϕ_k_1^(1)(Q_1) ⋯ϕ_k_L^(L)(Q_L)   ,


where each mode l is described by the N_l-dimensional basis set S_l defined as


    S_l = {ϕ_1^(l) (Q_l) , … , ϕ_N_l^(l) (Q_l) }  .


The many-body basis function ϕ_k_1^(1)(Q_1) ⋯ϕ_k_L^(L)(Q_L) can be encoded as an ONV as


    ϕ_k_1(Q_1) ⋯ϕ_k_L(Q_L)
                          ≡|0_1 ⋯ 1_k_1⋯ 0_N_1,
                                       0_1 ⋯ 1_k_2⋯ 0_N_2, 
                                       ⋯ , 
                                       0_1 ⋯ 1_k_L⋯ 0_N_L⟩  .

The full ONV is then given by the expression in Eq. (<ref>) where different ONV subspaces are separated by a comma.
For each ONV space, we sort the modals in decreasing order of energy.
Each mode is described by one and only one basis function, therefore the occupation of each ONV subspace is one.

Based on the representation given in Eq. (<ref>), we introduce a pair of creation and annihilation operators per mode l and per basis function k_l defined as:


    a_k_l^†|⋯, 0_1 ⋯ 0_k_l⋯ 0_N_l, ⋯⟩   =  |⋯, 0_1 ⋯ 1_k_l⋯ 0_N_l, ⋯⟩
    
        a_k_l^†|⋯, 0_1 ⋯ 1_k_l⋯ 0_N_l, ⋯⟩   =  0 
    
        a_k_l|⋯, 0_1 ⋯ 1_k_l⋯ 0_N_l, ⋯⟩   = |⋯, 0_1 ⋯ 0_k_l⋯ 0_N_l, ⋯⟩
    
        a_k_l|⋯, 0_1 ⋯ 0_k_l⋯ 0_N_l, ⋯⟩   =  0


with


    [ a_k_l^†, a_h_m^†]    = 0 
    [ a_k_l, a_h_m]    = 0 
    [ a_k_l^†, a_h_m]    = δ_l,m  , δ_k_l,h_m


This formalism was introduced for VCC <cit.> and later applied to the multi-configurational time-dependent Hartree method <cit.>.
The second quantization form of Eq. (<ref>) obtained by expressing the potential as in Eq. (<ref>) reads<cit.>


    ℋ_vib^SQ =   ∑_l=1^L 
        ∑_k_l,h_l^N_l⟨ϕ_k_l | T(Q_l) + V^[l](Q_l) | ϕ_h_l⟩ a_k_l^+ a_h_l
    
     +   ∑_l<m^L ∑_k_l,h_l^N_l∑_k_m,h_m^N_m⟨ϕ_k_lϕ_k_m | V^[l,m](Q_l, Q_m) | ϕ_h_lϕ_h_m⟩ 
        a_k_l^+ a_k_m^+ a_h_l a_h_m + ⋯


Unlike its electronic-structure counterpart, Eq. (<ref>) contains in general coupling terms higher than two-body. Therefore, the number of Pauli terms to be evaluated on the quantum computer scales as 𝒪(N^2n) for a n-body truncation, where N is the overall number of modals.
We highlight that any PES can be encoded in the n-mode second quantization format provided that the integrals of the PES over the modals are available.
Eq. (<ref>) is therefore not restricted to PESs expressed as a power series.



 §.§ Wave function parametrization


The second quantization formalism introduced in the previous section allows one to express the VCI expansion in terms of ONVs constructed with modals that do not rely on the harmonic approximation.
The encoding of such ONVs on a quantum computer is straightforward if based on a one-to-one correspondence between modals and qubits. 
This mapping extends the direct mapping of Ref.  beyond harmonic reference basis sets.
The N_l modals for a given mode l are represented by a N_l-qubit register. 
We sort the modals in decreasing order of energy.
Therefore, the lowest-energy configuration is represented by the ONV |0_1 ⋯ 1_N_1, 0_1 ⋯ 1_N_2, ⋯, 0_1 ⋯ 1_N_L⟩ and is obtained by applying an X gate on the first qubit of each mode register initialized in the vacuum state.
The correlated wave function is obtained from the reference state by applying a set of excitation operators defined by a given wave function Ansatz.

In VQE-based electronic-structure quantum-computing two main strategies are available to prepare the wave function.
The first is based on the CC method and, more precisely, on its unitary formulation (UCC). 
It provides an intuitive expansion of the wave function in terms of excitation operators controlled by an efficiently parametrized circuit. <cit.>
However, this circuit comprises a large number of 2-qubit gates (CNOT gates) and, hence, its practical use is limited by coherence time and the gate error rates.
The second approach does not have a classical equivalent and is tailored to quantum hardware.
A heuristic wave function Ansatz is built concatenating parametrized single-qubit rotations and entangling blocks. <cit.>
The number of parameters can be increased by repeating the same set of operations (but with independent parameters) d times (where d refers to the circuit depth) to reach the desired accuracy for the ground state energy.

The same strategies can be followed for preparing vibrational wave functions. 
The UVCC circuit can be obtained from the unitary version of the VCC <cit.> Ansatz:


    |Ψ⟩ = e^𝒯-𝒯^†|Ψ_ref⟩  ,


where |Ψ_ref⟩ is the reference ONV. 𝒯 is the cluster operator (and 𝒯^† its adjoint) expressed here up to second order as


    𝒯 = 𝒯_1 + 𝒯_2    ,


with


    𝒯_1    = ∑_l^L ∑_h_l,k_l^N_lθ_h_l,k_l a^†_h_l a_k_l
    𝒯_2    = ∑_l<m^L ∑_h_l,k_l^N_l∑_h_m,k_m^N_mθ_h_l k_l h_m k_m a^†_h_l a^†_h_m a_k_m a_k_l


and (h_l,k_l) and (h_m,k_m) label couple of modals for the modes l and m, respectively.
The bosonic a^†_k_l and a_k_l operators (see Eq. (<ref>)) are mapped to the Pauli operators σ^+_k_l=σ^x_k_l + i σ^y_k_l and σ^-_k_l=σ^x_k_l - i σ^y_k_l, respectively.
In this way, the exponential operator of Eq. (<ref>) can be factorized with a Trotter expansion and expressed as a product of quantum gates.
Compared to the fermion-to-qubit mappings used in electronic-structure calculations (where the antisymmetry of the wave function is encoded in the circuit by applying for instance the Jordan Wigner transformation <cit.>) the circuit depth for the implementation of the UVCC Ansatz is greatly reduced.

Among the heuristic circuits designed for electronic structure calculation <cit.>, we consider here the SwapRZ Ansatz defined as


    |Ψ'⟩ = e^i𝒯'|Ψ_ref⟩  ,


with

    𝒯' = ∑_i<j^N_qθ_i,j (X_iX_j + Y_iY_j)   ,


where we use the notation X_i and Y_i for the σ_x and σ_y Pauli matrices acting on qubit i and N_q is the overall number of qubits.
In addition, two layers of single-qubit RZ rotations parametrized by an extra set of 2N_q angles are applied before and after the entangler block required which encodes the expansion in Eq. (<ref>).
The SwapRZ circuit ensures that the expansion for |Ψ'⟩ in made of ONVs with L and only L occupied modals. However, since the circuit also entangles pair of qubits describing modals belonging to different modes, this procedure does not ensure that a single modal per mode will be occupied.
This is not the case for the UVCCS Ansatz (UVCC with 𝒯 = 𝒯_1), where single excitations are confined to the modal space of the same mode (see Eq. (<ref>)).
Simultaneous excitations of two different modes (as those included in 𝒯_2) are not explicitly captured by the SwapRZ Ansatz with depth 1 and, therefore, we expect that deeper circuits are required to accurately represent the wave function.

Another strategy proposed in the context of electronic structure <cit.> is to build the circuit from a layer of RY and RZ rotations on each qubit followed by a block of CNOT gates entangling all qubits. 
We refer to the resulting circuit as RYRZ.
Note that rotations around the Y axis of each qubit induce a change in both the overall modals occupation and the individual occupation of each mode.

To constrain the optimization to the correct symmetry subspace, both SwapRZ and RYRZ heuristic circuits must be combined to a modified Hamiltonian ℋ_vib' where a penalty function is added to increase the energy of the states with unphysical occupation<cit.>,

    ℋ_vib' = ℋ_vib 
        + μ∑_l=1^L ( ⟨Ψ|𝒩_l |Ψ⟩ - 1)^2   ,


where μ is an arbitrary parameter and the number operator 𝒩_l for mode l is defined as

    𝒩_l = ∑_k_l=1^N_l a_k_l^† a_k_l  .

In the next section we will show that the optimization of the wave function is much more efficient with the UVCC Ansatz than with the heuristic ones.
Therefore, it is desirable to derive a quantum circuit inspired by UVCC that involves a smaller number of CNOT gates, and hence is more suited for near-term quantum calculations.

Given the current coherence time and gate error rates, it is challenging to include double excitations within the UVCC circuit.
In fact, each element of 𝒯_2 can be decomposed as:

    σ^+_iσ^+_jσ^-_kσ^-_l -  c.c. 
        =  2i(    X_iY_jX_kX_l + Y_iX_jX_kX_l + Y_iY_jX_kY_l + Y_iY_jY_kX_l 
    
            -     X_iX_jX_kY_l - X_iX_jY_kX_l - X_iY_jY_kY_l - Y_iX_jY_kY_l)   .

Therefore, the corresponding quantum circuit obtained after exponentiation of the operator given in Eq. (<ref>) and its Trotterization contains 8× 6 CNOT gates per excitation.



Combining the considerations above, we propose to approximate the UVCCSD circuit with a more compact heuristic Ansatz that we name Compact Heuristic for Chemistry (CHC). This circuit exploits the fact that the relations


    |⟨a_i^+ a_m Ψ_ref| e^iθ_m^i   X_iX_m|Ψ_ref⟩|^2 
        = | ⟨a_i^+ a_m Ψ_ref| e^θ_m^i (σ^+_iσ^-_m -
            c.c.)|Ψ_ref⟩ |^2


and


    | ⟨a_i^+ a_j^+ a_m a_n Ψ_ref| e^iθ_m,n^i,j  X_iX_jX_mX_n|Ψ_ref⟩|^2
        = | ⟨a_i^+ a_j^+ a_m a_n Ψ_ref|
            e^θ_m,n^i,j (σ^+_iσ^+_jσ^-_mσ^-_n 
            - c.c.)|Ψ_ref⟩|^2


hold for indices m, n and i, j corresponding to occupied and unoccupied modals in the reference state, respectively.
However the relative phase of the configurations in the resulting state differs. 
To correct this phase difference, we introduce the compact circuits U_s^m,i(θ_m^i) and U_d^m,n,i,j(θ_m,n^i,j), presented in Figures <ref>a and <ref>b respectively, for which the above relations become


    U^m,i_s(θ_m^i) |Ψ_ref⟩ 
        = e^θ_m^i (σ^+_iσ^-_m - c.c.)|Ψ_ref⟩  ,


and


    U^m,n,i,j_d(θ_m,n^i,j) |Ψ_ref⟩ 
        = e^θ_m,n^i,j 
          (σ^+_iσ^+_jσ^-_mσ^-_n - c.c.)|Ψ_ref⟩


for variable parameters θ_m^i and θ_m,n^i,j.
The shortcomings of CHC are related to the fact that the U^m,i_s(θ_m^i) and U^m,n,i,j_d(θ_m,n^i,j) operators are applied sequentially for each excitation. Therefore, in general, the circuit that corresponds to a given excitation (shown in Fig <ref>) is not directly applied on the reference state |Ψ_ref⟩ but rather on a superposition state generated through the previous excitations.
This can produce unphysical configurations with the wrong number of particles for each mode.
For weakly correlated systems, we expect CHC to act very similarly to UVCC as the weight of the unphysical configurations will be negligible.
However, this approximation deteriorates when working with strongly correlated systems, i.e. with a VCI wave function that is not dominated by a single configuration.
The decisive advantage of CHC is that the number of CNOT gates is reduced by approximately one order of magnitude compared to UCC, allowing therefore the computation of larger systems.
We note that the unphysical configurations can also be efficiently projected out with a Monte-Carlo inspired scheme described in Ref. . 
The scaling and performance of the CHC circuit are presented in Section <ref>.



 §.§ Extension to excited states


The calculation of the ground state is not sufficient for most vibrational-structure calculations for which vibrational excitation energies must also be considered.
This is the case for the simulation of vibrational spectra, in which peaks are located at transition frequencies.
In quantum computing, the calculation of excited states can be performed with the quantum EOM (qEOM) algorithm<cit.> in which a vibrational excited state |n⟩ is expressed as |n⟩ = O_n^†|0⟩, where |0⟩ indicates the vibrational ground state and O_n^† is defined as


    O_n^† 
          = ∑_α∑_μ_α[ X_μ_α^(α) (n) ℰ_μ_α^(α) 
              - Y_μ_α^(α)(n) (ℰ_μ_α^(α))^†]   .


In Eq. (<ref>), α is the order of the excitation and the collective index μ_α runs over all modals involved in the excitation. The excitation operators ℰ_μ_α^(α) are sequences of creation and annihilation operators. For instance ℰ_μ_2^(2) = a_k_l^† a_h_m^† a_0_m a_0_l, with l<m running over the number of modes and k_l,h_m over the number of modals for modes l and m, respectively. In addition, X_μ_α^(α) is the expansion coefficient for the ℰ_μ_α^(α) excitation operator, and Y_μ_α^(α) is the coefficient for the corresponding de-excitation.
The excitation energy E_0n is found by solving classically the following pseudo-eigenvalue problem


    [  M  Q; Q* M* ][ X_n; Y_n ]
      = E_0n[   V   W; -W* -V* ][ X_n; Y_n ] ,


where the matrix elements


    M_μ_αν_β = ⟨0| [(ℰ_μ_α^(α))^†,ℋ, ℰ_ν_β^(β)]|0⟩, 
        Q_μ_αν_β = -⟨0| [(ℰ_μ_α^(α))^†, ℋ, (ℰ_ν_β^(β))^†]|0⟩,
        V_μ_αν_β = ⟨0| [(ℰ_μ_α^(α))^†, ℰ_ν_β^(β)]|0⟩,
       W_μ_αν_β = -⟨0| [(ℰ_μ_α^(α))^†, (ℰ_ν_β^(β))^†]|0⟩.


are evaluated with the approximated ground state wavefunction generated with the quantum circuit.
The accuracy of qEOM can be improved systematically by increasing the maximum excitation order α.
The inclusion of excitations with α > 2 may become important for PESs that require the inclusion of three- and higher-order terms in the n-body expansion of Eq. (<ref>).
This increases the size of the pseudo-eigenvalue problem and, therefore, the number of measurements to be performed on the quantum hardware.



§ RESULTS


As an example system, we choose the carbon dioxide molecule, which is well studied with traditional approaches<cit.>.
We also estimate the quantum computing resources needed for the simulation of two larger molecules, namely formaldehyde (H_2CO) and formic acid (HCOOC).



 §.§ Ground state calculations with state-of-the-art approaches



We study the UVCC, SwapRZ and RYRZ wave function approaches on the PES of CO_2 defined by the bending and the symmetric stretching modes.
We describe the system Hamiltonian in second quantization as in Eq. (<ref>) with two modals per mode where the reference modal basis (see Eq. (<ref>)) is obtained as eigenfunctions of the one-body Hamiltonian


    ( T(Q_l) + V^[l](Q_l) ) ϕ_k_l(Q_l)
        = ϵ_k_lϕ_k_l(Q_l)      (l=1,2)   .


One could include vibrational correlations in Eq. (<ref>) by introducing the averaged two- and higher-order couplings with the other modes in the potential operator, as done in VSCF <cit.>.
The resulting set of modals would lead to more compact VCI and VCC expansions than those based on the modals from Eq. (<ref>).
However, modals obtained from Eq. (<ref>) are sufficient to assess the quality of our algorithm and demonstrate that bases different from the harmonic ones can be used.
We optimize the ground-state equilibrium geometry of CO_2 using density functional theory with the B3LYP exchange-correlation functional <cit.> and the cc-pVTZ <cit.> basis set.
We approximate the PES with a quartic force field and calculate the anharmonic force fields by semi-numerical differentiation of the analytical Hessian as implemented in  <cit.>.

All VQE calculations are run with  <cit.>. For the SwapRZ and RYRZ circuits, we set the penalty term μ in Eq. (<ref>) to 10^5 (we find empirically that values below this threshold do not proscribe the convergence to the vacuum state).
For all VQE simulations, we apply the exact unitary matrix representation of the circuit on the reference ONV without taking into account sampling, decoherence and gate noise.
The results of the simulations are reported in Figure <ref>a.



We note that the convergence for both the SwapRZ and the RYRZ circuits is generally slower than with UVCCSD and convergence is reached only after 8000 iterations for SwapRZ, and after more than 10000 iterations for RYRZ. 
We recall that UVCC is the only occupation-conserving Ansatz, therefore this pronounced difference can be due to the modified Hamiltonian of Eq. (<ref>).
The RYRZ circuit cannot be used without the penalty term since it does not conserve the overall modal occupation, and therefore it would converge to the vacuum state.
Conversely, even if the SwapRZ does not ensure that only one modal per mode is occupied, the energy gap between two modals is high enough to prevent the algorithm to converge to states with the incorrect occupation such as, for instance, the state |11,00⟩.
Hence, we combine the SwapRZ Ansatz with the original non-modified Hamiltonian to assess the effect of the addition of the penalty term, and compare the results with the UVCC approach. 
We run a VQE with the SwapRZ circuit for depths 1, 2 and 3. The number of entangling gates is 56, 24, 48, and 72 for UVCC, SwapRZ1, SwapRZ2, and SwapRZ3, which corresponds to a number of variational parameters of 3, 14, 24, and 34, respectively.
Fig. <ref>b shows that, even if the optimization is faster with the unmodified Hamiltonian, UVCC still outperforms SwapRZ for all depths.
This suggests that the explicit inclusion of double-excitations in the wave function Ansatz is crucial in order to obtain reliable energies, even if this also leads to deeper circuits due to the presence of multiple CNOT gates.



 §.§ Performance of the CHC Ansatz


To reduce the circuit depth without sacrificing the accuracy, we approximate the UVCC wave function with the CHC Ansatz (Eqs. (<ref>) and (<ref>)).
The scaling of UVCC and CHC in terms of number of CNOT gates and number of parameters is shown on Table <ref> for three molecules: CO_2, H_2CO and HCOOC.
For all these cases, we constructed the reference modal basis as described for CO_2, i.e. from the eigenfunctions of the one-body Hamiltonian of Eq. (<ref>), based on geometries optimized with B3LYP/cc-pVTZ and calculating the anharmonic quartic force field by including all third-order and the semi-diagonal fourth-order terms.



We study the CHC circuit for CO_2 with different number of modes and modals per mode (see Figure <ref>).
The results are of reasonable accuracy, with deviations <15 cm^-1 compared to the exact diagonalization energies.
As expected, CHC works best with two modals per mode since the entangling blocks are always applied to a state of the subspace with the correct symmetry, and therefore it does not create configurations with the wrong number of particles.
Errors increase with the number of modals for each mode.





 §.§ Quantum computation of vibrational structures on existing hardware


In presence of typical hardware noise the loss in accuracy inherent to CHC will be balanced by the reduced circuit depth.
To study this effect, we use both UVCC and CHC in simulations with a noise model and in hardware calculations with the ibmq_almaden 20-qubit processor.
For this study, we fix the number of modes to two and progressively increase the number of modals per mode.
The gate angles, {θ_i}, defining the circuit gates are initiated randomly and restricted to angles -0.2≤θ_i≤0.2 [For small angles the reference configuration remains largely dominant. Hence, CHC is expected to lead to a reasonable approximation of UVCC for the same set of variational parameters.].
The noise model includes only depolarization errors for single- and two-qubit gates, for which the error rate is based on the the average gate depolarization error associated to all qubits in the ibmq_almaden 20-qubit device.
The error rate values are 7 × 10^-4, 1.4 × 10^-3 and 2.2 × 10^-2 for U2, U3 and CNOT gates, respectively.
For each set of randomly parametrized UVCC and CHC circuits we evaluate the underlying probability distributions by performing 10000 measurements on the final states.
The final states are then compared to the reference obtained with the exact (noise-free) simulation of the UVCC Ansatz (note that in both cases UVCC is the reference since we aim at approximating UVCC with CHC and assess the loss in accuracy due to the approximation versus the gain due to less noise).
This process is repeated 10 times renewing each time the parameter set. 
The fidelity of the resulting probability distribution is given by:


    ℱ^circ = 1 - ∑_i^2^N_q|C^circ_i-C^ref_i|/∑_i^2^N_qC^circ_i+C^ref_i


where C_i is the count for the occurrence of state i, with i ranging over all possible 2^N_q states (N_q being the number of qubits).
The fidelities obtained with both UVCC and CHC are shown as a function of system size in Fig. <ref>a.
The number of CNOTs is also given in Fig. <ref>a for two modes and two modals per mode, and for four modes and four modals per mode.
These results prove that in the presence of noise a better accuracy is reached with CHC rather than with UVCC.



We repeat the experiment with a quantum processors (ibmq_almaden, 20 qubits) and a system size corresponding to two modes and two modals per mode (four qubits). 
Figure <ref>b shows the histogram of the states distribution corresponding to a single trial (one set of gate parameters) confirming our observation that, with the present hardware, compact heuristic circuits such as the ones obtained with CHC outperform the original UVCC Ansatz. 



 §.§ Quantum computation of the vibrational excited states of CO_2


We calculate the excitation energies of a CO_2 molecule with (A) two modes, two modals per mode; (B) two modes, four modals per mode and (C) four modes, two modals per mode.
For all three cases, we limit the qEOM operators to 1- and 2-body excitation operators to restrict the number of measurements required to evaluate the matrix elements in Eq. (<ref>).
The ground state is approximated by running a VQE calculation with both UVCC and CHC.
The results of these (noise-free) simulations are presented in Table <ref>.
The reference values are obtained from the exact diagonalization of the system Hamiltonian.

In case (A) all excitation energies are found with an accuracy < 1 cm^-1.
For (B), the accuracy is lower for CHC compared to UVCC. 
This is expected since for a large number of modals the accuracy of CHC ground state calculations decreases.
Finally, case (C) shows that higher order excitations need to be included in ℰ_μ_α^(α) in order to reach an accuracy of about 1 cm^-1 for the highest excitation energies with both wave function Ansätze.





§ DISCUSSION


The proposed quantum algorithms for the calculation of the vibrational frequencies support PES representations and modals routinely employed in state-of-the-art traditional calculations. 
This enables a fair comparison of the scaling of the classical and quantum vibrational structure algorithms that provides an estimate of the resources required to reach a quantum advantage in vibrational-structure calculations.
As mentioned in the introduction, large-scale vibrational structure calculations are nowadays possible either with efficient VCI algorithms <cit.> or with non-linear wave function parametrizations, as is done in vDMRG <cit.> and VCC <cit.>.
The latter algorithm is the direct classical counterpart of our UVCC-based quantum algorithm and will be our reference for comparing the scaling with its quantum counterpart.

The scaling of VCC depends on the order of the highest degree of excitation that is included in the 𝒯 operator (Eq. (<ref>)).
For instance, VCC[2pt3] comprises all two-mode excitations and treats triple excitations perturbatively and can be applied, in its straightforward formulation, to molecules with up to seven atoms, such as ethylene oxide <cit.> including six modals for each vibrational mode.
The simulation of such molecules on quantum hardware would require 90 qubits.
The corresponding circuit that includes single and double excitations, and approximates triple excitations with the CHC Ansatz (which we denote as UVCCSD(T)) would contain about 10^6 CNOT gates.
By approximating also the single and double excitations with CHC, the number of 2-qubit gates drops to about 10^4.
Currently, the state-of-the-art quantum hardware comprises about 50 qubits and has a coherence time supporting circuits with no more than 10^2 CNOT gates.
The 2-qubit error rate of about 10^-2 is currently the limiting factor for running such circuits in the state-of-the-art hardware.
By improving the 2-qubit gate fidelity and according to the estimated evolution of the quantum volume in superconducting
quantum computers <cit.>, molecules of the dimensions of ethylene oxide will become accessible using the proposed algorithm within the next generation of quantum hardware. 

Molecules with up to seven atoms and described by a maximum of four-mode excitations can be studied with VCC by adopting a tensor-factorized representation of the amplitudes.<cit.>
The inclusion of three- and four-body coupling terms in the potential and in the VCC wave function leaves the number of qubits required unchanged.
However, it also induces a raise in both the circuit depth and the number of measurements i.e., number of terms in the Hamiltonian.
For the electronic Hamiltonian, Motta and co-workers showed <cit.> that tensor factorizations can be used to reduce the number of measurements and circuit depth for fermionic systems.
We expect that the same holds true also for the vibrational case, and we will consider such extension in future works.

The computational cost of VCC depends also on the choice of the coordinates used to describe the Hamiltonian.
It is known that for particular choices, such as local modes<cit.> or VSCF-optimized coordinates, the size of the off-diagonal anharmonic couplings can be significantly reduced.
This procedure has been already exploited to speed up traditional VCC calculations <cit.> including only excitations between modes localized on the same portion of the molecules. 
This simplification has made VCC calculation feasible for systems as large as the water hexamer <cit.>.
Our algorithm supports any choice of the reference coordinate system.
In the same way, the UVCC circuit can be adapted to a local mode representation by allowing only gates representing excitations for modes localized on nearby portions of the molecule.
For instance, in water clusters <cit.> one could apply the UVCC circuit for excitations localized on one water molecule and include an approximated treatment of inter-fragment correlations with CHC.



§ CONCLUSIONS


We designed and compared different quantum computing strategies to calculate the vibrational structure of molecular systems amenable to near-term quantum computers.
We represented the vibrational wave function and PES in the n-mode-based Fock space <cit.> that supports an arbitrary one-body reference basis and PES expressed as a generic many-body expansion.
This enabled us to overcome the limitations of recent algorithms <cit.> that rely on a harmonic-based reference and are, therefore, not flexible enough for strongly anharmonic systems.
We compared state-of-the-art circuits to prepare the wave function of a two-dimensional Hamiltonian modelling the nuclear dynamics of CO_2 and introduced the Compact Heuristic circuit for Chemistry (CHC).
On the one hand, the Unitary Vibrational Coupled Cluster (UVCC) delivers the most accurate vibrational energies, but at the price of very deep circuits that are difficult to implement on currently available quantum hardware already for three-atom molecules.
On the other hand, heuristic circuits provide less accurate vibrational energies, while being shallower.
In this work, we showed how CHC represents an optimal compromise combining the advantages of UVCC and heuristic wave function approaches.
However, the CHC wave function does not fulfill the symmetries of the vibrational Hamiltonian.
Therefore one first needs to project it onto the correct symmetry subspace before evaluating the vibrational energy.
This effect is only minor for the systems studied here, but we expect it becomes larger for strongly anharmonic molecules.
In those cases, a modal basis obtained from a VSCF calculation can significantly improve the accuracy.
A second limitation of our algorithm is that each modal is mapped to a different qubit.
Therefore, a large portion of the qubit Hilbert space does not correspond to a physically acceptable state.
Occupation number vectors based on alternative mappings, such as the ones introduced in Ref. , produce more compact representations of the vibrational states.
However, the reduction in the number of qubits comes at the cost of an increase in the circuit depth for the representation of the wave function.
Heuristic circuits inspired by the CHC strategy, but adapted to these more compact mappings, could enable the calculation of vibrational energies for molecules with more than three atoms on a state-of-the-art quantum computer.
Finally, we extended the quantum Equation of Motion (qEOM)<cit.> approach to calculate vibrational excitation energies and applied it to the CO_2 molecule.

This work also sets the fundamentals for the quantum computation of the ground-state energy of interacting fermions and bosons, such as polaronic <cit.> or quantum optics Hamiltonians <cit.>.
In the quantum-chemistry context, this is the case of the pre-Born-Oppenheimer molecular Hamiltonian <cit.> that has been studied so far with the quantum phase estimation algorithm.<cit.>
Moreover, the algorithm can be extended to the time domain (e.g., using the time-dependent Schrödinger formalism of Ref. ) to address quantum dynamics.


This work was supported by ETH Zürich through the ETH Fellowship No. FEL-49 18-1.

The authors acknowledge financial support from the Swiss National Science Foundation (SNF) through the grant No. 200021-179312.

IBM, IBM Quantum,  are trademarks of International Business Machines Corporation, registered in many jurisdictions worldwide. Other product or service names may be trademarks or service marks of IBM or other companies.






@ifundefinedendmcitethebibliography
  
72








f
subitem(mcitesubitemcount)

  
  

[Bowman(1978)]bowman1978
Bowman, J. M. Self-consistent field energies and wavefunctions for coupled
  oscillators. J. Chem. Phys. 1978, 68, 608–610
[Carney et al.(1978)Carney, Sprandel, and Kern]carney1978
Carney, G. D.; Sprandel, L. L.; Kern, C. W. Variational Approaches to
  Vibration-Rotation Spectroscopy for Polyatomic Molecules. Adv. Chem.
  Phys. 1978, 37, 305–379
[Gerber and Ratner(1979)Gerber, and Ratner]gerber1979
Gerber, R.; Ratner, M. A. A semiclassical self-consistent field (SC SCF)
  approximation for eigenvalues of coupled-vibration systems. Chem. Phys.
  Lett. 1979, 68, 195–198
[Roy and Gerber(2013)Roy, and Gerber]Gerber2013_Review
Roy, T. K.; Gerber, R. B. Vibrational self-consistent field calculations for
  spectroscopy of biological molecules: new algorithmic developments and
  applications. Phys. Chem. Chem. Phys. 2013, 15,
  9468–9492
[Christiansen(2003)]christiansen2003
Christiansen, O. Møller–Plesset perturbation theory for vibrational wave
  functions. J. Chem. Phys. 2003, 119, 5773–5781
[Barone(2005)]Barone2005_VPT2
Barone, V. Anharmonic vibrational properties by a fully automated second-order
  perturbative approach. J. Chem. Phys. 2005, 122,
  014108
[Bowman et al.(1979)Bowman, Christoffel, and Tobin]bowman1979
Bowman, J. M.; Christoffel, K.; Tobin, F. Application of SCF-CI theory to
  vibrational motion in polyatomic molecules. J. Phys. Chem.
  1979, 83, 905–912
[Thompson and Truhlar(1980)Thompson, and Truhlar]thompson1980
Thompson, T. C.; Truhlar, D. G. SCF CI calculations for vibrational eigenvalues
  and wavefunctions of systems exhibiting Fermi resonance. Chem. Phys.
  Lett. 1980, 75, 87–90
[Christoffel and Bowman(1982)Christoffel, and Bowman]christoffel1982
Christoffel, K. M.; Bowman, J. M. Investigations of self-consistent field, scf
  ci and virtual stateconfiguration interaction vibrational energies for a
  model three-mode system. Chem. Phys. Lett. 1982, 85,
  220–224
[Neff and Rauhut(2009)Neff, and Rauhut]Rauhut2009_VCI-Fast
Neff, M.; Rauhut, G. Toward large scale vibrational configuration interaction
  calculations. J. Chem. Phys. 2009, 131, 124129
[Christiansen(2004)]christiansen2004
Christiansen, O. Vibrational coupled cluster theory. J. Chem. Phys.
  2004, 120, 2149–2159
[Seidler and Christiansen(2009)Seidler, and
  Christiansen]Christiansen2009_Automatic-VCC
Seidler, P.; Christiansen, O. Automatic derivation and evaluation of
  vibrational coupled cluster theory equations. J. Chem. Phys.
  2009, 131, 234109
[Oschetzki and Rauhut(2014)Oschetzki, and
  Rauhut]Rauhut2014_LiFClusters
Oschetzki, D.; Rauhut, G. Pushing the limits in accurate vibrational structure
  calculations: anharmonic frequencies of lithium fluoride clusters (LiF)n, n
  = 2–10. Phys. Chem. Chem. Phys. 2014, 16,
  16426–16435
[Thomas et al.(2018)Thomas, Carrington, Agarwal, and
  Schaefer]Carrington2018_UracilNaphtalene
Thomas, P. S.; Carrington, T.; Agarwal, J.; Schaefer, H. F. Using an iterative
  eigensolver and intertwined rank reduction to compute vibrational spectra of
  molecules with more than a dozen atoms: Uracil and naphthalene. J.
  Chem. Phys. 2018, 149, 064108
[Yu and Bowman(2019)Yu, and Bowman]Bowman2019_EigenWater
Yu, Q.; Bowman, J. M. Classical, Thermostated Ring Polymer, and Quantum
  VSCF/VCI Calculations of IR Spectra of H_7O_3^+ and H_9O_4^+ (Eigen)
  and Comparison with Experiment. J. Phys. Chem. A 2019,
  123, 1399–1409
[Baiardi et al.(2017)Baiardi, Stein, Barone, and
  Reiher]Baiardi2017_vDMRG
Baiardi, A.; Stein, C. J.; Barone, V.; Reiher, M. Vibrational Density Matrix
  Renormalization Group. J. Chem. Theory Comput. 2017,
  13, 3764–3777
[Baiardi et al.(2019)Baiardi, Stein, Barone, and
  Reiher]Baiardi2019_ExcitedStates-vDMRG
Baiardi, A.; Stein, C. J.; Barone, V.; Reiher, M. Optimization of highly
  excited matrix product states with an application to vibrational
  spectroscopy. J. Chem. Phys. 2019, 150, 094113
[Baiardi and Reiher(2020)Baiardi, and Reiher]Reiher2020_Review
Baiardi, A.; Reiher, M. The density matrix renormalization group in chemistry
  and molecular physics: Recent developments and new challenges. J. Chem.
  Phys. 2020, 152, 040903
[Peruzzo et al.(2014)Peruzzo, McClean, Shadbolt, Yung, Zhou,
  Love, Aspuru-Guzik, and O'Brien]peruzzo2014
Peruzzo, A.; McClean, J.; Shadbolt, P.; Yung, M.-H.; Zhou, X.-Q.; Love, P. J.;
  Aspuru-Guzik, A.; O'Brien, J. L. A variational eigenvalue solver on a
  photonic quantum processor. Nat. Comm. 2014, 5
[Yung et al.(2014)Yung, Casanova, Mezzacapo, Mcclean, Lamata,
  Aspuru-Guzik, and Solano]yung2014
Yung, M.-H.; Casanova, J.; Mezzacapo, A.; Mcclean, J.; Lamata, L.;
  Aspuru-Guzik, A.; Solano, E. From transistor to trapped-ion computers for
  quantum chemistry. Sci. Rep. 2014, 4, 3589
[McClean et al.(2016)McClean, Romero, Babbush, and
  Aspuru-Guzik]mcclean2016
McClean, J. R.; Romero, J.; Babbush, R.; Aspuru-Guzik, A. The theory of
  variational hybrid quantum-classical algorithms. New J. Phys.
  2016, 18, 023023
[Wang et al.(2019)Wang, Higgott, and Brierley]wang2019
Wang, D.; Higgott, O.; Brierley, S. Accelerated variational quantum
  eigensolver. Phys. Rev. Lett. 2019, 122, 140504
[Kandala et al.(2017)Kandala, Mezzacapo, Temme, Takita, Brink,
  Chow, and Gambetta]kandala2017
Kandala, A.; Mezzacapo, A.; Temme, K.; Takita, M.; Brink, M.; Chow, J. M.;
  Gambetta, J. M. Hardware-efficient variational quantum eigensolver for small
  molecules and quantum magnets. Nature 2017, 549,
  242
[Nam et al.(2019)Nam, Chen, Pisenti, Wright, Delaney, Maslov,
  Brown, Allen, Amini, Apisdorf, Beck, Blinov, Chaplin, Chmielewski, Collins,
  Debnath, Ducore, Hudek, Keesan, Kreikemeier, Mizrahi, Solomon, Williams,
  David Wong-Campos, Monroe, and Kim]nam2019
Nam, Y.; Chen, J.-S.; Pisenti, N. C.; Wright, K.; Delaney, C.; Maslov, D.;
  Brown, K. R.; Allen, S.; Amini, J. M.; Apisdorf, J.; Beck, K. M.; Blinov, A.;
  Chaplin, V.; Chmielewski, M.; Collins, C.; Debnath, S.; Ducore, A. M.;
  Hudek, K. M.; Keesan, M.; Kreikemeier, S. M.; Mizrahi, J.; Solomon, P.;
  Williams, M.; David Wong-Campos, J.; Monroe, C.; Kim, J. Ground-state energy
  estimation of the water molecule on a trapped ion quantum computer.
  arXiv preprint arXiv:1902.10171 2019, 
[Ollitrault et al.(2019)Ollitrault, Kandala, Chen, Barkoutsos,
  Mezzacapo, Pistoia, Sheldon, Woerner, Gambetta, and
  Tavernelli]ollitrault2019
Ollitrault, P. J.; Kandala, A.; Chen, C.-F.; Barkoutsos, P. K.; Mezzacapo, A.;
  Pistoia, M.; Sheldon, S.; Woerner, S.; Gambetta, J.; Tavernelli, I. Quantum
  equation of motion for computing molecular excitation energies on a noisy
  quantum processor. arXiv preprint arXiv:1910.12890 2019,
  
[Colless et al.(2018)Colless, Ramasesh, Dahlen, Blok,
  Kimchi-Schwartz, McClean, Carter, De Jong, and Siddiqi]colless2018
Colless, J. I.; Ramasesh, V. V.; Dahlen, D.; Blok, M. S.; Kimchi-Schwartz, M.;
  McClean, J.; Carter, J.; De Jong, W.; Siddiqi, I. Computation of molecular
  spectra on a quantum processor with an error-resilient algorithm. Phys.
  Rev. X 2018, 8, 011021
[Chaban et al.(1999)Chaban, Jung, and Gerber]Gerber1999_nMode
Chaban, G. M.; Jung, J. O.; Gerber, R. B. Ab initio calculation of anharmonic
  vibrational states of polyatomic systems: Electronic structure combined with
  vibrational self-consistent field. J. Chem. Phys. 1999,
  111, 1823–1829
[Kongsted and Christiansen(2006)Kongsted, and
  Christiansen]Kongsted2006_nModeRepresentation
Kongsted, J.; Christiansen, O. Automatic generation of force fields and
  property surfaces for use in variational vibrational calculations of
  anharmonic vibrational energies and zero-point vibrational averaged
  properties. J. Chem. Phys. 2006, 125, 124108
[McArdle et al.(2019)McArdle, Mayorov, Shan, Benjamin, and
  Yuan]mcardle2019
McArdle, S.; Mayorov, A.; Shan, X.; Benjamin, S.; Yuan, X. Digital quantum
  simulation of molecular vibrations. Chem. Sci. 2019,
  10, 5725–5735
[Sawaya et al.(2019)Sawaya, Menke, Kyaw, Johri, Aspuru-Guzik,
  and Guerreschi]sawaya2019
Sawaya, N. P.; Menke, T.; Kyaw, T. H.; Johri, S.; Aspuru-Guzik, A.;
  Guerreschi, G. G. Resource-efficient digital quantum simulation of d-level
  systems for photonic, vibrational, and spin-s Hamiltonians. arXiv
  preprint arXiv:1909.12847 2019, 
[Wang and Thoss(2009)Wang, and Thoss]Wang2009_SQMCTDH
Wang, H.; Thoss, M. Numerically exact quantum dynamics for indistinguishable
  particles: The multilayer multiconfiguration time-dependent Hartree theory in
  second quantization representation. J. Chem. Phys. 2009,
  131, 24114
[Taube and Bartlett(2006)Taube, and Bartlett]Bartlett2006_UCC-Review
Taube, A. G.; Bartlett, R. J. New perspectives on unitary coupled-cluster
  theory. Int. J. Quantum Chem. 2006, 106,
  3393–3401
[Bartlett and Musiał(2007)Bartlett, and Musiał]bartlett2007
Bartlett, R. J.; Musiał, M. Coupled-cluster theory in quantum chemistry.
  Rev. Mod. Phys. 2007, 79, 291
[Lee et al.(2019)Lee, Huggins, Head-Gordon, and
  Whaley]Lee2019_Generalized-UCC
Lee, J.; Huggins, W. J.; Head-Gordon, M.; Whaley, K. B. Generalized Unitary
  Coupled Cluster Wave functions for Quantum Computation. J. Chem. Theory
  Comput. 2019, 15, 311–324
[Evangelista et al.(2019)Evangelista, Chan, and
  Scuseria]Evangelista2019_UCC
Evangelista, F. A.; Chan, G. K.-L.; Scuseria, G. E. Exact parameterization of
  fermionic wave functions via unitary coupled cluster theory. J. Chem.
  Phys. 2019, 151, 244112
[Wilson et al.(1980)Wilson, Declus, and Cross]WilsonBook
Wilson, E. B.; Declus, J. C.; Cross, P. C. In Molecular Vibrations: The
  Theory of Infrared and Raman Vibrational Spectra; Dover,, Ed.; Dover
  Publications; New edition edition (March 1, 1980), 1980
[Papousek and Aliev(1982)Papousek, and Aliev]PapousekBook
Papousek, D.; Aliev, M. R. In Molecular Vibrational-rotational Spectra:
  Theory and Applications of High Resolution Infrared, Microwave and Raman
  Spectroscopy of Polyatomic Molecules; Elsevier Science,, Ed.; 1982
[Hermes and Hirata(2014)Hermes, and Hirata]Hirata2014_StochasticVSCF
Hermes, M. R.; Hirata, S. Stochastic algorithm for size-extensive vibrational
  self-consistent field methods on fully anharmonic potential energy surfaces.
  J. Chem. Phys. 2014, 141, 244111
[Hirata and Hermes(2014)Hirata, and Hermes]Hirata2014_NormalOrdered
Hirata, S.; Hermes, M. R. Normal-ordered second-quantized Hamiltonian for
  molecular vibrations. J. Chem. Phys. 2014, 141,
  184111
[Nagalakshmi et al.(1994)Nagalakshmi, Lakshminarayana,
  Sumithra, and Prasad]Prasad1994_Original
Nagalakshmi, V.; Lakshminarayana, V.; Sumithra, G.; Prasad, M. D. Coupled
  cluster description of anharmonic molecular vibrations. Application to O_3
  and SO_2. Chem. Phys. Lett. 1994, 217, 279
[Banik et al.(2008)Banik, Pal, and
  Prasad]Prasad2008_VCC-Response
Banik, S.; Pal, S.; Prasad, M. D. Calculation of vibrational energy of molecule
  using coupled cluster linear response theory in bosonic representation:
  Convergence studies. J. Chem. Phys. 2008, 129,
  134111
[Faucheaux et al.(2018)Faucheaux, Nooijen, and
  Hirata]Hirata2018_SimilarityTransformedVCC
Faucheaux, J. A.; Nooijen, M.; Hirata, S. Similarity-transformed
  equation-of-motion vibrational coupled-cluster theory. J. Chem. Phys.
  2018, 148, 054104
[Faucheaux and Hirata(2015)Faucheaux, and
  Hirata]Hirata2018_HighOrderVCC
Faucheaux, J. A.; Hirata, S. Higher-order diagrammatic vibrational
  coupled-cluster theory. J. Chem. Phys 2015, 143,
  134105
[Bowman(1986)]bowman1986
Bowman, J. M. The self-consistent-field approach to polyatomic vibrations.
  Acc. Chem. Res. 1986, 19, 202–208
[Christiansen(2004)]Christiansen2004_nMode
Christiansen, O. A second quantization formulation of multimode dynamics.
  J. Chem. Phys 2004, 120, 2140–2148
[McArdle et al.(2018)McArdle, Endo, Aspuru-Guzik, Benjamin, and
  Yuan]mcardle2018
McArdle, S.; Endo, S.; Aspuru-Guzik, A.; Benjamin, S.; Yuan, X. Quantum
  computational chemistry. arXiv preprint arXiv:1808.10402
  2018, 
[Barkoutsos et al.(2018)Barkoutsos, Gonthier, Sokolov, Moll,
  Salis, Fuhrer, Ganzhorn, Egger, Troyer, Mezzacapo, Filipp, and
  Tavernelli]barkoutsos2018
Barkoutsos, P. K.; Gonthier, J. F.; Sokolov, I.; Moll, N.; Salis, G.;
  Fuhrer, A.; Ganzhorn, M.; Egger, D. J.; Troyer, M.; Mezzacapo, A.;
  Filipp, S.; Tavernelli, I. Quantum algorithms for electronic structure
  calculations: Particle-hole Hamiltonian and optimized wave-function
  expansions. Phys. Rev. A 2018, 98, 022322
[Lee et al.(2018)Lee, Huggins, Head-Gordon, and
  Whaley]lee2018
Lee, J.; Huggins, W. J.; Head-Gordon, M.; Whaley, K. B. Generalized Unitary
  Coupled Cluster Wave functions for Quantum Computation. J. Chem. Theory
  Comput. 2018, 15, 311–324
[Jordan and Wigner(1993)Jordan, and Wigner]jordan1993
Jordan, P.; Wigner, E. P. The Collected Works of Eugene Paul Wigner;
  Springer, 1993; pp 109–129
[Mazzola et al.(2019)Mazzola, Ollitrault, Barkoutsos, and
  Tavernelli]mazzola2019
Mazzola, G.; Ollitrault, P. J.; Barkoutsos, P. K.; Tavernelli, I. Nonunitary
  Operations for Ground-State Calculations in Near-Term Quantum Computers.
  Phys. Rev. Lett. 2019, 123, 130501
[Rodriguez-Garcia et al.(2007)Rodriguez-Garcia, Hirata, Yagi,
  Hirao, Taketsugu, Schweigert, and Tasumi]Hirata2007_CO2-PES
Rodriguez-Garcia, V.; Hirata, S.; Yagi, K.; Hirao, K.; Taketsugu, T.;
  Schweigert, I.; Tasumi, M. Fermi resonance in CO2: A combined electronic
  coupled-cluster and vibrational configuration-interaction prediction.
  J. Chem. Phys. 2007, 126, 124303
[Sode et al.(2013)Sode, Keceli, Yagi, and
  Hirata]Hirata2013_CO2-Pressure
Sode, O.; Keceli, M.; Yagi, K.; Hirata, S. Fermi resonance in solid CO2 under
  pressure. J. Chem. Phys. 2013, 138, 074501
[Wang and Bowman(2017)Wang, and Bowman]Bowman2017_CO2-Copper
Wang, Q. K.; Bowman, J. M. Two-component, ab initio potential energy surface
  for CO2—H2O, extension to the hydrate clathrate, CO2@(H2O)20, and VSCF/VCI
  vibrational analyses of both. J. Chem. Phys. 2017,
  147, 161714
[Becke(1993)]becke1993
Becke, A. D. A new mixing of Hartree–Fock and local density-functional
  theories. J. Chem. Phys. 1993, 98, 1372–1377
[Dunning(1989)]Dunning1989_Basis
Dunning, T. H. Gaussian basis sets for use in correlated molecular
  calculations. I. The atoms boron through neon and hydrogen. J. Chem.
  Phys. 1989, 90, 1007–1023
[Frisch et al.(2009)Frisch, Trucks, Schlegel, Scuseria, Robb,
  Cheeseman, Scalmani, Barone, Petersson, Nakatsuji, Li, Caricato, Marenich,
  Bloino, Janesko, Gomperts, Mennucci, Hratchian, Ortiz, Izmaylov, Sonnenberg,
  Williams-Young, Ding, Lipparini, Egidi, Goings, Peng, Petrone, Henderson,
  Ranasinghe, Zakrzewski, Gao, Rega, Zheng, Liang, Hada, Ehara, Toyota, Fukuda,
  Hasegawa, Ishida, Nakajima, Honda, Kitao, Nakai, Vreven, Throssell,
  J. A. Montgomery, Peralta, Ogliaro, Bearpark, Heyd, Brothers, Kudin,
  Staroverov, Keith, Kobayashi, Normand, Raghavachari, Rendell, Burant,
  Iyengar, Tomasi, Cossi, Millam, Klene, Adamo, Cammi, Ochterski, Martin,
  Morokuma, Farkas, Foresman, and Fox]g09
Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.;
  Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.;
  Li, X.; Caricato, M.; Marenich, A.; Bloino, J.; Janesko, B. G.; Gomperts, R.;
  Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.;
  Sonnenberg, J. L.; Williams-Young, D.; Ding, F.; Lipparini, F.; Egidi, F.;
  Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.;
  Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.;
  Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.;
  Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.;
  J. A. Montgomery, J.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.;
  Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.;
  Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.;
  Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.;
  Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.;
  Fox, D. J. Gaussian˜09 Revision C.01. 2009; Gaussian Inc. Wallingford
  CT
[Abraham et al.(2019)Abraham, Akhalwaya, Aleksandrowicz,
  Alexander, Alexandrowics, Arbel, Asfaw, Azaustre, Barkoutsos, et al.
  others]Qiskit
Abraham, H.; Akhalwaya, I. Y.; Aleksandrowicz, G.; Alexander, T.;
  Alexandrowics, G.; Arbel, E.; Asfaw, A.; Azaustre, C.; Barkoutsos, P. K.,
  et al.  Qiskit: An Open-source Framework for Quantum Computing.
  2019
[Scribano et al.(2010)Scribano, Lauvergnat, and
  Benoit]Benoit2010_FastVCI
Scribano, Y.; Lauvergnat, D. M.; Benoit, D. M. Fast vibrational configuration
  interaction using generalized curvilinear coordinates and self-consistent
  basis. J. Chem. Phys. 2010, 133, 094103
[Seidler et al.(2009)Seidler, Matito, and
  Christiansen]Seidler2009_ThreeModeCoupling-VCC
Seidler, P.; Matito, E.; Christiansen, O. Vibrational coupled cluster theory
  with full two-mode and approximate three-mode couplings: The VCC[2pt3] model.
  J. Chem. Phys. 2009, 131, 034115
[Cross et al.(2019)Cross, Bishop, Sheldon, Nation, and
  Gambetta]cross2019
Cross, A. W.; Bishop, L. S.; Sheldon, S.; Nation, P. D.; Gambetta, J. M.
  Validating quantum computers using randomized model circuits. Phys.
  Rev. A 2019, 100, 032328
[Godtliebsen et al.(2013)Godtliebsen, Thomsen, and
  Christiansen]Christianen2013_TensorDecomposition-VCC
Godtliebsen, I. H.; Thomsen, B.; Christiansen, O. Tensor Decomposition and
  Vibrational Coupled Cluster Theory. J. Phys. Chem. A 2013,
  117, 7267–7279
[Madsen et al.(2018)Madsen, Godtliebsen, Losilla, and
  Christiansen]Christiansen2018_VCC-TensorFactorized
Madsen, N. K.; Godtliebsen, I. H.; Losilla, S. A.; Christiansen, O.
  Tensor-decomposed vibrational coupled-cluster theory: Enabling large-scale,
  highly accurate vibrational-structure calculations. J. Chem. Phys.
  2018, 148, 024103
[Motta et al.(2018)Motta, Ye, McClean, Li, Minnich, Babbush,
  and Chan]Motta2018_LowRank
Motta, M.; Ye, E.; McClean, J. R.; Li, Z.; Minnich, A. J.; Babbush, R.;
  Chan, G. K.-L. Low rank representations for quantum simulation of electronic
  structure. arXiv preprint arXiv:1808.02625 2018, 
[Jacob and Reiher(2009)Jacob, and Reiher]Reiher2009_LocalModes
Jacob, C. R.; Reiher, M. Localizing normal modes in large molecules. J.
  Chem. Phys. 2009, 130, 084106
[Thomsen et al.(2014)Thomsen, Yagi, and
  Christiansen]Thomsen2014_OptimizedCoordinates-VCC
Thomsen, B.; Yagi, K.; Christiansen, O. Optimized coordinates in vibrational
  coupled cluster calculations. J. Chem. Phys 2014, 140,
  154102
[Macridin et al.(2018)Macridin, Spentzouris, Amundson, and
  Harnik]Macridin2019_PolaronModel
Macridin, A.; Spentzouris, P.; Amundson, J.; Harnik, R. Digital quantum
  computation of fermion-boson interacting systems. Phys. Rev. A
  2018, 98, 042312
[Paolo et al.(2019)Paolo, Barkoutsos, Tavernelli, and
  Blais]DiPaolo2019_VQE-LightMatter
Paolo, A. D.; Barkoutsos, P. K.; Tavernelli, I.; Blais, A. Variational Quantum
  Simulation of Ultrastrong Light-Matter Coupling. arXiv preprint
  arXiv:1909.08640 2019, 
[Bubin et al.(2013)Bubin, Pavanello, Tung, Sharkey, and
  Adamowicz]Bubin2013_Review
Bubin, S.; Pavanello, M.; Tung, W.-C.; Sharkey, K. L.; Adamowicz, L.
  Born–Oppenheimer and Non-Born–Oppenheimer, Atomic and Molecular
  Calculations with Explicitly Correlated Gaussians. Chem. Rev.
  2013, 113, 36–79
[Mátyus and Reiher(2012)Mátyus, and Reiher]Reiher2012_ECG
Mátyus, E.; Reiher, M. Molecular structure calculations: A unified quantum
  mechanical description of electrons and nuclei using explicitly correlated
  Gaussian functions and the global vector representation. J. Chem.
  Phys. 2012, 137, 024104
[Veis et al.(2016)Veis, Višňák, Nishizawa, Nakai, and
  Pittner]Veis2016_QC-PreBO
Veis, L.; Višňák, J.; Nishizawa, H.; Nakai, H.; Pittner, J. Quantum
  chemistry beyond Born–Oppenheimer approximation on a quantum computer: A
  simulated phase estimation study. Int. J. Quantum Chem. 2016,
  116, 1328–1336
[Li and Benjamin(2017)Li, and Benjamin]Benjamin2017_VQE-Dynamics
Li, Y.; Benjamin, S. C. Efficient Variational Quantum Simulator Incorporating
  Active Error Minimization. Phys. Rev. X 2017, 7,
  021050



