Digital quantum simulation of molecular vibrations

Abstract

Molecular vibrations underpin important phenomena such as spectral properties, energy transfer, and molecular bonding. However, obtaining a detailed understanding of the vibrational structure of even small molecules is computationally expensive. While several algorithms exist for efficiently solving the electronic structure problem on a quantum computer, there has been comparatively little attention devoted to solving the vibrational structure problem with quantum hardware. In this work, we discuss the use of quantum algorithms for investigating both the static and dynamic vibrational properties of molecules. We introduce a physically motivated unitary vibrational coupled cluster ansatz, which also makes our method accessible to noisy, near-term quantum hardware. We numerically test our proposals for the water and sulfur dioxide molecules.

---

I. Introduction

Simulating many-body physical systems enables us to study chemicals and materials without fabricating them, saving both time and resources. The most accurate simulations require a full quantum mechanical treatment—which is exponentially costly for classical computers. While many approximations have been developed to solve this problem, they are often not sufficiently accurate.¹ One possible route to more accurate simulations is to use quantum computers. Quantum computation can enable us to solve certain problems asymptotically more quickly than with a ‘classical’ computer.^2–4 While the quantum computers that we currently possess are small and error prone, it is hoped that we will one day be able to construct a universal, fault-tolerant quantum computer—widely expected to be capable of outperforming its classical counterparts on certain tasks. One example of such tasks is simulating quantum systems on quantum computers.^5–7 In particular, simulating chemical systems, such as molecules,⁸ has received significant attention.^84,85 This may stem from the commercial benefits of being able to investigate and design such systems in silico.⁹ The development of quantum computational chemistry has arguably echoed its classical counterpart. In both fields, the majority of investigations have focused on the electronic structure of molecules.¹⁰ This has resulted in a wealth of well established methods for solving problems of electrons. However, methods concerned with the nuclear degrees of freedom are comparatively less well established. Understanding vibrations is critical for obtaining the most accurate models of real physical systems.¹⁰ Unfortunately, the most detailed classical simulations of vibrations are limited to small molecules, consisting of a few atoms.¹¹ While approximations can be used to treat larger systems, these tend to be less accurate than experiments.¹² Although recently proposed analog quantum algorithms^13–19 are capable of simulating molecular vibrations using resources which scale polynomially with the size of the molecule, the long term scalability of these approaches has yet to be established.

In this work, we discuss a general method for efficiently simulating molecular vibrations on a universal quantum computer. Our method targets the eigenfunctions of a vibrational Hamiltonian with potential terms beyond quadratic order (‘anharmonic potentials’). These wavefunctions can then be used to efficiently calculate properties of interest, such as absorption spectra at finite temperatures. We can also use our method to perform simulations of vibrational dynamics, enabling the investigation of properties such as vibrational relaxation.

II. Vibrations

A consequence of quantum mechanics is that molecules are never at rest, possessing at least the vibrational zero-point energy correction to the ground-state energy.^20–22 As a result, vibrations affect all chemical calculations, to a greater or lesser extent. They are important in both time dependent and independent contexts. From a dynamics perspective, vibrational structure affects high frequency time-resolved laser experiments,²³ reaction dynamics,^24–26 and transport.^27,28 In a static context, vibrations underpin spectral calculations, such as: infrared and Raman spectroscopy²⁹ and fluorescence.³⁰ These calculations determine the performance of solar cells^31,32 and industrial dyes,^33,34 as well as the susceptibility of molecules to photodamage.^13,35

Despite their importance for accurate results, studying vibrations has proven difficult. There are several possible routes to obtaining an accurate description of vibrational behaviour. Real-space, grid based methods, which treat the electronic and nuclear degrees of freedom on an equal footing, are limited to systems of a few particles. While algorithms to efficiently solve this problem on a universal quantum computer exist,^36,37 it will take many years to develop a quantum computer with the required number of qubits.³⁸ Alternatively, one may separate the electronic and nuclear degrees of freedom. We can solve for the electronic energy levels of the system as a function of the nuclear positions, which enables us to map out potential energy surfaces for the system. A number of approximate classical methods have been developed to solve this problem,^1,39 as well as several quantum algorithms.^8,40–42 These electronic potential surfaces can then be viewed as the nuclear potential, determining the vibrational energy levels. This is known as the vibrational structure problem. The accuracy of the nuclear potential is determined by the accuracy of the electronic structure calculation, as well as the number of points obtained for the potential energy surface. Once this potential has been obtained, a number of classical methods can be used for solving both the time dependent and independent Schrödinger equations.

The most simple method uses the ‘harmonic approximation’. This treats the nuclear potential in the vicinity of the equilibrium geometry as a harmonic oscillator potential, resulting in energy eigenstates which are harmonic oscillator eigenfunctions.

Alternatively, one may consider higher order expansions of the nuclear potential, resulting in more accurate calculations.⁴³ One common route towards obtaining the nuclear potential is to first carry out many electronic structure calculations on the system, in the vicinity of the minimum energy configuration. Each of these electronic structure calculations is approximate, and so the cost of each one scales polynomially with the system size. However, if one proceeds to obtain the nuclear potential using this simple grid based method, then a number of grid points scaling exponentially with the number of modes is required.⁴⁴ In practice one can often instead construct an approximate nuclear potential by considering a reduced number of mode couplings, or using interpolation, or using adaptive methods. A review of these, and other state-of-the-art methods can be found in ref. 44. The requirement to first perform multiple electronic structure calculations to obtain the anharmonic nuclear potential makes calculating vibrational energy levels expensive,⁴⁵ even if only mean-field vibrational calculations are then performed. If the correlation between different vibrational modes is included in the calculation, then the simulation becomes even more expensive. While most of the existing classical vibrational simulation methods scale polynomially with the number of modes in the system (e.g. vibrational self-consistent field methods,¹² or vibrational coupled cluster theory⁴⁶), and are sufficiently accurate for some systems, they only provide approximations to the true full configuration interaction vibrational wavefunction, which can be exponentially costly to obtain. A similar hierarchy of accuracy also exists for dynamics simulations.

The computational difficulties described above make accurate vibrational calculations on large systems very challenging for classical computers. To overcome these challenges, quantum solutions have been suggested for the vibrational structure problem.^13–19,47 To date, the majority of suggestions have focused on analog quantum simulation of vibrations. In analog simulations, the simulator emulates a specific system of interest, but cannot in general be programmed to perform simulations of other, different systems. Huh et al. proposed using boson sampling circuits to determine the absorption spectra of molecules.¹³ These boson sampling circuits consist of photons passing through an optical network. This initial proposal relied on the harmonic oscillator approximation at zero temperature, but does take into account bosonic mode mixing due to nuclear structural changes that result from electronic excitation. This method has since been experimentally demonstrated,^15,16 and extended to finite temperature spectra.^14,19 The main limitation of these simulations is the use of the harmonic oscillator approximation for the vibrational wavefunction. It is in general difficult to engineer ground states of anharmonic Hamiltonians using an optical network, as non-linear operations, such as squeezing, are required. Optical networks have also been used for simulating vibrational dynamics.¹⁷ These simulations investigated vibrational transport, adaptive feedback control, and anharmonic effects.

The aforementioned schemes make use of the analogy between the vibrational energy levels in molecules in the Harmonic oscillator approximation, and the bosonic energy levels accessible to photons and ions. One advantage of this is that the bosonic modes are in principle able to store an arbitrary number of excitations. As these analog simulators are relatively simple to construct (when compared to a universal, fault-tolerant quantum computer), they will likely prove useful for small calculations in the near-term. However, it is not yet known how to suppress errors to an arbitrarily low rate in analog simulators. As a result, if we are to simulate the vibrational behaviour of larger quantum systems, we will likely require error corrected universal quantum computers. This motivates our work on methods for vibrational simulation on universal quantum computers.

The rest of this paper is organised as follows. In Section III, we introduce the vibrational structure problem for molecules and show how this problem can be mapped onto a quantum computer. In Section IV, we show how to solve both static and dynamic problems of molecular vibrations. Finally, in Section V, we present the results of numerical simulations of the H₂O and SO₂ molecules.

III. Encoding

A. Vibrational Hamiltonian

Under the Born–Oppenheimer approximation, nuclear variables are treated as parameters in the electronic structure problem and are restored as quantum nuclear variables at the level of the full problem. In the following, we neglect the rotational degrees of freedom from negligible rotational–vibrational couplings for rigid molecules (rigid rotator approximation). After diagonalising the electronic Hamiltonian and neglecting nonadiabatic couplings, the molecular Hamiltonian becomeswhere |ψ_s〉_e are the electronic energy eigenstates. The effective nuclear Hamiltonian H_s iswhere q = (q₁, q₂…) are nuclear coordinates, p = −i∂/∂q are nuclear momenta with ħ = 1, and V_s(q) is the effective nuclear potential. This potential is determined by the corresponding electronic potential energy surface of |ψ_s〉_e. As described in Appendix A, we work in mass-weighted normal coordinates and decouple the rotational and vibrational modes. The potential V_s(q) can be approximated aswhere ω_i is the harmonic frequency of the i^th vibrational normal mode. Thus, the nuclear Hamiltonian H_s can be approximated by a sum of independent harmonic oscillators,with a^†_i and a_i being the creation and annihilation operators of the i^th harmonic oscillator. This is the commonly used ‘harmonic approximation’. Even for accurate potentials and rigid molecules, the harmonic approximation is less accurate than modern spectroscopic techniques.¹² This approximation becomes inadequate for large and ‘floppy’ molecules.¹² Improved results can be obtained by including anharmonic effects which requires information of higher order potential terms in the Hamiltonian.⁴⁴ For example, we can expand the potential aswhere M is the number of modes, k_i1,i2,…,ij are the coefficients of the expansion q_i1q_i2…q_ij, and the harmonic frequencies are In general, the eigenstates of these Hamiltonians are entangled states, when working in a basis of harmonic oscillator eigenstates. Consequently, solving the higher order vibrational Hamiltonian is a hard problem for classical computers.

In contrast, we show below that it is possible to efficiently encode the k^th order nuclear Hamiltonian into a Hamiltonian acting on qubits. We can then use quantum algorithms to efficiently calculate the static and dynamic properties of the nuclear Hamiltonian.

B. Mapping to qubits

We first discuss mapping the molecular Hamiltonian into qubits. We work in the basis of harmonic oscillator eigenstates, as these can be easily mapped to qubits. The direct mapping presented below was originally suggested in the context of simulating general bosonic systems by Somma et al.⁴⁸ It has been used recently in the context of quantum simulation of nuclear physics to investigate the binding energy of a deuteron nucleus.⁴⁹ The compact mapping discussed below was proposed by Veis et al., in the context of using quantum computers to simulate ‘nuclear orbital plus molecular orbital (NOMO)’ theory, which uses Gaussian orbitals for the nuclei, and treats them on an equal footing to the electrons.⁵⁰ This differs from our work, which separates the nuclear and electronic degrees of freedom, and predominantly considers a harmonic oscillator basis for the vibrational modes. While this tailors our method for vibrational problems, it means we are limited to solving problems for which the Born–Oppenheimer approximation is valid, unlike ref. 50.

Focusing first on one harmonic oscillator, ĥ = ωa^†a, we consider the truncated eigenstates with the lowest d energies, |s〉 with s = 0,1,…,d − 1. A direct mapping of the space {|s〉} is to encode it with d qubits as

|s〉 = ⊗^s−1_j=0|0〉_j|1〉_s⊗^d−1_j=s+1|0〉_j, (6)

with creation operator

The annihilation operator can be obtained by taking the Hermitian conjugate of a^†. As an alternative to the direct mapping, we can use a compact mapping, using K = [log d]′ qubits,

|s〉 = |b_K−1〉|b_K−2〉…|b₀〉, (8)

with binary representation s = b_K−12^K−1 + b_K−22^K−2 + …b₀2⁰. The representation of the creation operator is

These binary projectors can then be mapped to Pauli operators;

when decomposing a^† and a into local Pauli matrices, there are O(d) and O(d²) terms for the direct and compact mapping, respectively. In Fig. 1, we show the number of qubits required to describe the vibrational Hamiltonians of several molecules, for both mappings.

Fig. 1 Number of qubits required for the direct and compact mappings with d = 4 energy levels for each mode.

As p and q can both be represented by a linear combination of creation and annihilation operators, we can thus map the nuclear vibrational Hamiltonian to a qubit Hamiltonian. If the molecule has n atoms, it has M = 3n − 6 vibrational modes for a nonlinear molecule, and M = 3n − 5 for a linear molecule. The vibrational wavefunction can then be represented with Md (direct mapping) or M log(d) (compact mapping) qubits. This can be contrasted with the exponentially scaling classical memory required to store the wavefunction. If the potential is expanded to k^th order (with M ≫ k), the Hamiltonian contains O(M^kd^k) (direct) or O(M^kd^2k) (compact) terms. These terms are strings of local Pauli matrices. In this work, we take d to be a small constant. This approximation constrains us to the low energy subspace of the Hamiltonian, which should be valid for calculations of ground and low-lying excited states. The applicability of this approximation to the simulation of dynamics is discussed in Section VI. We set k = 4 to investigate the Hamiltonian to quartic order. The resulting Hamiltonian has O(M⁴) terms.

IV. Simulating molecular vibrations

Once the vibrational modes have been mapped to qubits, we can use quantum algorithms to obtain the static and dynamic properties of the system. We can write the qubit Hamiltonian as where λ_i are coefficients determining the strength of each term in the Hamiltonian.

A. Vibrational energy levels

An important, but classically difficult problem, is to obtain accurate energy levels for the vibrational Hamiltonian. The spectrum of the vibrational Hamiltonian provides corrections to the electronic eigenstates used to predict reaction rates.²² Moreover, we will show how these energy levels can be used to calculate the absorption spectrum of molecules in Section IV.B. Of particular interest are the lowest lying energy levels at low temperature. Using a universal quantum computer, we can first prepare an initial state that has a large overlap with the ground state of the vibrational Hamiltonian. We can then use the phase estimation algorithm^7,51 to probabilistically obtain the ground state and ground state energy. A possible initial ground state is the lowest energy product state of the harmonic oscillator basis states, |ψ₀〉 = ⊗_m|s_m〉_m. However, we note that the overlap between this state and the true ground state may decrease exponentially with the size of the molecule. This so-called ‘orthogonality catastrophe’ has been discussed previously in the context of electronic structure calculations on a quantum computer.^52,53 As a result, for large systems it may be more efficient to use an initial state obtained from a classical vibrational self-consistent field (VSCF) calculation. VSCF is the vibrational analogue of the Hartree–Fock method in electronic structure theory. VSCF optimises the basis functions to minimise the energy of the Hamiltonian with a product state.

Another route to a state with a large overlap with the ground state, is to prepare the VSCF state, and then adiabatically evolve under a Hamiltonian that changes slowly from the VSCF Hamiltonian, to the full vibrational Hamiltonian H. This approach has received significant attention within quantum computing approaches to the electronic structure problem, since it was first proposed in the context of quantum computational chemistry in ref. 8. However, both adiabatic state preparation and phase estimation typically require long circuits, with a large number of gates. As a result, quantum error correction is required to suppress the effect of device imperfections. It is therefore helpful to introduce variational methods, which may make these calculations feasible for near-term, non-error corrected quantum computers. Variational methods replace the long gate sequences required by phase estimation with a polynomial number of shorter circuits.^42,54 This dramatically reduces the coherence time required. As a result, quantum error correction may not be required, if the error rate is sufficiently low, in the context of the number of gates required. The circuits used consist of a number of parametrised gates which seek to create an accurate approximation of the desired state. The parameters are updated using a classical feedback loop, in order to produce better approximations of the desired state. The circuit used is known as the ‘ansatz’ circuit.

Inspired by classical methods for the vibrational structure problem, we introduce the unitary vibrational coupled cluster (UVCC) ansatz. This is a unitary analogue of the VCC ansatz introduced in ref. 12 and 46. We note that a similar pairing exists for the electronic structure problem, where the unitary coupled cluster (UCC) ansatz^55,56 has been suggested as a quantum version of the classical coupled cluster method. The UVCC ansatz is given by

|Ψ([small theta, Greek, vector])〉 = exp(T̂ − T̂^†)|Ψ₀〉, (11)

where the initial state |Ψ₀〉 can be either the ground state |ψ₀〉 of the harmonic oscillators or the VSCF state |Ψ_VSCF〉, T̂ is the sum of molecular excitation operators truncated at a specified excitation rank, and [small theta, Greek, vector] are the parameters defined below. Similar to the unitary coupled cluster ansatz in electronic structure problems, the single and double excitation operators are

T̂ = T̂₁ + T̂₂ +…, (12)

withHere, we omit the subscript of modes for simplicity. θ_sm,tm and θ_sm,tm,pm,qn are real parameters, and [small theta, Greek, vector] = {θ_sm,tm,θ_sm,tm,pm,qn}. The T̂ operators can be mapped to qubit operators via either the direct or compact mapping.

The UVCC ansatz seeks to create a good approximation to the true ground state by considering excitations above a reference state. We note that the classical VCC ansatz is not a unitary operator. Correspondingly, the method is not variational, meaning that energies are not bounded from below. Moreover, we expect that the UVCC ansatz will deal better with problems of strong static correlation than the VCC ansatz, as the former can easily be used with multi-reference states. This echoes the way in which the UCC ansatz can be used with multi-reference states,⁵⁶ while it is typically more difficult when using the canonical CC method.¹

Once we have obtained the energy levels of the vibrational Hamiltonian using the methods discussed above, we can calculate the infrared and Raman frequencies, using the difference between the excited and ground-state energies.⁵⁷

It is often also the case that one is interested in the properties of a system in thermal equilibrium, rather than a specific eigenstate. We can also use established quantum algorithms with the Hamiltonians described above to construct these thermal states. On error corrected quantum computers, we can use the heuristic algorithms presented in ref. 58 and 59 to construct these thermal states. Alternatively, we can use near-term devices to implement hybrid algorithms for imaginary time evolution.^60–62

B. Franck–Condon factors

In addition to focusing on the eigenstates or thermal states of a single vibrational Hamiltonian, we can also consider vibronic (vibrational and electronic) transitions between the vibrational levels resulting from different electronic potential energy surfaces. Consider two electronic states, |i〉_e and |f〉_e. The molecular Hamiltonian is

H_mol = |i〉〈i|_e⊗H_i + |f〉〈f|_e⊗H_f, (14)

where H_i and H_f are vibrational Hamiltonians, with energy eigenstates |ψⁱ_vib〉 and |ψ^f_vib〉, respectively. Using Fermi's Golden rule, the probability of a photon-induced transition between two wavefunctions |ψⁱ〉 = |i〉⊗|ψⁱ_vib〉 and |ψ^f〉 = |f〉⊗|ψ^f_vib〉 is proportional to the square of the transition dipole moment P² = |〈ψ^f|[small mu, Greek, circumflex]|ψⁱ〉|², using first order time-dependent perturbation theory. Within the Condon approximation, [small mu, Greek, circumflex] = [small mu, Greek, circumflex]_e + [small mu, Greek, circumflex]_N, the transition probability becomes proportional to P² = |〈ψ^f_vib|ψⁱ_vib〉|²·|〈f|[small mu, Greek, circumflex]_e|i〉|². Here |〈ψ^f_vib|ψⁱ_vib〉|² are referred to as Franck–Condon integrals. Without the Condon approximation, the Franck–Condon integrals become |〈ψ^f_vib|[small mu, Greek, circumflex](q)|ψⁱ_vib〉|², with [small mu, Greek, circumflex](q) = |〈f|[small mu, Greek, circumflex]|i〉|².

In practice, |ψⁱ_vib〉 and |ψ^f_vib〉 are eigenstates of Hamiltonians with different harmonic oscillator normal modes q^f and qⁱ. These modes are related by the Duschinsky transform q^f = Uqⁱ + d.^45,63 According to the Doktorov unitary representation of the Duschinsky transform, the harmonic oscillator eigenstates are related by^13,14,64

|s_f〉 = Û_Dok|s_i〉 (15)

where |s_i〉 and |s_f〉 are harmonic oscillator eigenstates in the initial and final coordinates qⁱ and q^f, respectively. The Doktorov unitary can be decomposed into a product of unitary operators Û_Dok = Û_tÛ^†_s′Û_sÛ_r, which depend on the displacement vector d, the rotation matrix U, and matrices of the frequencies of the harmonic oscillators and The definitions of the unitary operators are shown in Appendix C.

If |Ψⁱ_vib〉 and |Ψ^f_vib〉 are the qubit wavefunctions resulting from diagonalisation of H_i and H_f using a quantum computer, they will be obtained in different normal mode bases |s_i〉 and |s_f〉, respectively. We cannot directly calculate the Franck–Condon integrals using |〈Ψ^f_vib|Ψⁱ_vib〉|², as this does not take into account the different bases. Instead, we must implement the Doktorov unitary to get the Franck–Condon integrals

|〈ψ^f_vib|ψⁱ_vib〉|² = |〈Ψ^f_vib|Û_Dok|Ψⁱ_vib〉|². (16)

The Franck–Condon integrals without the Condon approximation can be efficiently calculated via

|〈ψ^f_vib|[small mu, Greek, circumflex](q)|ψⁱ_vib〉² = |〈Ψ^f_vib|[small mu, Greek, circumflex](q^f)Û_Dok|Ψⁱ_vib〉|². (17)

They can be both efficiently calculated with the generalised SWAP-test circuit.⁶⁵

Alternatively, we can obtain the Franck–Condon integrals without realising the Doktorov transform. The qubit states |Ψⁱ_vib〉 and |Ψ^f_vib〉 are obtained from H_i(qⁱ) and H_f(q^f) with normal mode coordinates qⁱ and q^f, respectively. Instead, we can focus on one set of normal mode coordinates qⁱ and represent the Hamiltonian H_f in qⁱ, H′_f(qⁱ). By solving the energy eigenstates of H′_f(qⁱ), we can directly get and calculate the Franck–Condon integrals without realising the Doktorov transform. However, as the Hamiltonian H′_f(qⁱ) is not encoded in the correct normal mode basis, the ground state of the harmonic oscillators or the VSCF state |Ψ_VSCF〉 may not be an ideal initial state to start with. However, this effect may be negligible if the overlap between |Ψⁱ_vib〉 and is suitably large. In this case, the initial state |Ψ₀〉 for |Ψⁱ_vib〉 should also be an ideal initial state for The aforementioned transformation can be implemented by transforming the normal mode coordinates qⁱ as described in ref. 45.

C. Vibrational dynamics

In this section, we consider methods to investigate the dynamic properties of vibrational Hamiltonians. Vibrational dynamics underpin phenomena including energy and electron transport^27,28 and chemical reactions.^24–26 Dynamical behaviour can be studied by transforming to a single-mode basis of spatially localised vibrational modes, as described in ref. 17. The spatially localised vibrational modes a^L_i, are related to the normal modes a_ivia a basis transformationwith real unitary matrix U_i,j. We can obtain the corresponding localised Hamiltonian H^L, using the transformation of the normal coordinates and momentaGiven an initial state of the localised vibrations, the dynamics can be simulated by applying the time evolution operator e^−iHLt. This can be achieved in a number of ways, using different Hamiltonian simulation algorithms, including: Trotterization (also referred to as product formulae),^5,66 the Taylor series method,^67–69 and qubitization^70,71 in conjunction with quantum signal processing.^72,73 The product formula method is the most simple to realise. If H^L can be decomposed as the time evolution operator e^−iHLt can be realised using a product formula,where N is chosen to be sufficiently large to suppress the error in the approximation.

Alternatively, the vibrational dynamics can be realised using a recently proposed variational algorithm.⁷⁴ One could use either a UVCC ansatz, or a Trotterized ansatz.^75,76

V. Numerical simulations

In this section, we demonstrate how the techniques described above can be used to calculate the vibrational energy levels of small molecules. We focus on the polyatomic molecules H₂O and SO₂, which both have three vibrational modes. The coefficients were computed at MP2/aug-cc-pVTZ level of theory using Gaussian09 software.⁸³ in Appendix D. We consider the cases with two and four energy levels for each of the harmonic oscillator modes, yielding Hamiltonians acting on 6 and 12 qubits for the direct mapping, and 3 and 6 qubits for the compact mapping. We used the compact mapping in our numerical simulations, as it requires fewer qubits. There are 216 and 165 terms in the Hamiltonian for H₂O and SO₂, respectively.

Table 1 Coefficients of the potential energy surface of H₂O and SO₂. The coefficients are in atomic units, where the unit of length is a₀ = 1 Bohr (0.529167 × 10⁻¹⁰ m), the unit of mass is the electron mass m_e, and the unit of energy is 1 Hartree (1 Hartree = e²/4πε₀a₀ = 27.2113 eV)

k	H2O	SO2
k 1,1	0.275240 × 10^−4	0.252559 × 10^−5
k 2,2	0.151618 × 10^−3	0.125410 × 10^−4
K 3,3	0.161766 × 10^−3	0.176908 × 10^−4
K 1,1,1	0.121631 × 10^−6	0.316 646 × 10^−8
K 1,1,2	0.698476 × 10^−6	0.575325 × 10^−8
K 1,2,2	−0.266427 × 10^−6	0.197771 × 10^−7
k 2,2,2	−0.312538 × 10^−5	−0.668689 × 10^−7
K 1,3,3	−0.915428 × 10^−6	−0.370850 × 10^−9
k 2,3,3	−0.964649 × 10^−5	−0.284244 × 10^−6
K 1,1,1,1	−0.463748 × 10^−9	0.330842 × 10^−11
K 1,1,2,2	−0.449480 × 10^−7	−0.172869 × 10^−9
K 1,2,2,2	0.957558 × 10^−8	−0.215928 × 10^−9
k 2,2,2,2	0.433267 × 10^−7	0.225400 × 10^−9
K 1,1,3,3	−0.555026 × 10^−7	−0.356155 × 10^−9
K 1,2,3,3	0.563566 × 10^−7	−0.128135 × 10^−9
k 2,2,3,3	0.269239 × 10^−6	−0.220168 × 10^−8
K 3,3,3,3	0.462143 × 10^−7	0.458046 × 10^−9
k 2,3,3,3	0	−0.720760 × 10^−11

---

We first calculate the energy levels under the harmonic approximation. We compare this to the energy levels obtained with a fourth order expansion of the potential. The results for H₂O are shown in Fig. 2. We can see that although the ground state can be well approximated by the harmonic oscillators, the excited states deviate from the harmonic oscillators at higher energy levels. The results for SO₂ can be found in the Appendix. These calculations highlight the importance of anharmonic terms in the potential for even small molecules.

Fig. 2 Vibrational spectra of H₂O with two and four energy levels for each mode. The solid lines are the energy levels of the harmonic oscillator eigenstates and the dashed lines are the vibrational spectra of the Hamiltonian with a fourth order expansion of the potential.

Next, we implemented the UVCC ansatz to obtain the vibrational energy levels of H₂O using the variational quantum eigensolver.⁴² For simplicity, we considered two energy levels for each mode. To implement the UVCC ansatz, we first calculate the imaginary part of T̂ and encode it into a linear combination of local Pauli terms, i.e.,

Then, as for the UCC ansatz, we realise exp(T̂ − T̂^†) by a first order Trotterisation via For example, the UVCC ansatz of H₂O with two energy levels can be prepared by the circuit in Fig. 3.

Fig. 3 The UVCC ansatz of three modes each with two energy levels. There are nine gates with six parameters (joined gates share the same parameters). The single qubit gate on the i^th qubit is and the two qubit gate on the i^th and j^th is

Using the UVCC ansatz, we can obtain the vibrational ground state with a variational procedure. As the ground state is close to the initial state |0〉^⊗3, we start with parameters slightly perturbed from zero. We then use gradient descent to find the minimum energy of the system. The results are shown in Fig. 4.

Fig. 4 Solving the vibrational ground state of H₂O with the UVCC ansatz. Here, we consider two energy levels for each mode.

VI. Discussion

In this work, we have extended many of the techniques developed for quantum simulation to the problem of simulating vibrations. Understanding vibrations is an important problem for accurately modelling chemical systems, yet one that is difficult to solve classically.

We have discussed ways to map between vibrational modes and qubit states. This is only possible when we consider a restricted number of harmonic oscillator energy levels. This approximation is appropriate when investigating the low energy properties of the Hamiltonian, such as the ground state energy. However, it may not always be possible to use this approximation when considering time evolution, as this requires the exponentiation of our truncated Hamiltonian. One possible route to overcome this challenge is to repeat simulations which consider an increasing number of energy levels, and then to extrapolate to the infinite energy level result. This technique was used in a similar context to this work in ref. 49.

Once the vibrational Hamiltonian has been mapped to a qubit Hamiltonian, much of the existing machinery for quantum simulation can be applied. Static properties, such as energy levels, can be calculated using phase estimation or variational approaches. To aid variational state preparation, we proposed a unitary version of the powerful VCC method used in classical vibrational simulations. The resulting energy eigenstates can be used as an input for SWAP-test circuits. These calculate the Franck–Condon factors for the molecules, which are related to the absorption spectra. Alternatively, one may investigate dynamic properties, using methods for Hamiltonian simulation to time evolve a specified state.

Compared with analog algorithms,^13,14,16,17 our method can easily take into account anharmonic terms in the nuclear Hamiltonian. Moreover, it could be used to simulate large systems by protecting the quantum computer with error correction. As our technique is tailored for simulating vibrational states, it makes it simple to investigate interesting vibrational properties, such as the Franck–Condon factors. However, as our approach uses the Born–Oppenheimer approximation, it is not suitable for all problems in chemistry, such as problems including relativistic effects⁷⁷ or conical intersections.^78–80 Future work will address whether restricting our vibrational modes to low-lying energy levels poses a significant challenge for problems of practical interest.

Appendix A: encoding vibrational Hamiltonians into qubits

The molecular Hamiltonian in atomic units iswhere M_I, R_I, and Z_I are the mass, position, and charge of nuclei I, and r_i is the position of electron i. Given the location of the nucleus, the electronic Hamiltonian isand the total Hamiltonian is

Under the Born–Oppenheimer approximation, we assume the electrons and nuclei are in a product state,

|ψ〉 = |ψ〉_n|ψ〉_e. (A4)

To get the ground state of the Hamiltonian, one can thus separately minimise over |ψ〉_n and |ψ〉_e,As only H_e(R_I) depends on |ψ〉_e, the minimisation over |ψ〉_e is equivalent to finding the ground state of H_e(R_I). Denotethen the ground state of H_mol can be found by solving the ground state of H₀,

In general, considering a spectral decomposition of the molecular Hamiltonian is

Here, |ψ_s〉_e are eigenstates of the electronic Hamiltonian andand

Finding the spectra of H_e(R_I) is called the electronic structure problem, which can be efficiently solved using a quantum computer.⁸ One approach is to consider a subspace that the ground state lies in and transform the Hamiltonian H_e(R_I) into the second quantised formulation, with a basis determined by the subspace. As electrons are fermions, the obtained Hamiltonian is a fermionic Hamiltonian. By using the standard encoding methods, such as Jordan–Wigner and Bravyi–Kitaev,⁸² the fermionic Hamiltonian is converted into a qubit Hamiltonian, whose spectra can be efficiently computed.

Focusing on the ground state of the electronic structure Hamiltonian, we show how to redefine H₀ in the mass-weighted basis and how to encode it with qubits. Denote then one can obtain the mass-weighted normal coordinates q_i by minimising the coupling between the rotational and vibrational degrees of freedom and diagonalising the Hessian matrix,

In the mass-weighted basis, the potential can be expanded via a Taylor series truncated at fourth orderand the total Hamiltonian becomes

If we consider the higher order terms as a perturbation, one can get the normal modes by solving the harmonic oscillator

We denote the eigenbasis for ĥ_i as , then the nuclear wave function can be represented byThe Hamiltonian under the normal mode basis becomes,Suppose the basis is truncated to the lowest d energy levels, the space of H_{s¹s²…s^M,t¹t²…t^M} is equivalent to the space of M d-level systems, or equivalently M log₂ d qubits.

Appendix B: variational quantum simulation with the unitary vibrational coupled cluster ansatz

We can make use of variational methods to find the low energy spectra of the vibrational Hamiltonian. Inspired by classical computational chemistry, we introduce the unitary vibrational coupled cluster (UVCC) ansatz

|VCC〉 = exp(T̂ − T̂^†)|Φ₀〉 (B1)

where the reference state |Φ₀〉 is a properly chosen initial state, and T is the sum of molecular excitation operators truncated at a specified excitation rank (often single and double excitations)⁵⁶

T̂ = T̂₁ + T̂₂ +…

with

The initial state can be the product of the ground-state of each mode

|Φ₀〉 = |ψ¹₀ψ²₀…ψ^M₀〉. (B2)

Or we can also run a vibrational self-consistent field (VSCF) to get the Hartree–Fock initial state

|Φ₀〉 = |ϕ¹ϕ²…ϕ^M〉, (B3)

which is obtained by minimising the energy of the Hamiltonianby solving the self-consistent equation

H_i|ϕ_i〉 = E_i|ϕ_i〉, (B5)

with H_i = 〈ϕ¹…ϕⁱ⁻¹ϕⁱ⁺¹…ϕ^M|H₀|ϕ¹…ϕⁱ⁻¹ϕⁱ⁺¹…ϕ^M〉.

Appendix C: Duschinsky transform

The relation between the initial coordinates q₁ and final coordinates q₂ is

q₁ = Uq₂ + d, (C1)

where U is the Duschinsky rotation matrix and d is the displacement vector. The harmonic oscillator Hamiltonian with unit mass iswith The operators q₁ and p₁ can be represented by the creation and annihilation operators

The transformation for p is p₁ = U^†p₂. The transformation for the creation operators are

where and .

It was shown by Doktorov et al.⁶⁴ that the Duschinsky transform can be implemented using a unitary transform inserted into the overlap integral

〈ν_f|ν_i〉 = 〈ν′_f|Û_Dok|ν_i〉 (C5)

where |v〉 is a harmonic oscillator eigenstate in the initial coordinate q and |v′〉 is a harmonic oscillator eigenstate in the final coordinate q′. The Doktorov unitary can be decomposed into a product of unitaries, which depend on the displacement vector d⃑, the rotation matrix U, and matrices of the eigenenergies of the harmonic oscillator states ; and . It is given by

Û_Dok = Û_tÛ^†_s′Û_sÛ_r (C6)

wherewhere a⃑ = (a₀,…,a_M)^T, andandand

These exponentials could be expanded into local qubit operators using Trotterization. It is important to note that these relations are only valid when the single-mode basis functions are chosen to be harmonic oscillator eigenstates.

Appendix D: numerical simulation

The simulation results of the vibrational energy levels of SO₂ are shown in Fig. 5.

Fig. 5 Vibrational spectra of SO₂ with two and four energy levels for each mode. The solid lines are the energy levels of the harmonic oscillator eigenstates and the dashed lines are the vibrational spectra of the Hamiltonian with a fourth order expansion of the potential.

Note added

After this work was released as a preprint, a notable related work was released online.⁸¹ That paper provides a quantum algorithm for efficiently calculating the entire converged Franck–Condon profile, as opposed to the method we present herein, which calculates a single Franck–Condon factor. However, the techniques we have discussed for finding vibrational ground states and thermal states provide a way to realise the crucial first step of their algorithm. As such, our two works are highly complementary, and can be compared.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge insightful comments from Joonsuk Huh. This work was supported by BP plc and the EPSRC National Quantum Technology Hub in Networked Quantum Information Technology (EP/M013243/1). X. Shan acknowledges the use of the University of Oxford Advanced Research Computing (ARC) facility http://dx.doi.org/10.5281/zenodo.22558.

References

1. T. Helgaker, P. Jorgensen and J. Olsen, Molecular electronic-structure theory, John Wiley & Sons, 2014.
2. R. P. Feynman, Int. J. Theor. Phys., 1982, 21, 467–488.
3. P. W. Shor, 1994 Proceedings, 35th Annual Symposium on Foundations of Computer Science, 1994, pp. 124–134.
4. L. K. Grover, Proceedings of the twenty-eighth annual ACM symposium on Theory of computing, 1996, pp. 212–219.
5. S. Lloyd, Science, 1996, 273, 1073–1078.
6. D. S. Abrams and S. Lloyd, Phys. Rev. Lett., 1997, 79, 2586–2589.
7. D. S. Abrams and S. Lloyd, Phys. Rev. Lett., 1999, 83, 5162–5165.
8. A. Aspuru-Guzik, A. D. Dutoi, P. J. Love and M. Head-Gordon, Science, 2005, 309, 1704–1707.
9. A. Aspuru-Guzik, R. Lindh and M. Reiher, ACS Cent. Sci., 2018, 4, 144–152.
10. O. Christiansen, Phys. Chem. Chem. Phys., 2007, 9, 2942–2953.
11. J. M. Bowman, T. Carrington and H.-D. Meyer, Mol. Phys., 2008, 106, 2145–2182.
12. O. Christiansen, J. Chem. Phys., 2004, 120, 2140–2148.
13. J. Huh, G. G. Guerreschi, B. Peropadre, J. R. McClean and A. Aspuru-Guzik, Nat. Photonics, 2015, 9, 615.
14. J. Huh and M.-H. Yung, Sci. Rep., 2017, 7, 7462.
15. W. R. Clements, J. J. Renema, A. Eckstein, A. A. Valido, A. Lita, T. Gerrits, S. W. Nam, W. S. Kolthammer, J. Huh and I. A. Walmsley, arXiv preprint arXiv:1710.08655, 2017.
16. Y. Shen, Y. Lu, K. Zhang, J. Zhang, S. Zhang, J. Huh and K. Kim, Chem. Sci., 2018, 9, 836–840.
17. C. Sparrow, E. Martín-López, N. Maraviglia, A. Neville, C. Harrold, J. Carolan, Y. N. Joglekar, T. Hashimoto, N. Matsuda and J. L. OBrien, et al. , Nature, 2018, 557, 660.
18. S. Chin and J. Huh, arXiv preprint arXiv:1803.10002, 2018.
19. L. Hu, Y.-C. Ma, Y. Xu, W.-T. Wang, Y.-W. Ma, K. Liu, H.-Y. Wang, Y.-P. Song, M.-H. Yung and L.-Y. Sun, Sci. Bull., 2018, 63, 293–299.
20. R. P. de Tudela, F. J. Aoiz, Y. V. Suleimanov and D. E. Manolopoulos, J. Phys. Chem. Lett., 2012, 3, 493–497.
21. S. Karazhanov, M. Ganchenkova and E. Marstein, Chem. Phys. Lett., 2014, 601, 49–53.
22. A. Gross and M. Scheffler, J. Vac. Sci. Technol., A, 1997, 15, 1624–1629.
23. T. Seideman, Forming Superposition States, in Computational Molecular Spectroscopy, Wiley, Chichester, 2000, pp. 589–624.
24. F. F. Crim, Proc. Natl. Acad. Sci., 2008, 105, 12654–12661.
25. D. Antoniou and S. D. Schwartz, J. Phys. Chem. B, 2001, 105, 5553–5558.
26. D. L. Proctor and H. F. Davis, Proc. Natl. Acad. Sci., 2008, 105, 12673–12677.
27. R. Borrelli, M. D. Donato and A. Peluso, Theor. Chem. Acc., 2006, 117, 957–967.
28. H. Hwang and P. J. Rossky, J. Phys. Chem. A, 2004, 108, 2607–2616.
29. C. Zhu, K. K. Liang, M. Hayashi and S. H. Lin, Chem. Phys., 2009, 358, 137–146.
30. T.-W. Huang, L. Yang, C. Zhu and S. H. Lin, Chem. Phys. Lett., 2012, 541, 110–116.
31. S.-Y. Yue, X. Zhang, G. Qin, J. Yang and M. Hu, Phys. Rev. B, 2016, 94, 115427.
32. L. Debbichi, M. C. Marco de Lucas, J. F. Pierson and P. Krger, J. Phys. Chem. C, 2012, 116, 10232–10237.
33. S. Dhananasekaran, R. Palanivel and S. Pappu, J. Adv. Res., 2016, 7, 113–124.
34. N. Biswas and S. Umapathy, J. Phys. Chem. A, 1997, 101, 5555–5566.
35. K.-W. Choi, J.-H. Lee and S. K. Kim, J. Am. Chem. Soc., 2005, 127, 15674–15675.
36. I. Kassal, S. P. Jordan, P. J. Love, M. Mohseni and A. Aspuru-Guzik, Proc. Natl. Acad. Sci., 2008, 105, 18681–18686.
37. I. D. Kivlichan, N. Wiebe, R. Babbush and A. Aspuru-Guzik, J. Phys. A: Math. Theor., 2017, 50, 305301.
38. N. C. Jones, J. D. Whitfield, P. L. McMahon, M.-H. Yung, R. Van Meter, A. Aspuru-Guzik and Y. Yamamoto, New J. Phys., 2012, 14, 115023.
39. A. Szabo and N. S. Ostlund, Modern quantum chemistry: introduction to advanced electronic structure theory, Courier Corporation, 2012.
40. R. Babbush, D. W. Berry, I. D. Kivlichan, A. Y. Wei, P. J. Love and A. Aspuru-Guzik, New J. Phys., 2016, 18, 033032.
41. R. Babbush, D. W. Berry, Y. R. Sanders, I. D. Kivlichan, A. Scherer, A. Y. Wei, P. J. Love and A. Aspuru-Guzik, Quantum Sci. Technol., 2017, 3, 015006.
42. A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik and J. L. Obrien, Nat. Commun., 2014, 5, 4213.
43. O. Christiansen, Phys. Chem. Chem. Phys., 2007, 9, 2942–2953.
44. O. Christiansen, Phys. Chem. Chem. Phys., 2012, 14, 6672–6687.
45. J. Huh, PhD thesis, Frankfurt institute for advanced studies, Johann Wolfang Goethe Universitat, Frankfurt, 2011.
46. O. Christiansen, J. Chem. Phys., 2004, 120, 2149–2159.
47. S. Joshi, A. Shukla, H. Katiyar, A. Hazra and T. S. Mahesh, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 90, 022303.
48. R. D. Somma, G. Ortiz, E. H. Knill and J. Gubernatis, Proc. SPIE, 2003, 5105.
49. E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean and P. Lougovski, Phys. Rev. Lett., 2018, 120, 210501.
50. L. Veis, J. Vik, H. Nishizawa, H. Nakai and J. Pittner, Int. J. Quantum Chem., 2016, 116, 1328–1336.
51. A. Y. Kitaev, preprint at http://arxiv.org/abs/quant-ph/9511026, 1995.
52. J. R. McClean, R. Babbush, P. J. Love and A. Aspuru-Guzik, J. Phys. Chem. Lett., 2014, 5, 4368–4380.
53. N. M. Tubman, C. Mejuto-Zaera, J. M. Epstein, D. Hait, D. S. Levine, W. Huggins, Z. Jiang, J. R. McClean, R. Babbush, M. Head-Gordon and K. B. Whaley, arXiv:1809.05523, 2018.
54. J. R. McClean, J. Romero, R. Babbush and A. Aspuru-Guzik, New J. Phys., 2016, 18, 023023.
55. M.-H. Yung, J. Casanova, A. Mezzacapo, J. McClean, L. Lamata, A. Aspuru-Guzik and E. Solano, Sci. Rep., 2014, 4, 3589.
56. J. Romero, R. Babbush, J. McClean, C. Hempel, P. Love and A. Aspuru-Guzik, Quantum Sci. Technol., 2018, 4, 014008.
57. E. B. Wilson, J. C. Decius and P. C. Cross, Molecular vibrations: the theory of infrared and Raman vibrational spectra, Courier Corporation, 1980.
58. K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin and F. Verstraete, Nature, 2011, 471, 87.
59. M.-H. Yung and A. Aspuru-Guzik, Proc. Natl. Acad. Sci., 2012, 109, 754–759.
60. S. McArdle, T. Jones, S. Endo, Y. Li, S. Benjamin and X. Yuan, arXiv:1804.03023, 2018.
61. X. Yuan, S. Endo, Q. Zhao, S. Benjamin and Y. Li, 2018, arXiv:1812.08767.
62. M. Motta, C. Sun, A. T. K. Tan, M. J. O. Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandao and G. K.-L. Chan, 2019, arXiv:1901.07653.
63. H. Kupka and P. H. Cribb, J. Chem. Phys., 1986, 85, 1303–1315.
64. E. Doktorov, I. Malkin and V. Man'ko, J. Mol. Spectrosc., 1977, 64, 302–326.
65. L. Cincio, Y. Suba, A. T. Sornborger and P. J. Coles, 2018, arXiv:1803.04114.
66. H. F. Trotter, Proc. Am. Math. Soc., 1959, 10, 545–551.
67. D. W. Berry and A. M. Childs, Quantum Inf. Comput., 2012, 12, 29–62.
68. D. W. Berry, A. M. Childs, R. Cleve, R. Kothari and R. D. Somma, Phys. Rev. Lett., 2015, 114, 090502.
69. D. W. Berry, A. M. Childs and R. Kothari, 2015 IEEE 56th Annual Symposium on Foundations of Computer Science, 2015, pp. 792–809.
70. G. H. Low and I. L. Chuang, arXiv:1610.06546, 2016.
71. G. H. Low, arXiv preprint arXiv:1807.03967, 2018.
72. G. H. Low, T. J. Yoder and I. L. Chuang, Phys. Rev. X, 2016, 6, 041067.
73. G. H. Low and I. L. Chuang, Phys. Rev. Lett., 2017, 118, 010501.
74. Y. Li and S. C. Benjamin, Phys. Rev. X, 2017, 7, 021050.
75. J. O'Gorman and E. T. Campbell, Phys. Rev. A, 2017, 95, 032338.
76. T. Jones and S. C. Benjamin, arXiv:1811.03147, 2018.
77. M. Reiher and A. Wolf, Relativistic quantum chemistry: the fundamental theory of molecular science, John Wiley & Sons, 2014.
78. W. Domcke, D. R. Yarkony and H. Köppel, Conical Intersections, World Scientific, 2004.
79. W. Domcke and D. R. Yarkony, Annu. Rev. Phys. Chem., 2012, 63, 325–352.
80. I. G. Ryabinkin, L. Joubert-Doriol and A. F. Izmaylov, Acc. Chem. Res., 2017, 50, 1785–1793.
81. N. P. D. Sawaya and J. Huh, 2018, arXiv:1812.10495.
82. J. T. Seeley, M. J. Richard and P. J. Love, J. Chem. Phys., 2012, 137, 224109.
83. M. Frisch, http://www.gaussian.com/, 2009.
84. S. McArdle, S. Endo, A. Aspuru-Guzik, S. Benjamin and X. Yuan, 2018, arXiv:1808.10402.
85. Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis and A. Aspuru-Guzik, 2018, arXiv:1812.09976.

---